**Assessment of Landslide-Induced Geomorphological Changes in Hítardalur Valley, Iceland, Using Sentinel-1 and Sentinel-2 Data**

**Zahra Dabiri 1,\*, Daniel Hölbling 1, Lorena Abad 1, Jón Kristinn Helgason 2, Þorsteinn Sæmundsson <sup>3</sup> and Dirk Tiede <sup>1</sup>**


Received: 17 July 2020; Accepted: 20 August 2020; Published: 24 August 2020

**Abstract:** Landslide mapping and analysis are essential aspects of hazard and risk analysis. Landslides can block rivers and create landslide-dammed lakes, which pose a significant risk for downstream areas. In this research, we used an object-based image analysis approach to map geomorphological features and related changes and assess the applicability of Sentinel-1 data for the fast creation of post-event digital elevation models (DEMs) for landslide volume estimation. We investigated the Hítardalur landslide, which occurred on the 7 July 2018 in western Iceland, along with the geomorphological changes induced by this landslide, using optical and synthetic aperture radar data from Sentinel-2 and Sentinel-1. The results show that there were no considerable changes in the landslide area between 2018 and 2019. However, the landslide-dammed lake area shrunk between 2018 and 2019. Moreover, the Hítará river diverted its course as a result of the landslide. The DEMs, generated by ascending and descending flight directions and three orbits, and the subsequent volume estimation revealed that—without further post-processing—the results need to be interpreted with care since several factors influence the DEM generation from Sentinel-1 imagery.

**Keywords:** object-based image analysis; Sentinel-1; Sentinel-2; digital elevation model; InSAR; landslide; landslide-dammed lake; river; Iceland

#### **1. Introduction**

Landslide mapping and analysis are essential aspects of hazard and risk analysis, and the accurate detection of land surface changes is crucial for understanding processes and interactions between human and natural phenomena [1,2]. Landslides can be triggered by earthquakes, snowmelt, severe rainfall, human activities, or a combination of these factors. Landslides can partially or entirely block rivers and subsequently create landslide-dammed lakes, whereby potential dam failures pose a substantial risk for downstream areas [3]. Large, rapid mass movements are a common geomorphological process in Iceland and represent a significant threat to people and infrastructure [4–6]. Examples of such landslides are a rock avalanche onto the Morsárjökull outlet glacier in 2007 [7], a debris slide at the Móafellshyrna Mountain in northern Iceland in 2012 [8], and the landslide at the Askja caldera in July 2014, which caused a displacement wave [9–11].

Earth Observation (EO) data from optical and synthetic aperture radar (SAR) sensors have proven valuable for mapping and monitoring geomorphological features and, in particular, different types of landslides [1,12–14]. Optical imagery from different sources such as satellites, or unmanned aerial vehicles (UAVs) [15,16] has been used for landslide mapping and analysis [15,17]. However, mapping landslides using only optical imagery is challenging because of the spectral, spatial, and temporal characteristics of landslides [18]. SAR imagery has the potential to be used for hazard assessments by providing large-scale two-dimensional high spatial and temporal resolution images of the Earth's surface [19]. Radar pulses can penetrate through clouds (nearly weather independent), and they can provide data during the day and night (sun independent) [20]. According to Lee and Pottier [19], the surface reflectivity measured by radar imagery (also known as the radar backscatter coefficient σ0) is a function of the radar system parameters (such as the frequency, polarization, incident angle) and the surface parameters (such as the topography, roughness, dielectric properties of the medium, moisture). These parameters can be used to extract features from SAR imagery. Moreover, different polarizations provide different information about features of interest on the ground, and depending on the structure and position of the feature of interest, they might appear differently with varying polarizations [21]. The combined interpretation of optical and SAR imagery provides promising opportunities for landslide investigation while existing approaches still need to be improved [22–25].

Several studies have shown that object-based image analysis (OBIA) can achieve better classification accuracies than pixel-based techniques [26,27]. OBIA has been used for more than two decades as a framework for feature extraction from different imagery [28,29]. OBIA allows the use of spectral, textural and spatial information and thus provides a suitable framework for landslide mapping [30]. The integration of optical and SAR imagery within an OBIA approach for mapping geomorphological features is promising since the properties of both data types can be exploited in one interlinked workflow [24].

The estimation of landslide volumes is another key parameter for hazard analysis. There are several methods for estimating the volume of a landslide, of which the use of pre- and post-event digital elevation models (DEMs) can provide a good estimation [31]. However, accurate post-event DEMs are rarely freely available. Thus, the straightforward generation of post-event DEMs based on freely available data is important. There are several techniques to create DEMs, such as using optical imagery and stereoscopy or using SAR imagery and interferometric SAR (InSAR) techniques. The InSAR technique is based on the processing of at least two complex SAR images covering the same area and acquired from slightly different points of view [32]. InSAR allows the extraction of three-dimensional information using the phase difference and the along- and cross-track location of targets on the image and can be used for measuring the topography of a surface and its changes over time [33–35]. The generation of DEMs from SAR imagery using InSAR techniques is not new; for example, the well-known Shuttle Topography Mission (SRTM) provides near-global DEM data produced by InSAR [36]. However, there is a need for higher temporal and spatial resolution DEM data, particularly for assessing landslide events.

In this study, Sentinel-2 optical data and Sentinel-1 SAR data were used for geomorphological feature mapping. Sentinel-2 is a multispectral optical sensor from the European Union's Earth Observation programme, Copernicus, with 13 spectral bands combining different spatial resolutions (up to 10 m for Red, Green, Blue (RGB) and Near-infrared (NIR) bands) and a repeat cycle of 5 days [37]. The Sentinel-1 A SAR instrument was the first satellite of the Copernicus Sentinel missions [38], which provides radar data following the retirement of the ERS-2 and Envisat missions. The Sentinel-1 mission comprises a constellation of two polar-orbiting satellites (A and B). The Sentinel-1 instrument provides high-resolution C-band (wavelength of 5.6 cm) data with a short revisit time of 12 days for each satellite, and six days with Sentinel-1 A and B [39]. The satellite has a new type of ScanSAR mode, known as Terrain Observation with Progressive Scan (TOPS) SAR [40]. It has a similar pixel resolution (about 3.5 m in the range direction and 14 m in the azimuth direction) to ScanSAR but with an enhanced signal-to-noise ratio distribution [41]. The high temporal and spatial resolution of the Sentinel-1 A and B SAR has been used in various InSAR applications [42]. However, using Sentinel-1 data for generating DEMs and its application in landslide analysis is still in the early stages [43,44]. In this research, we focused on a large landslide that occurred on the 7 July 2018 on the eastern side of

the Fagraskógarfjall mountain in Hítardalur valley, western Iceland, leading to a landside-dammed lake and changes in the watercourse [45].

The main objectives of this research are to semi-automatically assess the landslide-induced geomorphological changes in Hítardalur valley by jointly using Sentinel-1 and Sentinel-2 data within an OBIA workflow and to assess the potential of post-event DEMs generated from Sentinel-1 image pairs for the estimation of the Hítardalur landslide volume.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Our study area is the Hítardalur valley and, in particular, the area surrounding the Hítardalur landslide, which is located on the eastern side of the Fagraskógarfjall mountain in the Vesturland region, western Iceland (Figure 1). The elevation in the Hítardalur valley ranges from 89 to 743 m above sea level. The valley floor is partly covered by a lava field (Hagahraun) [46]. Thick sedimentary layers can be found within the stratigraphic sequence. The rock is rather dense as the cavities have been filled with minerals. Quaternary volcanic formations, created during the last glaciation and Holocene, can be found within the vicinity of the landslide. Faults and fractures are prominent within the region and can be divided into three groups, i.e., NE–SW, N–S and NW–SE oriented [47]. The faults and fractures that belong to the NE–SW and N–S groups formed mainly during the formation of the lava pile and are older than 8 million years. NW–SE oriented fractures belong to the so-called Snæfells fracture zone which extends from Kerlingarskarð in the north to Borgarfjördur valley in the south. The NW–SE fracture zone was very active 8 to 4.5 million years ago and is considered to be still active [48].

**Figure 1.** (**a**) The study area around the Hítardalur landslide shown on a Sentinel-2 image acquired on 17 July 2018 and (**b**) its location in western Iceland; background data © ESRI.

On the morning of the 7 July 2018, a large landslide blocked the Hítará river and led to the creation of a landslide-dammed lake. Water from the landslide-dammed lake found a new way through the highly porous lava field. The landslide originated in an area of the Fagraskógarfjall mountain that shows evidence of earlier displacements [45]. The year 2018 had more rainfall than usual in Iceland [49], likely leading to a higher water pressure in cracks and fractures that further weakened the rocks, and thus contributed to landslide initiation [50]. Helgason et al. [45] estimated that approximately 7 million m<sup>3</sup> of material was released and that a total of 10–20 million m<sup>3</sup> of material was displaced. The landslide deposition covers an area of approximately 1.5 km2 on the valley floor and is up to 30 m thick. The fall height was about 450 m over a runout length of approximately 2.3 km at a runout angle of approximately 12◦. The length of the riverbed covered with debris was estimated to be around 1.6 km. The Hítardalur landslide is considered to be one of the largest landslides in Iceland in recorded history [51]. Figure 2 shows field photographs of the Hítardalur landslide, the landslide-dammed lake and the river finding a new path after emerging from the dammed lake.

**Figure 2.** (**a**) Photograph showing the Hítardalur landslide and the landslide-dammed lake soon after the event occurred, i.e., taken on the afternoon of 7 July 2018 (photograph: © S. Asgeirsson). (**b**) View from the landslide deposition area towards the landslide source area (photograph: © J. K. Helgason; 7 July 2018). (**c**) Water from the dammed lake finds a new way (photograph: © J. K. Helgason; 7 July 2018). (**d**) View from the Fagraskógarfjall mountain towards the Hítardalur valley showing the landslide deposition area and the new watercourse flowing out from the lake into the Stekka river (photograph: © M. Olafsson; 23 November 2018).

#### *2.2. Data and Data Preparation*

#### 2.2.1. Optical Data

We used pre- and post-landslide event Sentinel-2 data (Table 1). The data selection was made based on the cloud cover estimates (<30%) provided by the European Space Agency (ESA) for pre- and post-event datasets, using the Google Earth Engine (GEE) platform.


**Table 1.** Sentinel-2 datasets used in this research. MSI stands for Multispectral Instrument.

#### 2.2.2. Radar Data

We used the Sentinel-1 Interferometric Wide Swath (IWS) Level-1 Single Look Complex (SLC) product. The SLC is the main product of the IWS mode. It is mainly used for interferometric applications [38,52]. The Sentinel-1 IWS SLC product includes three sub-swaths and each of them includes one intensity image per polarization channel, and several bursts over each of the sub-swaths [40,52]. However, for mapping purposes, the SLC data should be post-processed and converted to the Ground Range Detected Geo-referenced product (GRD). For mapping geomorphological features and related changes, backscatter coefficients (also known as sigma naught, or sigma 0) were created by converting the Sentinel-1 SLC images with slant-range geometry to Sentinel-1 GRD products. The GRD data consist of focused SAR data that have been corrected using SLC data by applying a calibration and speckle noise filtering using the Lee refined filter [53,54], and a ground geometry correction using an Earth ellipsoid model. The thermal noise of the GRD products was removed using a noise look-up table provided for each image to derive the calibrated noise profiles matching the calibrated GRD data and to improve the quality of the TOPSAR images [55,56]. A visual inspection of several Sentinel-1 scenes from different orbit tracks and the two available flight directions (ascending with the satellite moving north, and descending with the satellite moving south) showed that ascending images were more suitable than the descending images for mapping the geomorphological features and related changes in Hítardalur valley due to their geometry. The whole landslide, including both the landslide source and the deposition areas, is usually well visible on the ascending image, whereas only the deposition area is visible on the descending image. Therefore, one pre-event (5 July 2018) and two post-event (17 July 2018 and 5 August 2019) ascending Sentinel-1 TOPSAR scenes (SLC type) in the interferometric wide (IW) mode, with vertical-vertical (VV) and vertical-horizontal (VH) dual polarizations from track 16 were selected for the geomorphological mapping.

Figure 3 shows the combination of Sentinel-1 sigma naught VV polarization and the Sentinel-2 datasets for each year, which were used for geomorphological feature mapping.

Moreover, we used six additional Sentinel-1 IWS SLC datasets for generating post-event DEMs using InSAR. The quality of DEMs generated from InSAR analysis essentially depends on the perpendicular component of the spatial baseline (known as the perpendicular baseline or Bperp) between the radar antenna locations, divided by the distance to the ground. The distance to the ground is similar for the orbiting satellites, so the spatial baseline is the important variable. If the perpendicular baseline is too small, the topographic effects on the differential phase are not pronounced enough. On the other hand, if the baselines are too big, the coherent phase is increasingly different, leading to decorrelation. For example, for the ERS satellite, a suitable Bperp for DEM generation is recommended to be between 150 and 300 m at the time of the image acquisition [57]. This grants an angle between both acquisitions, which allows the retrieval of topographic variations from parallel-like effects. This Bperp was originally recommended for the ERS data, while for Sentinel-1 a lower range is usually used according to the literature. For example, Geudtner et al. [58] described that due to the orbit maintenance strategy, the baselines for the Sentinel-1 A and B are mostly in the range of ±150 m. Kyriou et al. [59] used a Bperp between 96 and 170 m for landslide mapping with Sentinel-1. For this study, we downloaded the available Sentinel-1 images with a Bperp range according to these recommendations for 2018 and 2019 (only from June to September to avoid potential snow cover) and processed them to generate post-event DEMs (Section 2.4).

**Figure 3.** Sentinel-1 and Sentinel-2 datasets used for geomorphological features mapping. The upper row shows the intensity layer for Sentinel-1 vertical-vertical (VV) polarization for the following image acquisition dates: (**a**) 5 July 2018; (**b**) 17 July 2018; (**c**) 5 August 2019. The lower row shows the true color composites of the Sentinel-2 images with the following acquisition dates: (**d**) 20 June 2018; (**e**) 17 July 2018; (**f**) 1 August 2019.

Likewise, the height ambiguity (*ha*) was used as an indicator of the accuracy of the topographic height, which generates the 2-π InSAR phase change (Equation (1)), whereby the higher the Bperp, the more accurate the altitude measurement [57,60]:

$$
\hbar \hbar a = \frac{\lambda}{2B\_{\rm perp}} \tag{1}
$$

where *ha* is the height ambiguity, λ is the radar wavelength, *R* is the range from satellite to ground, and θ is the look angle.

From the above equation, it is clear that *ha* is proportionally related to the Bperp, and a higher Bperp theoretically should result in a more accurate altitude measurement. However, for a longer baseline, it is more difficult to solve the 2π phase ambiguity by phase unwrapping than for a shorter baseline [61].

The temporal interval between two SAR images was another factor to consider for the image selection. It influences the coherence, which determines how accurately the phase can be measured. A low coherence indicates a higher phase noise, which can cause larger elevation errors and phase unwrapping errors. A maximum relative coherency is recommended to be considered for the selection of SAR data pairs for DEM generation. Moreover, the interferometric phase depends on the orbit geometry, scene topography, line-of-sight surface displacements, and atmospheric path delays [62].

We used the Alaska Satellite Facility Baseline Tool to identify and select image pairs. From the available orbit tracks (i.e., number 16, number 118, and number 155), and flight directions (descending and ascending), we selected three post-event Sentinel-1 (A and B) image pairs (Table 2) for DEM generation and volume estimation. We used the InSAR stack overview function implemented in the SNAP (Sentinel Application Platform) toolbox to estimate the expected coherency (called "modeled coherency"), and the height ambiguity between two TOPSAR pair images.

**Table 2.** Sentinel-1 data and their characteristics used to generate the digital elevation models (DEMs). The Bperp stands for the spatial perpendicular baseline, which was used as the main criterion for the selection of Sentinel-1 A and B image pairs. The modeled coherency and height ambiguity measurements were derived using the interferometric synthetic aperture radar (SAR) (InSAR) stack overview function implemented in the Sentinel Application Platform (SNAP) toolbox.


All SAR processing was completed using the open-source SNAP toolbox provided by ESA. Moreover, we used existing global earth topography and sea surface elevation data at 30 arc-second resolution (GETASSE30), which was freely available and accessible in the SNAP toolbox in the DEM generation process.

#### 2.2.3. Topographic Data

We used the freely available ArcticDEM (spatial resolution of 2 m), provided by the Polar Geospatial Center, for the validation of the generated DEMs. The ArcticDEM dataset is generated based on overlapping very high-resolution (VHR) optical satellite images with sub-meter resolution and stereoscopic imagery [63,64].

#### *2.3. Geomorphological Features Mapping*

In the mapping of geomorphological features, we focused on the Hítardalur landslide, the landslide-dammed lake, and the riverbed with flowing water. Additionally, we differentiated between the landslide source and deposition area and identified changes in the watercourse based on the pre-and post-event images. We used an OBIA approach and created a knowledge-based classification ruleset in eCognition (Trimble) software, integrating Sentinel-1 and Sentinel-2 datasets. We used spectral indices derived from the Sentinel-2 data, including a brightness layer (average of the three visible spectral bands blue (B2), green (B3), red (B4)) and the normalized difference vegetation index (NDVI). For the detection of riverbeds, we calculated an edge extraction layer based on the Sentinel-2 NIR (B8) spectral band using the Lee sigma [65] algorithm. We also used the intensity information from the pre- and post-event Sentinel-1 VV, and VH polarization intensity images to create subtraction layers for 2018 and 2019. The Sentinel-1 intensity subtraction layers and the edge extraction layers derived from band B8 were used for the object-based classification for 2018 and 2019 and are shown in Figure 4.

**Figure 4.** Examples of the layers extracted from the Sentinel-1 and Sentinel-2 datasets. (**a**) and (**b**) show the subtraction layers using pre- and post-event Sentinel-1 VV polarization intensity images for 2018 and 2019, respectively. The Hítardalur landslide and the landslide-dammed lake are distinguishable. (**c**) and (**d**) show the edge extraction layers derived from the post-event Sentinel-2 images (Near-infrared (NIR) band (B8)) for 2018 and 2019, respectively.

The Hítardalur landslide and landslide-dammed lake are distinguishable in the Sentinel-1 subtraction layer of the year 2018 (Figure 4a), but they are less obvious on the Sentinel-1 subtraction layer from 2019 (Figure 4b). Due to the SAR amplitude changes, the Hítardalur landslide appears very bright, and the landslide-dammed lake very dark on the Sentinel-1 subtraction layer of 2018 (Figure 4a). The white path visible on the Sentinel-1 VV polarization subtraction layer of 2018 (Figure 4a) is the Hítará riverbed (i.e., towards the south-west of the landslide). In the Sentinel-1 subtraction layer of the year 2019, the Hítará river change is hardly visible (Figure 4b). The Sentinel-2 edge extraction layer of the year 2018 clearly shows the new watercourse between the dammed lake and Stekka river (i.e., the river towards the south-east of the landslide) (Figure 4c). This connection, however, is not readily distinguishable on the edge extraction layer of the year 2019. In both images (Figure 4c,d), the outlines of the lake are partly visible, but, it is difficult to differentiate them from the edges of other features. The white circular area in Figure 4d is a cloud mask.

The OBIA classification of the landslide was based on the pre- and post-event VV polarization subtraction layer, the brightness layer was used for the classification of the landslide-dammed lake, and the riverbed with water was classified using the edge extraction layer and refined based on spectral and spatial parameters. The first step in the object-based mapping was the creation of image objects using the multiresolution segmentation algorithm. Then, we created a knowledge-based classification ruleset which includes thresholds based on the backscatter information from Sentinel-1 and the spectral information from the Sentinel-2 images for each year. The size of the segmentation-derived image segments is controlled by the scale parameter (SP), which is directly related to the local variations of the bands used in the multiresolution segmentation algorithm [66,67]. In general, the higher the value of the SP, the larger the resulting segments, and vice versa. First, multiresolution segmentation was applied at the pixel level to create suitable objects for the subsequent classification of the landslide, using the Sentinel-1 VV intensity change layer, and the landslide-dammed lake, using the Sentinel-1 VV intensity change and NDVI layers. The same features were used for classification in both years, whereby some adaptations in the thresholds for the 2019 mapping were required. The difference in the values of the Sentinel-1 VV polarization change layers between the year 2018 and 2019 could be related to the different time of the Sentinel-1 image acquisitions, and potential changes in the surface structure of the landslide, as well as different dielectric properties (e.g., moisture content) of the landslide. Another multiresolution segmentation was performed based on the unclassified segments from the first segmentation and using only the Sentinel-2 post-event images to create suitable objects for the

classification of the riverbeds with water. We excluded the Sentinel-1 datasets from this classification since riverbeds were barely visible in the datasets.

Moreover, we classified the original riverbeds with water using only the pre-event Sentinel-2 image (20 June 2018). Therefore, we used the edge extraction layer sigma >5 followed by a manual refinement, merging, and using a threshold for the object length of >470 pixels.

Figure 5 gives an overview of the OBIA mapping workflow and the features and thresholds which were used for the post-event geomorphological feature mapping for 2018 and 2019.

**Figure 5.** The classification framework used for the post-event geomorphological feature mapping for 2018 and 2019.

Moreover, we semi-automatically delineated the part of the landslide where the material was deposited in the year 2018 as an input for the subsequent landslide volume estimation. Therefore, we re-segmented the landslide class. Then the slope information derived from the Arctic DEM was used to differentiate the landslide deposition area from the landslide source area with a slope threshold of <12◦. The value of 12◦ was selected based on the work of Helgason et al. [45], impressions from the field, and considering the visual comparison of the field photographs with the remote sensing data.

#### *2.4. DEM Generation*

DEMs are an important component of landslide modelling, volume calculation, and landform classification. In this research, we performed an exploratory assessment of the usability of Sentinel-1 SLC TOPSAR datasets for the generation of post-event DEMs within an automated workflow. The three SLC TOPSAR image pairs were used to create post-event DEMs using InSAR. Each TOPSAR image was calibrated using precise orbit ephemerides (POD) products provided by the ESA. The POD products contain auxiliary information about the position of the satellite during the data acquisition. For computational efficiency, the split TOPSAR function was applied to each TOPSAR image to select the swath sub-bursts covering the area of interest. The TOPSAR images require a precise coregistration to ensure each ground target contributes to the same (range and azimuth) pixel in both images [57]. The coregistration of SLC TOPSAR images can be completed on a pixel by pixel basis, with an accuracy in the order of one-tenth of the resolution. Therefore, the process of DEM-assisted coregistration was applied using the back-geocoding function. Due to the steep azimuth spectrum ramp in each burst [42] the correction of the shift in the azimuthal direction was performed by enhanced spectral diversity [68]. Then, an interferogram generation was carried out by cross multiplying the complex conjugate of two coregistered TOPSAR images. During this process, the amplitude of both images is multiplied while the phase represents the phase difference between the two TOPSAR images. The flat-Earth phase component was then subtracted to remove the phase present in the interferometric signal due to the curvature of the reference surface. The flat-Earth phase component was estimated using the orbital- and metadata information [69]. The resulting interferometric fringes represent a full 2 π cycle, with arbitrary colors representing half the sensor's wavelength. In the case of repeat-pass acquisitions, such as in the case of Sentinel-1, the interferograms are rather noisy, mainly due to temporal decorrelation. To improve the quality of the interferograms and reduce the phase noise of the adjacent pixels, the Goldstein [70] filtering with fast Fourier transformation (FFT) was used. We tried different settings and window sizes of the FFT filter and finally selected 64 for the FFT filter definition and 3 × 3 for the window size. The multi-look speckle filtering was applied to increase the radiometric accuracy by reducing the local variability of the signal [71]. The multi-look speckle filtering improves the signal-to-noise ratio at the expense of the spatial resolution [72]. Moreover, applying a multi-look speckle filtering on the interferogram increases the computational efficiency for phase unwrapping. The phase unwrapping is necessary to resolve the unknown multiple-of-wavelength ambiguity in the interferometric phase. We used the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU), which is a two-dimensional phase unwrapping algorithm [73–75]. The SNAPHU algorithm is written in *C*, can be used within the SNAP toolbox, or as a standalone program on a Linux platform. The unwrapped phase values were converted into elevation values with respect to a reference ellipsoid [76]. This process specifies each pixel in the unwrapped phase image with respect to a Cartesian reference system using radar coordinates (range, azimuth, phase variation). The last step in the DEM generation was data geocoding to allow for the comparison of the results with a reference DEM.

The OBIA-delineated landslide deposition area for 2018 was used for the calculation of the landslide volume. The volume estimation was conducted by subtracting each of the three post-event DEMs from the pre-event ArcticDEM.

#### *2.5. Validation*

The accuracy of the geomorphological feature mapping was assessed by comparing it to mapping results available from the National Land Survey of Iceland (NLSI) [77], and information provided by Helgason et al. [45].

As for the accuracy assessment of the DEMs generated using Sentinel-1 image pairs, we used the high-resolution ArcticDEM as a reference. This is in line with the literature, which recommends that the resolution of the reference data should be at least three times higher than the resolution of the DEM elevation being evaluated [78,79]. The ArcticDEM has known biases of several meters due to errors in the sensor models due to the satellite position and different sensors in the ArcticDEM constellation. However, its vertical accuracy is corrected by registering it to the Ice, Clouds, and Land Elevation Satellite mission (ICESat) altimetry information and the results can be used without further corrections [80]. The DEM quality assessment was conducted by comparing the statistical measures, such as the minimum, maximum, mean, and standard deviation, of the generated DEMs with those of the reference ArcticDEM. The vertical quality assessment was completed using the root mean square error (RMSE), whereas the horizontal quality assessment was completed using autocorrelation and the Moran's I index [81,82].

The RMSE includes both random and systematic errors introduced during data production [83], and is expressed by the following equation:

$$RMSE = \sqrt{\frac{\sum \left(y\_i - y\_{ti}\right)^2}{N}} \tag{2}$$

where *yi* refers to the ith interpolated elevation, *yti* refers to the ith known or measured elevation of a sample point in a reference dataset, and *N* is the number of sample points.

While the RMSE is a valuable quality-control statistic, it does not provide an accurate assessment of how well each cell in the DEM represents the true elevation. The RMSE only provides an assessment of how well the DEM corresponds to the data to which it is compared [84]. The Moran's I index expresses local homogeneity by comparing the difference between neighboring pixels to the standard deviation. The Moran's I index ranges between +1 and −1, where +1 indicates a strong spatial autocorrelation, 0 a spatially uncorrelated data, and −1 a strong negative spatial autocorrelation. The Moran's I was calculated by adopting the "Queen's" rule, which takes into account the eight neighboring pixels. The generated DEMs from the InSAR analysis were first subtracted from the ArcticDEM, and then the Moran's I index was calculated on the difference map. We consider only the "stable area" for applying the Moran's I and RMSE measures, i.e., the area without the landslide and the landslide-dam lake.

#### **3. Results**

#### *3.1. Geomorphological Features Mapping Using OBIA*

The results of the geomorphological features mapping are shown in Figure 6. The pre-event result shows the original Hítará riverbed with water (Figure 6a). The post-event results (Figure 6b for 2018 and Figure 6c for 2019) show the landslide, the landslide-induced changes in the watercourse, and the landslide-dammed lake. The reference map from NLSI [77], which was used for comparison, is shown in Figure 6d. The water flowing out of the dammed lake found a new route through a lava field towards the south and then merged with the Stekka river, resulting in more water flowing in this riverbed than prior to the landslide. However, the upper part of the connection between the lake and the Stekka river is not easily distinguishable in the Sentinel-2 image taken in August 2019, probably due to less water flowing there at the time of the image acquisition one month later in the year 2019 (July 2018 versus August 2019). It seems that the Hítará river partly passes through the landslide and appears as a spring in the river bed downstream, i.e., approximately 2.5 km downstream from the lake and approximately 1 km downstream from the closest part of the landslide deposition. This is confirmed by the manually created reference data from NLSI [77]. The landslide area was estimated to be approximately 2000 ha in 2018 and 2019. Slight differences between the landslide areas likely result from variations in the segmentation. We estimated a lake area of 58 ha in 2018 and 47.1 ha in 2019. The reference values for the lake area are reported to be approximately 47 ha by Helgason et al. [45], but, the lake area from the shapefile provided by NLSI is 39.8 ha [77]. A comparison of our results with these numbers reveals an overestimation of the OBIA result for the year 2018. This could be associated with a high moisture content near the lakeshore and shallow water areas with high sediment load, which introduces uncertainty in the segmentation and classification [17]. The lower value for 2019 is probably due to the warm summer in 2019, which caused the discharge of some rivers to be lower for a few weeks [85].

The maximum width of the landslide deposition is approximately 1.5 km, and the overall runout length of the landslide is approximately 2.3 km. These results are in line with the numbers reported by Helgason et al. [45].

**Figure 6.** Object-based geomorphological feature mapping using Sentinel-1 and Sentinel-2 datasets. (**a**) shows the riverbed with water before the Hítardalur landslide using the pre-event Sentinel-2 image only; (**b**) shows the landslide, the landslide deposition area, the landslide-dammed lake, and the (now partly interrupted) riverbeds with water after the landslide event in 2018; (**c**) shows the landslide, the landslide-dammed lake, and the riverbeds with water in 2019; (**d**) shows the reference map from National Land Survey of Iceland IS 50V 24 September 2019 Vatnafar Flakar (© NLSI) [77]. The hillshade derived from the ArcticDEM was used as a background layer.

#### *3.2. DEMs from Sentinel-1 and Landslide Volume Estimation*

Figure 7 shows the results of the post-event DEMs generated using the Sentinel-1 image pairs from different orbit tracks, and the pre-event ArcticDEM for comparison.

**Figure 7.** Illustration of the three post-event DEMs derived from the Sentinel-1 image pairs. (**a**) for orbit track number 16 (11 July 2018 and 17 July 2018), ascending; (**b**) orbit track number 118 (5 August 2018 and 11 August 2018), ascending; (**c**) orbit track number 155 (4 July 2019 and 10 July 219), descending; (**d**) ArcticDEM. For visualization purposes, each DEM is overlaid with the hillshade derived from it.

For the volume estimation of the landslide, each of the three generated DEMs was subtracted from the ArcticDEM. Then a spatial subset for the landslide was created using the OBIA delineation of the landslide deposition area for 2018 (Figure 8). The DEM from track 16 shows maximum elevation differences to the ArcticDEM of more than 80 m, whereas the other two DEMs (track number 118 and 155) show lower, more homogenous and realistic elevation differences in the landslide deposition area. However, both DEMs from tracks 118 and 155 partly show negative values upslope, which indicate the removal of material and potential issues in the delineation of the deposition area.

**Figure 8.** Elevation differences of the three Sentinel-1 DEMs compared to the ArcticDEM for the landslide deposition area used for calculation of the landslide volume. (**a**) elevation differences of DEM from track 16 (11 July 2018 and 17 July 2018); (**b**) elevation differences of DEM from track 118 (5 August 2018 and 11 August 2018); (**c**) elevation differences of DEM from track 155 (4 July 2019 and 10 July 2019). The respective hillshade derived from the Sentinel-1 DEMs was overlaid for each result.

Table 3 shows the landslide deposition volumes. The DEM from track 16 gives a much higher volume than the other two DEMs. The landslide volume estimations using the DEMs from track 118 and 155 are in a similar range.



#### *3.3. DEM Validation*

Table 4 shows the measures used for the accuracy assessment of the generated DEMs. The comparison of the statistical measures shows an underestimation of the minimum DEM value for the track 16 DEM, while the other two DEMs show similar minimum statistics as the ArcticDEM. In contrast, the maximum value for the track 16 DEM is the closest compared to the ArcticDEM while the other two DEMs (tracks 118 and 155) show lower values. The mean and standard deviation values for the track 16 DEM are closer to the values of the ArcticDEM.

The RMSE for the track 16 DEM is the lowest of the three DEMs, which indicates less vertical variation between track 16 and the ArcticDEM, but at the same time it has the lowest Moran's I autocorrelation (0.92), which could be a result of the existence of more outliers compared to the track 118 and 155 DEMs, both of which have a Moran's I of 0.97. The latter DEMs have better autocorrelation due to more homogenous values.

It is worth mentioning that both quality assessment indicators used (i.e., RMSE and Moran's I) are global indicators, and they don't take into account local variations.


**Table 4.** Validation of the Sentinel-1 DEMs with respect to the ArcticDEM. Comparison is based on the area not influenced by the landslide. The Moran's I statistic is calculated on the difference map.

From the statistics in Table 4, it is challenging to choose the most suitable DEM for volume estimation. Despite a clear indication of which DEM should give the best results, we see major differences in the measured volumes for the three DEMs (track 16, 118, and 155). According to the available literature [45], approximately 7 million m3 of material was released from the landslide source area, and the volume of the debris in the depositions area is about 10–20 million m3. According to these statistics, the track 16 DEM shows an unrealistic value of 109 million m3, whereas the landslide volumes calculated from the other two DEMs are closer to the estimated reference value.

#### **4. Discussion**

During the last two decades, OBIA has been used for many applications [28,29]. This research shows the advantage of integrating backscatter information from Sentinel-1 SAR data and spectral information from Sentinel-2 optical data within an expert knowledge-based classification approach in an OBIA environment to assess landslide-induced geomorphological changes. We used Hítardalur as a case study, where a large landslide blocked the river, and, consequently, a lake formed and the river changed its direction. The OBIA classification results show no significant change in the landslide area in the year after the event, but changes in the area of the landslide-dammed lake were detected. However, this might be partly explained by uncertainties in image segmentation and classification and seasonal variations due to different acquisition dates (July 2018 versus August 2019). Moreover, the spatial resolution of Sentinel-1 and Sentinel-2 images are limiting factors for river classification. Therefore, we first identified the existing riverbeds based on an edge detection layer and then assessed whether water was present in the riverbeds.

The Hítardalur landslide has several interesting characteristics that make it suitable to be considered as a case study for creating DEMs using Sentinel-1 datasets. It is a large landslide with a deposition area spreading over the rather flat valley. Therefore, both ascending and descending Sentinel-1 datasets can be used. The resolution and accuracy of DEMs have a direct influence on subsequent analyses such as the landform classification and hydrological modeling [86]. We applied a straightforward workflow for generating DEM using Sentinel-1 image pairs from three different tracks. Our results are heterogeneous and indicate that the generation of DEMs has several limitations in this situation. InSAR techniques have been extensively used for DEM generation using SAR datasets, but the accuracy of the resulting DEMs depends on several factors related to the data and processing workflows and techniques. For selecting suitable SAR datasets, the orbit indetermination and Bperp range, atmospheric conditions, and temporal decorrelation (due to scene changes between the two passes) need to be considered. To achieve a good vertical accuracy, the SAR image pair requires a large Bperp [76]. The Sentinel-1 mission was mainly designed for the retrieval of surface deformations using differential InSAR (DInSAR) which requires a low Bperp between consecutive SAR images, and DEM generation was not the primary goal of the mission [58]. The recommended Bperp for the

DEM generation using the ERS satellite is between 150 and 300 m [57,87], but the ERS satellite is known to have high Bperp, whereas the common Bperp for Sentinel-1 datasets is below 30 m, which is very suitable for interferometry analysis but challenging for DEM generation [88,89]. In this study, we selected datasets in the Bperp range of 130 to 160 m, with a short temporal baseline of 6 days between each image pair. The temporal baseline is important to reduce phase noise. Although radar imagery is known for its all-weather measurement capability, changes in the atmospheric conditions (mainly water vapor) cause a variable path delay which results in atmospheric distortions and is the main error source in repeat-pass SAR interferometry [62]. In Iceland, the atmospheric and weather conditions make it difficult to find Sentinel-1 image pairs with a suitable Bperp and short temporal interval [90]. Therefore, we limited our search to the summer months (June to September), also to avoid snow cover. We also considered the coherency and the height ambiguity in our data selection. We selected three SAR image pairs from different orbit tracks, i.e., track 16 (11 July 2018 and 17 July 2018) and track 118 (5 August 2018 and 11 August 2018) for the ascending flight direction and track 155 (4 July 2019 and 10 July 2019) for the descending flight direction. One of the selected image pairs was in the range of the recommended Bperp, i.e., track number 155 with a Bperp of 159. This image pair also showed the lowest height ambiguity (98) compared to the others. The other two image pairs had a slightly lower Bperp (134 for track 16 and 142 for track 118) than recommended, but when considering the other selection criteria, they seemed to be the best datasets available.

We faced challenges in validating our classification and DEM-generated results. One issue was the availability of reference data. We used the values provided by Helgason et al. [45] as a reference for our classification and landslide volume estimation. Our volume estimations need to be considered with care since there are several uncertainties associated with them, and certain limitations make the comparison difficult. For example, the area used for the calculation of the landslide volume is not entirely clear and most likely differs from the delineation we used. We used the landslide deposition area, which was semi-automatically extracted using OBIA. Therefore, a one to one comparison of our results with other results reported in literature is not possible. Another challenge is directly related to the reliability of using the generated DEMs for volume calculations. We directly used the generated DEM from the InSAR processing workflow for volume estimation without any further post-processing, while the literature suggests that the InSAR-produced DEMs should be further post-processed to be ready for the end-user [42]. However, since we aimed to test the suitability of DEMs for landslide volume estimation created through a straightforward workflow that could be directly transferred to other areas, we did not apply any post-pressing to improve the volume estimation results. In this study, we demonstrated the potential and limitations inherent to the approach and the Sentinel-1 datasets for DEM generation. Thus, we did not apply any vertical correction on the DEMs generated from Sentinel-1 datasets. We did not use any pre-event DEM generated with this workflow. This was mainly to avoid introducing further errors to the volume estimation. However, there are several possibilities to deal with error distribution in the InSAR-generated DEMs, such as the fusion of the final products from the ascending and descending flight directions [32,91] and the horizontal and vertical adjustment of pixels in the post- and pre-event DEMs, to match them better and make the direct comparison more reliable [92]. Further work is needed to fully evaluate such approaches and the potential of using Sentinel-1 data to generate post-landslide event DEMs for volume estimation.

#### **5. Conclusions**

The combined use of Sentinel-1 and Sentinel-2 data offers opportunities for assessing landslides and landslide-induced geomorphological changes. Although many applications use EO datasets for landslide mapping [18,93], studies on integrating intensity information derived from Sentinel-1 SAR data with spectral information from Sentinel-2 optical data are still rare. In this study, we presented a semi-automated workflow to map the landslide-induced geomorphological changes in Hítardalur valley, Iceland, by jointly using Sentinel-1 and Sentinel-2 data within an OBIA framework and assessed the potential of post-event DEMs generated from Sentinel-1 image pairs from different orbit tracks for the estimation of the Hítardalur landslide volume.

It is plausible that we will face more frequent and large landslides in the future due to climate change. The high spatial and temporal resolution, free availability, and the integrated and effective use of Sentinel-1 and Sentinel-2 datasets can contribute enormously to analyzing such events and their consequences. However, the quality of post-event DEMs derived from Sentinel-1 data depends on several factors. These include the SAR image pair selection considering the perpendicular baseline, high coherency, and atmospheric conditions, as well as the process for DEM generation, for example, precise coregistration, phase unwrapping, and conversion of the unwrapped phase into elevation values. These are still active research areas [42]. Further research needs to be conducted to systematically assess the accuracy of the DEMs generated using Sentinel-1.

**Author Contributions:** Conceptualization: Z.D., D.H. and L.A.; Methodology: Z.D., D.H., L.A. and D.T.; Validation: Z.D., D.H., L.A. and J.K.H.; Formal Analysis: Z.D.; Investigation: Z.D., D.H., L.A., J.K.H. and Þ.S.; Data Creation, Z.D., D.H. and L.A.; Writing—Original Draft Preparation: Z.D. and D.H.; Writing—Review and Editing: Z.D., D.H., L.A., D.T., J.K.H. and Þ.S.; Visualization: Z.D., D.H., L.A.; Supervision: D.H.; Project Administration: D.H.; Funding Acquisition: D.H. and D.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by the Austrian Science Fund (FWF) through the project MORPH (Mapping, monitoring and modelling the spatio-temporal dynamics of land surface morphology; FWF-P29461-N29) and the Doctoral Collage GIScience (DK W 1237-N23), as well as by the Austrian Academy of Sciences (ÖAW) through the project RiCoLa (Detection and analysis of landslide-induced river course changes and lake formation).

**Acknowledgments:** The authors would like to thank S. Asgeirsson and M. Olafsson for providing field photographs.

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

**Data Availability:** Figure 6d shows the reference map from the National Land Survey of Iceland (NLSI) vector layer, available online at https://www.lmi.is (accessed on 15 May 2020); © NLSI (license: https://www.lmi.is/en/ licence-for-national-land-survey-of-iceland-free-data (accessed on 15 May 2020)).

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## **Bayesian Updating of Soil–Water Character Curve Parameters Based on the Monitor Data of a Large-Scale Landslide Model Experiment**

**Chengxin Feng 1, Bin Tian 1,\*, Xiaochun Lu 1, Michael Beer 2,3,4,\*, Matteo Broggi 2, Sifeng Bi 2, Bobo Xiong <sup>1</sup> and Teng He <sup>1</sup>**


Received: 30 June 2020; Accepted: 6 August 2020; Published: 10 August 2020

**Abstract:** It is important to determine the soil–water characteristic curve (SWCC) for analyzing landslide seepage under varying hydrodynamic conditions. However, the SWCC exhibits high uncertainty due to the variability inherent in soil. To this end, a Bayesian updating framework based on the experimental data was developed to investigate the uncertainty of the SWCC parameters in this study. The objectives of this research were to quantify the uncertainty embedded within the SWCC and determine the critical factors affecting an unsaturated soil landslide under hydrodynamic conditions. For this purpose, a large-scale landslide experiment was conducted, and the monitored water content data were collected. Steady-state seepage analysis was carried out using the finite element method (FEM) to simulate the slope behavior during water level change. In the proposed framework, the parameters of the SWCC model were treated as random variables and parameter uncertainties were evaluated using the Bayesian approach based on the Markov chain Monte Carlo (MCMC) method. Observed data from large-scale landslide experiments were used to calculate the posterior information of SWCC parameters. Then, 95% confidence intervals for the model parameters of the SWCC were derived. The results show that the Bayesian updating method is feasible for the monitoring of data of large-scale landslide model experiments. The establishment of an artificial neural network (ANN) surrogate model in the Bayesian updating process can greatly improve the efficiency of Bayesian model updating.

**Keywords:** large-scale landslide model experiment; soil–water characteristic curve; Bayesian updating; Markov chain Monte Carlo; artificial neural networks

#### **1. Introduction**

The prediction of a landslide and its stability performance have always been key difficulties of engineering. Particularly for the sensitive and landslide-prone areas in reservoir areas or mountain areas, the water migration inside a landslide is very prominent for its stability [1–3]. The soil–water characteristic curve (SWCC) is the cornerstone of the seepage calculation of unsaturated soil. In the analysis of saturated-unsaturated seepage, it is very important to fully describe the function of the slope of the moisture content [4,5]. The preliminary prediction based on the survey data often deviates from the actual monitoring value. As the real conditions are very complicated, the mechanism cannot be fully grasped through analytical models. Therefore, the trial algorithm is often used to repeatedly correct the soil parameters to obtain results that match the monitoring value [6,7]. The Bayesian method provides an effective way to solve such problems. By using the detection data, a priori information, and model assumptions, the parameter update is realized and the parameter uncertainty is quantified [8,9]. This technique uses a set of output measurement data (such as a point in a multi-dimensional output space) [10,11]. The purpose of this paper is to solve the problem of how to check the parameters of a SWCC model of a landslide using the Bayesian updating method.

The SWCC is a crucial input for modeling the geotechnical problems with unsaturated soil. It can be describing the relationship between matric suction and volumetric water content in unsaturated soils [12]. The main feature of SWCC is that the slope of the residual water content or saturation function in the pore water pressure range is from negative to positive [13]. As a basic equation in unsaturated soil mechanics, it plays an important role in the application of unsaturated soil mechanics for studying the strengths, deformations, and permeabilities of soils [14]. Past studies have shown that the shear strength of unsaturated soils can be estimated by the SWCC [15–17]. It can be seen that these parameters play important roles in the safety of geotechnical engineering [18–21]. Therefore, the SWCC parameters are evaluated as a method to determine the unsaturated soil properties of landslides. In engineering practice, the SWCC is usually measured by field or laboratory tests. The typical methods used are: the tensiometer method, the pressure plate meter method, the cold mirror hygrometer method, and the filter paper method. The limited number of data points are then used to estimate the SWCC by best fitting them with some parametric SWCC models, such as the van Genuchten–Mualem model (VGM), the van Genuchten-Burdine model (VGB), the van Genuchten model (VG), and the Fredlund and Xing model (FX) [14,22].

This paper adopts the Bayesian updating method, which is essentially a Bayesian statistical reasoning problem [23], to study the uncertainty of the parameters involved in the model. The Bayesian method is a probabilistic method based on the application of Bayesian rules. The resulting solution is called the posterior distribution, which is a multi-dimensional probability distribution function in the update parameter space. The posterior distribution is usually a complex distribution, and in principle it is integral in the parameter space. It can be approximated by statistical methods and using modern algorithms, namely, the Markov chain Monte Carlo (MCMC) method [24]. The distribution model, also known as the maximum posterior probability (MAP), defines the updated model, and the distribution around the model provides a measure of confidence in the model. The peak distribution often means a high degree of confidence in the updated model, and shallowness may make uncertain the choice of parameters [25,26]. Bayesian updating and prediction methods are widely used in geotechnical engineering, such as pile bearing capacity analysis, soft soil consolidation inversion analysis, jack-up drilling platform bearing capacity analysis, foundation pit support under normal use protection, and excavation problems and slope stability problems [27]. The key difficulty of the Bayesian method lies in the calculation of the MAP distribution. Because of the complexity of the MAP distribution and the possibility of involving high-dimensional integration problems, the direct numerical integration method is less efficient. To overcome this problem, the MCMC method avoids the computational complexity of directly solving the MAP distribution by generating a large number of random samples that obey the MAP distribution, and then it uses these random samples to characterize the MAP distribution [28,29]. Metropolis–Hastings (M–H) algorithm, as one of the commonly used MCMC methods, has been widely used in geotechnical engineering [30]. In this paper, the M–H method is used to generate a large number of MCMC random samples that follow the MAP distribution parameters, such as the posterior mean, posterior standard deviation, and posterior distribution. These statistics are the basic inputs for reliability analysis of geotechnical engineering [31].

Therefore, the purpose of this article is to show how to use MCMC simulations to effectively describe SWCC parameter uncertainty involving random variables. The structure of this article is as follows. First, the theory of the SWCC and the Bayesian theory of parameter uncertainty are introduced. Second, the hybrid Markov chain with emphasis on the uncertainty characterization of

the geotechnical model is introduced, and a framework for updating the posterior distribution of the parameters of SWCC is introduced. Then, a large-scale landslide model experiment is used to illustrate how to obtain the posterior distribution of the parameters of the SWCC through Bayesian updating. Finally, we analyze the influencing factors of posterior distribution and characterize the uncertainty of SWCC parameters.

#### **2. Theories and Methods**

#### *2.1. Soil–Water Characteristic Curve*

Among the existing empirical equations for curve-fitting of the SWCC, the most commonly used VG model was adopted in the present study, which is given by [12]:

$$
\theta\_w = \theta\_r + \frac{\theta\_s - \theta\_r}{[1 + \frac{\lambda}{a}^n]^{1 - \frac{1}{n}}} \tag{1}
$$

where *θ<sup>w</sup>* is the volumetric water content; *θ<sup>r</sup>* is the residual volumetric water content; *θ<sup>s</sup>* is the saturated volumetric water content; *λ* is matric suction; *a* and *n* are non-negative curve fitting parameters.

To obtain the prior distribution of SWCC parameters, a simple sensitivity analysis of the modeling parameters sheds light on the possible range of the numeric values of the parameters [23]. The morphology of the SWCC between *a* from 0.1 to 10 and *n* from 1.1 to 5 was studied. Figure 1a describes the influence of the value of *a* on the SWCC model when *n* is 5. Figure 1b describes the effect of the value of *n* on the SWCC model when *a* is 2. Each model parameter has a certain physical meaning, except *a* and *n*. From Figure 1, it becomes clear that the parameters a and n are critical for SWCC, determining the shape of the curve.

**Figure 1.** (**a**) The influence of parameters *a* on the soil–water characteristic curve (SWCC), and the *n* is 5 in this figure; (**b**) The influences of parameters *n* on the SWCC equation, and the *a* is 2 in this figure.

#### *2.2. Bayesian Model Updating*

The principle of Bayesian model updating is the Bayes' theorem. The pioneering work in the late 1990s attracted the attention of the engineering research community and gave detailed and up-to-date instructions [32,33]. Then, the obstacle to computing in practical engineering applications was the high level of computer resources required. This problem still exists, but there has been great progress in solving it.

The Bayesian model's updating process is based on Bayes' theorem. In our case, likelihood is a function of the parameters in the statistical model, which represents the likelihood of model parameters. The prior distribution is based on existing research assumptions. The posterior distribution can be

calculated by combining the prior distribution and the likelihood function. Its general formula is given by:

$$P(\theta|\mathcal{X}\_{\text{cbs}}) = \frac{P(\mathcal{X}\_{\text{cbs}}|\theta)P(\theta)}{P(\mathcal{X}\_{\text{cbs}})} \tag{2}$$

where *P*(*θ*) is the prior distribution representing the initial knowledge about the parameters *θ*; *P*(*θ*|**X**obs) is the posterior distribution representing the updated knowledge based on the observation data; *P*(**X**obs) is the normalization factor; *P*(**X**obs|*θ*) is the likelihood function of **X**obs for an instance of the parameters *θ*.

If it is possible to provide an update a prior distribution according the available knowledge, then the MAP distribution can be obtained. The Bayesian updating is a challenging problem, since it is very difficult to directly calculate the posterior probability density over the entire parameter space. However, the well-known M–H algorithm is an effective update tool [15]. The algorithm is essentially an iterative method, sampling from a series of intermediate partial differential equations that gradually converge to the MAP distribution. The *j*th intermediate probability density function (PDF) is expressed as a function of the likelihood probability function *PL* as:

$$P\_{\rangle} = P\_L(\mathfrak{X}\_{\text{obs}} | \theta)^{\beta\_{\uparrow}} P(\theta) \tag{3}$$

where the index of likelihood *β<sup>j</sup>* is the so-called reduction factor. Its value gradually increases from *βo*= 0 in the first iteration until it reaches *β<sup>j</sup>* = 1 in the last step [24]. *β<sup>j</sup>* is adaptively calculated based on the sample from the previous step.

The likelihood is a key component of the Bayesian update framework because it quantifies how relevant the model is to a given parameter instance by representing the likelihood of observation. Under the assumption of independence between observations, the probability *PL* is theoretically defined as:

$$P\_L(\mathbf{X}\_{\text{obs}}|\theta) = \prod\_{k=1}^{N\_{\text{obs}}} P(\mathbf{X}\_k|\theta) \tag{4}$$

where *P*(**X***k*|*θ*) is the probability density value at **X***<sup>k</sup>* for a given set of parameters *θ*; **X***<sup>k</sup>* is the *k*th **X**obs. *θ* is a set of parameters—in this analysis, the parameters *a* and *n*. In this process, the PDF should be estimated separately for each instance. However, in a complete Bayesian updating process, the required computation is costly, due to the large number of evaluations that are necessary to obtain to define a statistical sample.

#### **3. The Bayesian Updating Procedure of Monitor Data in the Large-Scale Landslide Model Experiment**

In this paper, the above-mentioned Bayesian updating method is used to calculate the SWCC model parameters. Bayesian updating is a method of statistical inference, in which as more evidence or information (in this paper is the observation samples) is obtained, Bayes' theorem (Equation (2)) is used to update the probability of a hypothesis. First, the prior distribution of the SWCC model parameters is determined through reference research data or data from a SWCC experiment, as done in this work. Then, one uses the M–H method to generate a defined (large) number of random MCMC samples following a iterative calculation until *β<sup>j</sup>* = 1, and the posterior distribution is obtained. The M–H method is a powerful Markov chain method with which to simulate multivariate distributions. The main steps include initialization and iteration [34]. Figure 2 shows a flowchart for implementation of the proposed Bayesian approach schematically. The work done covers both the experimental results to feed the model (left branch on Figure 2), and the analysis methodology (right branch on Figure 2). The experimental data were obtained after basic soil mechanical experiments, to characterize the material of the landslide body, and a large-scale landslide model experiment, including water content sensors to measure in real-time the hydrodynamic process. The setup and results of the experiments are described in Sections 4 and 5. Each step and its detailed description are summarized as follows.

	- (a) Infiltration experiments and a triaxial experiment were carried out, and the material ratio of landslide body was determined to be 7:3 for soil:stone through these experiments.
	- (b) Large-scale landslide models were built according to the layer-by-layer compaction method. The place where the sensor was buried was carefully compacted to avoid damage to the sensor.
	- (c) After the landslide model and sensor were ready, the moisture content sensor was tested to determine the reliability of its data.
	- (d) The SWCC model was selected (VG model), and we analyzed the two model parameters within its range to observe their influences on the curve shape.
	- (e) We assumed the prior distribution as normal or lognormal according to existing research. Then we determined the prior distribution of SWCC for the model parameters *a* and *n*, based on experimental data and research results.
	- (f) A finite element (FE) numerical model was developed according to the large-scale landslide test; we evaluated the seepage flow under the conditions of the reservoir water level.
	- (g) In the Bayesian updating process, it is necessary to call the FE model, taking a long time to obtain the large number of output values. Therefore, this paper used an artificial neural network (ANN) as a surrogate model to speed up the analysis.
	- (h) The MCMC algorithm was used to generate *N* random samples that followed the posterior distribution of SWCC parameters and iteratively continued until *β<sup>j</sup>* = 1. We got the posterior distribution of SWCC model parameters (posterior most likely value, posterior mean, posterior standard deviation, etc.).
	- (i) We analyzed the factors influencing the posterior distribution results, such as the number of observation datums, the prior distribution, and the observation data that changed over time.
	- (j) We characterized the uncertainty of the SWCC model. For the parameters of the SWCC, the *N* samples generated from the MCMC algorithm were ranked, and a 95% confidence interval was obtained.

These steps were programmed, and MATLAB 2019b was used in this process. The programming of the program is divided into three steps: firstly, import the simulated and observed values; then use the MCMC algorithm to generate a large number of samples (100,000 sets of data in this program) to iterate for Bayesian updating; finally, the posterior distribution of the parameters of the SWCC model is visualized. In a different analysis case, the same analysis tools can be used so that geotechnical practitioners only need to provide laboratory or in-situ test results, prior knowledge (e.g., reasonable ranges of model parameters), and candidate SWCC models as inputs. The main advantage of the above method is that it can quickly check distributions of the posterior parameters of the SWCC.

**Figure 2.** Flowchart of the Bayesian updating procedure.

#### **4. The Large-Scale Landslide Model Experiment and Bayesian Updating**

A large-scale landslide experiment was conducted using a soil–rock mixture and by monitoring the changes in water content to clarify the water migration during a landslide. The principle of numerical modeling, using an earlier developed model, is briefly explained. Simulations are presented and their results are compared to the experimental ones.

#### *4.1. Experiment Setup*

The experiments have been carried out in a large-scale landslide model experiment at The Key Laboratory of Geological Hazards and the Ministry of Education of the Three Gorges Reservoir Area of the Three Gorges University (China); the setup is shown in Figure 3. The model test platform was 6 m long, 0.8 m wide, and 1.5 m high. The door at the right end of the test platform could control the rise and fall of the reservoir water level. The landslide model was connected to external water supply pipes and drainage pipes. The speed at which the water level of the reservoir changed was adjusted by setting a flow valve.

**Figure 3.** The side view of the landslide model and configuration of instruments in the experiment. A drawing of the configuration of the experimental model, showing the bedrock region (within blue lines), the landslide body region (below the green line), and the embedded sensors to monitor the water content.

According to the nature of the landslide body in the Three Gorges reservoir area, the material of the landslide body was designed. The material used in the landslide body was composed of soil and small stones, and the ratio was 7:3. The landslide model test was completed according to the principle of landslide model forming. Considering the size of the model and the need to embed sensors, it was decided to choose the tamping method to build the landslide model. The sensors were embedded in the model in advance. The moisture content sensors were arranged in four sections, with 3 sensors in each section. Then we increased the water level of the reservoir to monitor the change of its water content. The overall landslide model is shown in Figure 4. In this figure, instruments for displacement monitoring and high-density electrical methods are included (the results are not used in this article.) In order to obtain the material parameters of the landslide model, a triaxial test and a constant head permeability coefficient experiment were carried out. A pressure plate instrument was used to obtain the characteristic curve of unsaturated soil water, as shown in Figure 5a,b. In this figure, the pressure pump and pressure plate instruments are displayed.

**Figure 4.** A picture of the landslide model built indoors, from a frontal view: the model was 1.5 m (height) ∗ 0.8 m (width). The water content sensor is inside the landslide body (not visible; the actual water content sensor is shown in Figure 6a).

**Figure 5.** Pressure plate instrument: (**a**) pump; (**b**) pressure gauge.

Sensors were used to obtain water content data in real time. Figure 6a shows the water content sensor. The power supply voltage of the sensor was 5 V, and the output signal type was RS485. The monitoring data of the water content sensor were saved when the water level was from 0.3 to 1 m. Figure 6b shows the sensor data interface. Its function was to convert RS485 signals into USB signals. Based on the Modbus protocol, the data acquisition software was developed using C# to save the water content data to the database. All sensors conducted data sampling and recording every 10 s. Figure 7 is a histogram representation of moisture sensor monitoring data. In this figure, the abscissa represents different moisture content, and the ordinate represents its frequency. These water content data were used as observation samples in Bayesian updating (specifically, as shown in Figure 2).

**Figure 6.** (**a**) Water content sensor (the rated voltage was 5 V); (**b**) signal conversion device (converting the RS485 signal to USB signal).

#### *4.2. FEM Simulation*

In the Bayesian updating process, a numerical model of the experimental setup is needed to provide simulation results. The SEEP/W software, a part of Geo Studio is a numerical simulation software for seepage analysis. In this study, Geo Studio software (2007) was used. It can mathematically simulate the real physical process of water flowing through the soil medium based on the finite element model. The equation used in the process is a steady-state seepage equation. The work of this section aims to model the seepage analysis of the landslide model in Figure 3 through the FEM. For this reason, a geometric model of a large-scale landslide model was built in FEM software. The purpose of this research was to study the seepage characteristics of landslides through FEM to simulate the changes in water content of landslides when the water level rises. We used SEEP/W software to develop and calibrate computer models of landslide model, and compare observed data with simulated data.

**Figure 7.** Monitoring data of water content.

In the FE model process, the main steps are to establish a landslide model, set material parameters, divide calculation units, and set boundary conditions [35]. Firstly, SEEP/W was used for geometric modeling of the landslide model according to the dimensions in Figure 3. The material models and properties of the soil were put into saturated-unsaturated seepage model, defined with a seepage constitutive relationship requiring two functions: the permeability coefficient function and water content function (SWCC) [35]. In this process, the saturated permeability coefficient was set to 1.6 × <sup>10</sup>6, and the SWCC parameters were set according to the prior distribution. Figure 8 shows the geometric model and the mesh used in this work. In this figure, the yellow area represents the landslide body. The unit accuracy of the finite element mesh was set to 0.1 m. The red part represents the boundary conditions. Two boundary conditions were set for the numerical model based on water level changes: water level boundary conditions and pressure boundary conditions. The boundary conditions of the model were a key component of the numerical analysis, and the ability to control the boundary conditions was also a direct reason for the powerful numerical analysis. In this process, the calculation method was the steady-state calculation, and the solver adopted the implicit solution method.

The FEM can evaluate water content values, depending on the input parameters (a,n) of the SWCC model. This results were the simulation samples used in Bayesian updating. The main purpose of this was to obtain the complete results of the FE model evaluation by SEEP/W software calculation, including water pressure, moisture content, etc. The results calculated by the SEEP/W software are analyzed and compared with the monitoring data of the landslide model in Section 4.1 to check for errors. The calculation of the SEEP/W software in this section also provides a model for the establishment of the surrogate model in the next section.

**Figure 8.** Mesh generation of finite element model (FEM).

#### *4.3. Surrogate Model*

Bayesian updating requires two sets of data: one of experimental data, and the other of simulation data. The experimental data were obtained through the above-mentioned large-scale landslide model experiment, and the data are shown in Figure 7. The simulation data were calculated according to the FEM in Section 4.2. It took 690 h (about 29 days) to perform a complete Bayesian update calculation with the FEM, and 2 h to use the surrogate model (Intel(R) Core(TM) i5-5200U CPU @ 2.20 GHz, RAM 8.00 GB). Therefore, a surrogate model was used in this study to obtain the simulation value. The surrogate model is a commonly used optimization method in engineering problems. When the actual problem is very computationally intensive and difficult to solve, a simplified model with less computational complexity and a rapid solution can be used to replace the original model to speed up the optimization process. In this study, an ANN was used to build the surrogate model.

In order to improve the efficiency of Bayesian updating, we adopted the method of the surrogate model. We used an ANN to train the surrogate model. The ANN is an algorithm whose design is driven by the function of the human brain and its components. ANNs use nonlinear techniques, so they can be used to model highly complex and nonlinear systems. In this study, a fully connected feedforward multi-layer structure using a backpropagation momentum learning algorithm was adopted. This kind of ANN architecture usually consists of an input layer, some hidden layers, and an output layer. For us there were 10 neuron nodes used in each layer. Each node was interconnected to the nodes of the upper layer through adjustable weights, and horizontal, self-direction, or backward connection was not allowed. The calculation and expression ability of a single neuron is limited. However, when they are connected to each other, the entire network can model complex functions. The training pattern consists of a set of matched input and output vectors. The logical sigmoid function was used as the activation function. In order to develop the artificial neural network model, Neural Network Fitting toolbox and MATLAB R2019b were used. The Neural Network toolbox provides various parameters for neural network development, which can be selected flexibly. We set the distribution of the parameters of the SWCC model, and then used the SEEP/W software to calculate the value of the water content. We saved the SWCC model parameters and the generated water content data for the establishment of the surrogate model. The input data selected for this part were a array (2 × 1000) of SWCC parameters *a* and *n*, and the target array was water content data. In the training process, 70% of the data were used for training, 15% of the data for verification, and 15% of the data for testing. The training data were presented to the network during the training process, and the network was adjusted according to its errors. The verification data were used to measure the generalization of the network and stop training when the generalization stopped improving. The test data had no effect on training, so they provide an independent measurement of network performance during and after training. In this section, the number of hidden neurons was 10. Finally, we checked the error of the surrogate model built using an ANN.

Figure 9 shows the regression diagram of the water content and SWCC parameters at specific points of the landslide model. In this figure, we can observe that all datasets were correctly fitted to the row. This shows that our neural network structure was correct, and it can also be used to predict the outputs of other input datasets. Figure 9a is the regression diagram of 70% of the training data. Figure 9b is the regression diagram of the 15% group of the validation data. Figure 9c is the regression diagram of the 15% group of test data. Figure 9d is the regression diagram of all data. The results of ANN prediction quality are detailed in Table 1. The goodness of fit (R) is close to 1 and the expected value of the mean square of the error (MSE) is small. The model results show that the ANN model can be used to predict the relationship between the water content at specific points of the landslide model and the parameters of the SWCC.

**Figure 9.** Neural network training regression: (**a**) The regression diagram of the 70% group of the training data; (**b**) The regression diagram of the 15% group of the validation data; (**c**) The regression diagram of the 15% group of test data; (**d**) The regression diagram of all data.

**Table 1.** The results of the prediction with the artificial neural networks (ANN).


#### **5. Results**

#### *5.1. Effect of Dataset Size on the Outcome of Bayesian Updating*

Using the surrogate model and the experimental data, the iterative Bayesian updating was done until the simulated value was close to the observed value. Bayesian updating is an iterative process. First, it sets the prior distribution of SWCC parameters and uses the MCMC algorithm to generate a large number of samples that obey the prior distribution. Then it performs iterative calculations, using surrogate models and experimental water cut monitoring data. Bayesian updating is performed until the simulated value is close to the observed value. For each iteration, a set of samples of SWCC parameters are generated. All the data of the iteration are presented into a plotmatrix, as shown in Figure 10a. In this figure, the top-left figure is the frequency diagram of parameter *a*, the bottom-right figure is the frequency diagram of parameter *n*, and the bottom-left and top-right figures represent the correlation. Figure 10b shows the distribution of SWCC parameters obtained by Bayesian updating, converged after 5 iterations (*j* = 5). This is the plotmatrix of the SWCC parameters after the fifth iteration.

**Figure 10.** Results of using 1000 sets of water content monitoring data to update: (**a**) prior distribution; (**b**) posterior distribution.

In order to study the changes in the posterior distribution when monitoring the increase in the amount of data, we used 200, 500, 1000, 2000, and 5000 sets of monitoring data to calculate the posterior distribution and analyze the results. As the amount of data increased, the trends of Bayesian updating results are shown in Figure 11, and the values of Bayesian updating results are shown in Table 2. In that figure, the red line represents the median, the upper and lower lines of the blue rectangle represent the upper and lower quartiles, the short red lines represent outliers, and the short black lines represent the minimum and maximum values. The trend of the change of parameter *a* is shown in Figure 11a. It can be seen that the median value of parameter *a* fluctuates around 1.46 and has a tendency to increase slightly. Outliers are distributed below the minimum value, indicating that the posterior distribution is left skewed. The trend of the change of parameter *n* is shown in Figure 11b. It can be seen that the median value of the parameter *n* fluctuated around 3.98. Outliers are distributed over the maximum value, indicating that the posterior distribution skewed right. When the number of data was 1000, the data of the posterior distribution were relatively stable. Therefore, in the following analysis, 1000 sets of observation data were used as samples to study the influence of the distribution type on the posterior distribution.

**Figure 11.** Impact of data volume on update results: (**a**) Parameter *a*; (**b**) parameter *n*.


**Table 2.** Impact of dataset size on update results.

#### *5.2. The Effect of Prior Distribution on the Outcome of Bayesian Updating*

In this Bayesian updating process, the prior distribution of SWCC parameters *a* and *n* needs to be set first. In the previous analysis (Section 5.1) the prior distribution of variables a and n was adjusted to a lognormal distribution. The converged posterior distribution was analyzed and the most likely value, mean, and standard deviation of the parameter were evaluated. In Table 3, the results for the analysis with dataset size value 1000 are shown.

To investigate the effect that a different prior distribution has on the results, the prior distribution was changed to a normal distribution. The results are plotted in Figure 12, and the distribution parameters are shown in Table 3. Comparing the results of Table 3, the values of the parameters *a* and *n* are compatible within the deviations obtained.

**Figure 12.** Update results when the prior distribution is a normal distribution: (**a**) prior distribution; (**b**) posterior distribution.

**Table 3.** Impact of prior distribution on update results.


#### *5.3. The Effects of Six Stages on the Outcome of Bayesian Updating*

During the experiment (see Section 4.1 for details), we set the data collection time interval of the moisture content sensor to 10 s. When the water level changed from 0.3 m to 1 m, six stages were selected to study the trend of the posterior distribution of SWCC parameters when the water level changed. With the monitoring data obtained over time, the trend of Bayesian updating results is shown in Figure 13, and the value of Bayesian updating results is shown in Table 4. In this boxplot, the red line represents the median, the upper and lower lines of the blue rectangle represent the upper

and lower quartiles, the red short lines represent the outliers, and the black short lines represent the minimum and maximum values. The trend of the change of parameter *a* is shown in Figure 13a. It can be seen that the median value of parameter *a* fluctuates around 1.47 and has a tendency to increase slightly. Outliers are distributed above the maximum value. The trend of the change of parameter *n* is shown in Figure 13b. It can be seen that the median value of parameter *n* fluctuates around 3.98 and has a slightly lower trend.

**Figure 13.** The influence of the monitoring value at different time stages on the update result: (**a**) parameter *a*; (**b**) parameter *n*.


**Table 4.** Impacts of different stages on update results.

#### *5.4. Characterization of Uncertainty*

According to the previous results, it is possible to characterize the uncertainty of the model parameters. For a sufficiently large number of samples, and considering the minor effects caused by the selected prior distribution type and the time stage, the uncertainties obtained are reliable. Considering the VG model for the SWCC, the parameters *a* and *n* were evaluated with 100,000 MCMC samples and 95% confidence intervals (for 97.5% and 2.5% upper and lower bounds). The results are plotted on the curves in Figure 14. As can be seen in this figure, the most likely SWCC models were all within the 95% confidence interval. Therefore, the proposed method can reasonably determine the SWCC model based on limited experimental data and quantify the uncertainty of the SWCC model estimate.

**Figure 14.** The 95% confidence interval of VG obtained from the Bayesian updating.

#### **6. Discussion**

Although slight changes of SWCC parameters have minor effects on the stability of a landslide, they improve the accuracy of landslide reliability calculation. It is very necessary to conduct in-depth research on SWCC parameters. In previous studies, the main focus was on the SWCC model. We used indoor test data to select SWCC model [36]. Compared with the deterministic SWCC model, this paper took into account the uncertain factors caused by soil variability and unstable monitoring data. This paper mainly checked the parameters of SWCC based on the monitoring data of large-scale model tests.

This paper mainly describes the impacts of monitoring data volume, prior distribution type, and different stages on posterior distribution. Figure 11 and Table 2 show the effect of data volume from 200 to 5000 on the posterior distribution of SWCC model parameters. The left figure represents the prior distribution, and it can be clearly seen that the prior distribution is the same. The figure on the right represents the posterior distribution. As the number of observation samples increases, it has little effect on the posterior distribution of SWCC parameters. It can be seen from Table 2 that the results of the posterior distribution are very stable. Therefore, it is feasible to perform Bayesian updating of parameters using monitoring data. This makes it possible to realize the online SWCC parameter monitoring and analysis system. In the research, this study selected 1000 sets of data for Bayesian updating. Figures 10 and 12 compare the effects of the prior distribution type on the results of the SWCC parameters' posterior distributions. It can be clearly seen that the coefficient of variation is smaller when the prior distribution is lognormal. Therefore, when using monitoring data to perform Bayesian updating of SWCC parameters, it is more reasonable to set the prior distribution to the lognormal distribution to obtain the posterior distribution. Figure 13 and Table 4 show the changes in the posterior distributions of SWCC parameters with different stages. In this figure, it can be seen that the mean value of parameter *a* in the SWCC posterior distribution increases from 1.462 to 1.480, while the mean value of parameter *n* fluctuates from 3.985 to 3.983 (the fluctuation range is within 0.1). Due to the short monitoring time, the change to the posterior distribution is very subtle. In the follow-up research, the monitoring time could be extended to obtain the SWCC parameter change rule. Figure 14 shows the 95% confidence interval of the posterior distribution of SWCC. The most likely SWCC falls within the 95% confidence interval of its posterior distribution.

This paper studied the method of using landslide monitoring data to study SWCC parameters. However, due to the influences of many uncertain factors, the simulation model's prediction and the real structural response often have certain errors. Therefore, it is necessary to verify the prediction accuracy of the simulation model in future research.

#### **7. Conclusions**

This paper presents a SWCC parameter calculation framework based on Bayesian landslide monitoring data theory. This method makes full use of prior information and limited monitoring data to determine the posterior information of SWCC model parameters. It characterizes the uncertainty of the SWCC model and is of great significance for reliability analysis. The influences of the dataset size, prior distribution type, and time-varying monitoring data on the result of posterior distribution were studied. Finally, the uncertainty of the SWCC model parameters was characterized. Taking the water content monitoring data of the large-scale landslide model test as an example, the effectiveness of this method is explained. The main conclusions from this study are as follows:


**Author Contributions:** Conceptualization, C.F. and M.B. (Michael Beer); data collection, C.F., X.L. and T.H.; methodology, C.F., M.B. (Matteo Broggi), and B.T.; software, C.F. and M.B. (Matteo Broggi); validation, C.F., B.X., and B.T.; data curation, C.F.; writing—original draft preparation, C.F.; writing—review and editing, S.B.; project administration, B.T.; funding acquisition, B.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the National Key Research and Development Program of China (2017YFC1501100), and the Research Fund for Excellent Dissertation of China Three Gorges University (2019SSPY005).

**Acknowledgments:** Thanks to the National Research Council and Three Gorges University for their support.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **3D Dilatometer Time-Series Analysis for a Better Understanding of the Dynamics of a Giant Slow-Moving Landslide**

### **Jan Blah ˚ut 1,\*, Jan Balek 1, Michal Eliaš <sup>2</sup> and Stavros Meletlidis <sup>3</sup>**


Received: 30 June 2020; Accepted: 3 August 2020; Published: 7 August 2020

#### **Featured Application: Analysis of monitoring data from very slow-moving landslides.**

**Abstract:** This paper presents a methodological approach to the time-series analysis of movement monitoring data of a large slow-moving landslide. It combines different methods of data manipulation to decrease the subjectivity of a researcher and provides a fully quantitative approach for analyzing large amounts of data. The methodology was applied to 3D dilatometric data acquired from the giant San Andrés Landslide on El Hierro in the Canary Islands in the period from October 2013 to April 2019. The landslide is a creeping volcanic flank collapse showing a decrease of speed of movement during the monitoring period. Despite the fact that clear and unambiguous geological interpretations cannot be made, the analysis is capable of showing correlations of the changes of the movement with increased seismicity and, to some point, with precipitation. We consider this methodology being the first step in automatizing and increasing the objectivity of analysis of slow-moving landslide monitoring data.

**Keywords:** slow-moving landslide; landslide monitoring; time-series analysis; San Andrés Landslide; El Hierro; Canary Islands

#### **1. Introduction**

Analysis of monitoring data from very slowly moving landslides faces very difficult issues considering the precision and accuracy of the used monitoring instruments, the influence of the different environmental factors on the measured data, the frequency of the measurements and the overall complexity of the phenomena [1–4]. In the case of large landslides (of more than several km2) another problem arises with choosing the right placement of the monitoring devices on the ground, as remote sensing methods such as LiDAR or InSAR, despite their enormous applicability, still do not provide comparable sub-mm accuracy [5,6].

When analyzing the landslide movement monitoring results, the time-series analysis should be used. Time-series analysis comprises different statistical methods to understand the underlying context of data points or to make forecasts [7]. It can capture seasonal behaviors, trends and changes [8]. Time-series analysis is routinely applied in econometrics [9,10]. In geosciences, there is a constant use of time-series analysis within GNSS observations [11] or climatology [12]. Thus, it might be a bit surprising that in landslide science, time-series analysis studies are mostly related to remote sensing

(SAR) data analysis [13–15]. To some point, this is not surprising as SAR data are available in defined time-steps and form an undisrupted time series, which is the basic assumption of the analysis. In the case of landslides studies, time-series analysis is usually focused on trend behavior, sudden changes and seasonality analysis [16–18]. The big advantage of time-series analysis methods lies in the possibility to make predictions of future behavior based on previous activity. This is especially important in landslide disaster management. Complex non-linear time series are ubiquitous in geosciences [19]. The recently developed adaptive regression models (used in this work) appear to offer good solutions in multivariate non-stationary time-series analyses [20].

In the era of Big Data, it is more and more possible to harvest/gather a large amount of environmental and geophysical data related to slope stability using a variety of sensors and systems with high precision [21,22], which can be used for interpreting and monitoring the results and finding correlations to understand the landslide activity. However, a recent study [23] highlighted that the observations are prohibitively site-specific, highly subjective and fail to account for the spatial variability and correlations in slope movements and common trigger factors.

The landslide movement is usually neither smooth nor continuous. There are changes in the nature of the observed monitoring data at certain moments (changes of speed or direction of movement, jump events). We assume that these changes (events) during the development of landslide movements are caused by extreme states of factors (triggers) affecting the landslide movement (i.e., seismicity, precipitation) [24].

In this paper, we analyzed data from precise 3D dilatometers installed on the detachment planes of a giant complex landslide in a volcanic environment. The time-series analysis was performed on 1D, 2D and 3D data projections using adaptive regression. The main aim of this work was to find a suitable methodology for evaluating the time series of micro-movements to find the changes in landslide movement (events) in them and then compare them with the time moments of extreme external conditions including available climatic and seismic data to find correlations with the possible triggers of the landslide activity. The aim was fulfilled by applying known methods and combining them to understand better the landslide behavior and its responses to possible influencing factors.

#### **2. Study Area**

This study takes place on the island of El Hierro, Canary Islands, Spain. El Hierro is a composed oceanic volcano with the oldest subaerial rocks having an age of 1.12 Ma. Its youngest subaerially exposed rocks are represented by the Rift Series, with a maximum age of 0.16 Ma [25]. Over the past 33,000 years, onshore eruptions have reoccurred approximately once every 1000 years [26]. The latest, ongoing phase of volcanism began around 2.5 ka [27]. Recently a period of intense seismic activity spanned from July 2011 [28] to 2014 [29]. During this period, an offshore eruption started on 10 October 2011 and finished in March 2012 [30]. A more detailed description of the geology of the island has recently been presented elsewhere [31,32].

The characteristic three-point star morphology of El Hierro is a result of several enormous flank collapses (Figure 1). Until now, seven debris avalanches have been identified: Tiñor (<880 ka), Las Playas I (545–176 ka), Las Playas II (176–145 ka), El Julan (>158 ka), El Golfo A (176–133 ka), El Golfo B (87–39 ka), and Punta del Norte (unknown age) [27,33–42].

In addition to the previously mentioned flank collapses, a giant San Andrés Landslide (SAL) has been identified on the northeastern part of the island. The SAL is a deep-seated gravitational slope deformation [43] or a slump [44]. It has been supposed that the landslide underwent a single-event development at some point between 545 ka and 176 ka [45] and was inactive since then. Recent research shows, however, that there were at least two distinguishable separate slip events. The first one between 545 ka and 430 ka, and the second one between 183 ka and 52 ka [46]. Moreover, it has been proven that there is a continuous creep with a rate up to 0.5 mm·a−1, suggesting that the landslide mass may be moving steadily to the east and southeast [32,47,48]. Numerical stability modelling [49] showed that creep movements are possible due to unconsolidated sediments at the sea bottom, and destabilization of the landslide might be possible during periods of intense seismic activity.

**Figure 1.** Location of the study area: (**a**) the El Hierro Island with a simplified geological map showing the location of the TM-71 3D dilatometers and showing the direction of the axis (X, Y, Z) of the measurement: grey rectangle shows the location of Figure 8c; (**b**) dilatometer HIE2 located underground on the upper scarp of the SAL; (**c**) dilatometer HIE3 located in the open air on the left side of the SAL.

This creep has been monitored by high-precision (±0.007 mm) 3D dilatometer gauges [3] with automatic data processing [50]. There are two devices installed on the landslide detachment plane in the upper scarp and on the left side of the landslide, respectively (Figure 1). The instrument records data based on optical−mechanical interferometry via the generation of moiré patterns, which result from the bending and interference of light rays as they pass through two specially designed optical grids [51].

#### **3. Monitoring Data Collection and Preparation**

The data used in this analysis can be split into three categories: (i) data from the 3D dilatometers placed on the SAL detachment plane; (ii) climatic variables and (iii) seismicity (Figure 2).

**Figure 2.** Data used in the analysis: (**a**) HIE2 dilatometer: *X*-axis—black line, *Y*-axis—red line, *Z*-axis—blue line; (**b**) HIE3 dilatometer: *X*-axis—black line, *Y*-axis—red line, *Z*-axis—blue line; (**c**) daily average air temperatures; (**d**) daily precipitation; (**e**) registered earthquakes; (**f**) calculated PGA.

#### *3.1. 3D Dilatometer Data*

Two 3D dilatometers named HIE2 and HIE3 are placed on the SAL's detachment plane. They are about 2.5 km apart (Figure 1) and the measurements are taken every 24 h at 00:00 UTC. The records (Figure 2) are analyzed in a Cartesian coordinate system with *X*-axis showing compression/extension perpendicular to the detachment plane; *Y*-axis showing strike-slip movement along the detachment plane and *Z*-axis showing vertical movement perpendicular to the detachment plane. As the records are taken at 00:00 UTC the records were correlated with climatic and geophysical variables from the previous day. Automated monitoring started on 16/10/2013 and covers a period of 2010 d till 17/04/2019. Due to technical difficulties, several instrument failures resulting in missing records happened. In the case of HIE2 device, there are 494 missing measurements (24.6%), in the case of HIE3 there are 358 missing records (17.8%).

#### *3.2. Climatic Variables*

Climatic variables are strongly influencing any landslide behavior. For that purpose, temperature, precipitation and insolation data from the El Hierro airport [52] were downloaded. The location of this station is not directly on the SAL. However, it is in its vicinity and provides daily undisturbed data for the whole monitoring period. Precipitation can play a principal role in landslide activation. For that purpose, the daily rainfall amounts were analyzed. Additionally, the 30-day antecedent precipitation index (API30) [53,54] was calculated to better depict the influence of retrospective precipitation history and its influence on soil moisture [55]:

$$API\_n = \sum\_{i=0}^{n} \left( c^i \cdot P\_{n-i} \right) \tag{1}$$

where *n* is the number of antecedent days; *c* is decay constant and it depends on watershed and seasonal parameters. The value of *c* ranges between 0.80 and 0.98 [56]. In this study, a *c* of 0.95 was used; *i* is the number of days counting backwards from the date on which the API is determined and *Pn* is the amount of precipitation *i* days before the causal rainfall (mm).

Mean and maximum daily temperatures can influence both the monitoring instrument and the rock on which it is attached by thermal dilation [57]. For that purpose, the temperature data were used to observe influence on the monitoring records and assess seasonal trends.

#### *3.3. Seismic Data*

Seismicity plays a major role in the stability of volcanic flanks. During the analyzed period, a total of 6537 seismic events were recorded around El Hierro Island by the Instituto Geográfico Nacional network [58]. We used this data to calculate the PGAs (Peak Ground Accelerations) of these events. The Canary Islands does not have developed local ground-motion attenuation relationship and there are no accelerometer network data available. For that reason, we applied an attenuation relationship developed for the Hawaiian Islands (to some point, a similar volcanic archipelago) by [59] and suggested to be used on the Canary Islands by [60]:

$$
\log\_{10}{PGA} = 0.518 + 0.387 \left( M - 6 \right) - \log\_{10}{r} - 0.00256r + 0.335S \tag{2}
$$

where *PGA* is measured in units of g, *M* is the magnitude, *r* = (*d*<sup>2</sup> + 11.292)1/<sup>2</sup> and *S* is 0 for lava sites and 1 at ash sites. The distance parameter *d* is the closest distance from the recording site to the surface projection of the fault rupture area.

For this analysis, we used the maximum recorded PGA for the particular day and also calculated the sums of daily PGAs (named PGA dose). The maximum PGAs and PGA doses were calculated for both HIE2 and HIE3 instruments considering their geographic location.

#### **4. Methodological Approach**

#### *4.1. Input Data Pre-Processing*

The relative motion on the detachment planes of the SAL is represented by time series which are not completely homogeneous. There are many missing values due to technical problems. This is not necessarily a problem as the reading frequency of the device (every 24 h) to the speed of the movement is very high. However, the methodology of time-series processing needs to be adapted to this fact.

The time series were decomposed into trend, residual and seasonal components (as is visible on Figure 3). The periodicity of the seasonal components reflects the effect of the thermal dilation of the anchoring elements holding the TM-71 device, but it can also be partially influenced by the volume changes of the monitored blocks. A Fourier model was fitted to estimate the characteristics of the periodic component, specifically the 1st harmonics. The period of the harmonic function was fixed at 365 d (Figure 3b). The applicability of the harmonic function 0.15 mm corresponds to the amplitude of the air temperature (4 ◦C) and the length of the steel anchoring elements (2 m) at the coefficient of thermal expansion <sup>α</sup> <sup>=</sup> 16.3 <sup>×</sup> <sup>10</sup><sup>−</sup>6.

**Figure 3.** Example of the time-series decomposition: (**a**) close to a linear decrease of the velocity on HIE3 *Y*-axis; (**b**) significant seasonal periodicity visible on HIE2 *X*-axis.

The trend component is very well described by the approximation by a second-order polynomial. This indicates a linear change in velocity over time. This trend is clearly visible on the HIE3 *X*-axis (Figure 3).

#### *4.2. Detection of Change-Points in Time Series of Individual Components*

The time series of the measured movements are not stationary but often express sudden changes. We assume that these changes are caused by external factors. Over time, several methods have been developed to detect such changes such as CUSUM (Cumulative Sum) and its modifications (firstly introduced by [61]). Limiting factors of a wide range of methods are very often the computational demands in combination with a large volume of data or the requirement for time-series homogeneity. The problem of geotechnical measurements (field measurements in general) is the incompleteness of the time series or its inhomogeneity—different time spacing between recorded measurements. In this work, the PARCS (Paired Adaptive Regressors for Cumulative Sum) method was used, a method for detecting a change in the mean value of the multidimensional time series [20]. The time series of measurements are burdened with random errors. The mean value is then close to the value of the monitored variable in a large number of measurements.

The PARCS method looks for a predefined number of change-points in the time series. The input data to the algorithm are an array containing the time-series data and the scalar value indicating the number of searched change-points (N). The advantage over other methods is that the time series may not be completely consistent and may contain holes. The PARCS method was modified for data analysis so that the input parameter N was used only as an a priori (initial) value and the target mean error parameter was added. After the detection of change-points, the time series is divided into segments, which are approximated by linear regression. If the standard deviation on one of the sections exceeds the given limit value, the parameter N is increased (Figure 4).

The cut-off target value of the standard deviation was derived from the assumed accuracy of the dilatometer measurement as twice the standard deviation of the measurement found during laboratory testing [62].

A relatively large number of change-points was detected. Taking a closer look at the trajectory of motion, it is clear that some of these change-points do not lead to a change in the trajectory of motion but are the result of chaotical behavior. Therefore, we decided to directly analyze the trajectory of motion in the plane (2D) and space (3D) to find the key points of its change and their time stamps (Figure 5).

**Figure 4.** Modified PARCS method applied to the time series of residual values (HIE3 *z*-axis): (**a**) time stamps of detected change-points marked in green, target standard deviation of 0.05 mm; (**b**) comparison of detected time stamps and seismic events.

**Figure 5.** Example of the processing of the movement trajectory on a Z-Y plane on: (**a**) HIE2 device; (**b**) HIE3 device: raw trajectory (blue), smoothed trajectory using 30-day moving average (red), simplified trajectory (green).

The trajectory of the movement is represented by an array of coordinates with a removed periodic component. At first, the trajectory was smoothed using a 30-day moving average of the position to remove random noise. Then, the key breaking points of the motion curve were found, the curve was expressed by the minimum number of linear segments so that the accuracy criterion was preserved (Figure 5). For this purpose, the Ramer–Douglas–Peuckner (RDP) simplification algorithm was used [63]. The advantage of this approach is that all three components of motion can be analyzed together. Compared to other change-point detection methods applied to components, the RDP algorithm is simpler and enables very fast sampling of the curve and detection of change-points.

Importance of the change-points was rated to highlight the most significant changes in the pathway of movement of the monitored SAL detachment plane. We assume that the importance of change-point is increasing with its shortest distance from the connection of two adjacent change-points (Figure 6). The most significant changes in the trajectories were highlighted by red stars (Figures 10 and 11).

**Figure 6.** (**a**) Importance of change-points was quantified using minimal distance di; (**b**) change-point defined using error ellipse.

Note that in the figure above (Figure 6b), the change-point is displayed as one single moment, but a closer view to the trajectory of the movement shows, that in case of such slow movement (daily displacement is smaller than random noise) there should be some temporal buffer around the change-point where the trajectory is changing. This enabled to decide if the change of the trajectory is the reaction on a given event (for example seismic).

The time span of possible change-point location is in two- or three-dimensional space represented by error ellipse or ellipsoid, where its center is located in the change-point. The size of the major and minor axes of the ellipse/ellipsoid is derived from the magnitude of the standard deviations in the respective directions. All points of the trajectory located within the ellipse/ellipsoid are detected as possible change-points, and, therefore, each turn of the trajectory is represented by a set of relevant time stamps. As a consequence, one turn of the trajectory is now represented with various numbers of dates, which can be confronted with the time stamps of other events.

Both approaches: change-point detection in time series of components and trajectory simplification leads to detecting time stamps, which were then confronted with time stamps of above normal events in the time series of potential landslide triggers.

We assume that significant events in the movement are caused by significant events in the seismicity or precipitation. Statistics describe this relation as causality. We used Granger causality to quantify the relationship between sets of events with different extents, which allows determining if one set of values is caused by another on a given significance level [64].

There was a huge number of seismic events during the analyzed period—most of them with very low magnitude. To compare seismic events with significant events in the movement, it was needed to highlight above normal seismic events. For this reason, the average value and standard deviation of daily PGA maxima were computed. Then, only those seismic events were selected whose standard deviation of maximum PGA average was three times higher than the standard deviation. To highlight episodes with higher seismicity, five-day sums of PGA maxima were computed (Figure 7).

**Figure 7.** Cumulated PGA five-day sums and its moving 30-day average (red line).

#### **5. Results and Discussion**

Figure 8 shows the total displacement distance from the beginning of the monitoring. There is a visible decrease in movement speed on both devices. Although, in the case of HIE3 this is less pronounced.

**Figure 8.** Cumulated spatial displacement since the beginning of the monitoring: (**a**) HIE2 dilatometer; (**b**) HIE3 dilatometer; (**c**) spatial interpretation: red line—horizontal plane, blue line—vertical plane.

Both 3D dilatometers HIE2 and HIE3 showed permanent very slow creep movement during the analyzed period, generally in the downslope sense. It can be noted from Figure 8c, that on HIE3 device the horizontal component of the movement is parallel to the general slope. On HIE2 device the horizontal component of the movement is almost perpendicular to the general slope (sinistral strike slip). This is probably a long-term behavior, which was already documented by [47] and it is recognizable on a step-like shape of the gully below the device. On both dilatometric devices the movement direction corresponds to the El Hierro earthquake and volcanic activity epicenters, which are concentrated to the west and southwest from the SAL landslide. The speed of the movement and its direction changed over time. Movement trend on the HIE2 device is well described by a second order polynomial (Figure 9a–c). The X and Y axes show a close-to-linear decrease in the speed of the movement.

**Figure 9.** Decomposition of movements on the different axes showing the general decrease of the speed of the movement, black line—raw data; blue line—detrended data; red line—trend; green line—velocity.: (**a**) HIE2 *X*-axis; (**b**) HIE2 *Y*-axis; (**c**) HIE2 *Z*-axis; (**d**) HIE3 *X*-axis; (**e**) HIE3 *Y*-axis; (**f**) HIE3 *Z*-axis.

Analysis of GNSS and InSAR data from the beginning of the 3D dilatometer monitoring period [29] shows large absolute movements of the island in the period from July 2012 to March 2014. Although there are no permanent GNSS receivers on the SAL body, the nearest permanent GNSS stations show northeast shift and uplift of several centimeters. Similar results are shown from InSAR data, which are

available for the island. The results, however, show largest surface movements in the western and southwestern parts of the island, while the SAL body on the northeast was much more stable at that time. It has to be noted, however, that HIE2 and HIE3 automatic devices are measuring relative movements on the detachment plane, while GNSS and InSAR are showing absolute movements of the island.

PS InSAR analysis performed for the period from October 2013 to October 2015 [65] shows rather different behavior on different parts of the SAL body, suggesting that the landslide moves as separate creeping blocks (as also shown in [32]). The speed of the movement during the two-year analyzed period reaches few mm.a−<sup>1</sup> downslope.

It has to be noted that a principal road connecting the La Estaca port with the island's capital Valverde is crossing the SAL detachment plane just below the HIE3 dilatometer. This road has been damaged by cracks during 2012–2014 with general direction along the detachment plane, and the tarmac had to be repaired.

Movement trends on the HIE3 device are much more fluctuating (Figure 9d–f). This is not surprising considering the fact, that this device is installed in the open air and thus more prone to climatic influence (i.e., temperature changes, humidity). In the case of the *Z*-axis, the movement speed seems to slightly increase. This is most probably due to the fact, that this device is measuring a separated block within the entire landslide body, which exhibits slightly different behavior than the rest of the landslide [32].

It can be seen, that most significant changes in the direction of the movement happened in December 2013, at the very beginning of the monitoring. Comparison of the change-points in the trajectory on both HIE2 and HIE3 with the PGA and API are shown on Figures 10 and 11. There is a recognizable correlation of seismic activity and changes of movements on both devices. However, HIE2 shows more movement changes correlating with PGA than the HIE3 device. The most notable change-point probably influenced by the seismic activity occurred at the end of December 2013, when the strongest earthquake during the monitoring period (M = 5.1) also occurred. Other change-points presumably influenced by seismic activity happened in 14-Nov-2013; 23-Dec-2013; 24-Dec-2013; 25-Dec-2013; 28-Dec-2013; 16-Mar-2014; 04-Sep-2014; 30-Mar-2015; 08-Nov-2015; 13-Nov-2015; 21-Apr-2016; 08-Jun-2017; 14-Mar-2018; 14-Jun-2018; 15-Jun-2018 and 17-Oct-2018.

In case of API there is not a clearly visible correlation with change-points, except in December 2013; July 2016 and April 2018 affecting the HIE3 device and July 2016 affecting the HIE2 device. This suggests that the nature of the measured movement has, more probably, an endogenous cause. These observed results are in good accordance with the analysis performed by [32]. In their study, a period from October 2013 to June 2014 was analyzed on the HIE3 device. Three major impulses were recognized (rainfall event in December 2013 and seismic events in December 2013 and March 2014). All these three impulses were also identified using this proposed methodology.

**Figure 10.** HIE2—Time stamps of change-points (green lines) detected using the Ramer–Douglas–Peuckner (RDP) algorithm vs. PGA (**a**) and vs. antecedent precipitation index (API) (**b**). Red stars show the time stamps of five most significant changes of the trajectory.

**Figure 11.** HIE3—Time stamps of change-points detected using the RDP algorithm vs. PGA (**a**) and vs. API (**b**). Red stars show the time stamps of five most significant changes of trajectory.

However, during most of the analyzed period, there is no clear correlation between change-points and seismic or rainfall activity. This can be caused by two reasons. Firstly, the seismic activity (PGA), which is presumed to be the main triggering factor of SAL movement [49], has ceased significantly since 2013, as is decreasing the average speed of movement of the SAL documented in this work; changes in the residuals are negligible. Secondly, the influence of precipitation (API) is probably negligible on the SAL behavior. Only three events were identified, mostly affecting the HIE3 device. This is probably caused by the presence of a shallower landslide block influenced by precipitation (see [32] for more details).

#### **6. Conclusions**

The proposed methodology allowed us to process a large amount of monitoring non-stationary data and to find critical time stamps—change-points, when the nature of the monitored landslide movement changed significantly. Moreover, it suppresses the influence of the scientist, who is always subjective and to some point tends to find what he is looking for.

The PARCS method detects change-points very well but is relatively hardware-demanding. It has to be noted, however, that the magnitude of the detected changes is close to the confidence limit of the measurements determined in laboratory conditions [62]. It is worth noting the significant change in the trajectory of motion detected on both measuring devices in December 2013 was influenced by both precipitation (at the beginning of December) and strong seismicity (24–26 December).

Despite not acquiring unambiguous results in terms of the geological interpretation, we believe that the proposed methodological approach can significantly contribute to the landslide monitoring data analysis of any slow-moving landslides in the future by both:


Future research should focus more on defining the influencing magnitude of the triggering factors to determine better the behavior of the SAL in the future.

**Author Contributions:** Conceptualization, J.B. (Jan Blah ˚ut) and J.B. (Jan Balek); methodology, J.B. (Jan Blah ˚ut) and J.B. (Jan Balek); formal analysis, J.B. (Jan Balek) and M.E.; investigation, J.B. (Jan Blah ˚ut), S.M.; data curation, J.B. (Jan Blah ˚ut), J.B. (Jan Balek) and M.A., writing—original draft preparation, J.B. (Jan Blah ˚ut), J.B. (Jan Balek); visualization, J.B. (Jan Blah ˚ut), J.B. (Jan Balek), writing—revision, J.B. (Jan Blah ˚ut), J.B. (Jan Balek). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Czech Science Foundation, grant number GJ16-12227Y and by the long-term conceptual development research organization RVO: 67985891.

**Acknowledgments:** The authors are indebted to María José Blanco and her staff at the Centro Geofísico de Canarias (IGN) for their help with installation and maintenance of the 3D dilatometer monitoring points on El Hierro and IGN for the open seismic data availability. OGIMET is acknowledged for providing daily meteorological data from the El Hierro Airport. Jana Šreinová and Monika Hladká are thanked for processing and analyzing the 3D dilatometer data. We also thank to the two anonymous reviewers who helped to improve the clarity of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Integrated Field Surveying and Land Surface Quantitative Analysis to Assess Landslide Proneness in the Conero Promontory Rocky Coast (Italy)**

### **Francesco Troiani 1, Salvatore Martino 1, Gian Marco Marmoni 1, Marco Menichetti 2, Davide Torre 1, Giulia Iacobucci <sup>1</sup> and Daniela Piacentini 2,\***


Received: 16 June 2020; Accepted: 8 July 2020; Published: 13 July 2020

#### **Featured Application: Assessment of landslide proneness, with a particular emphasis on catastrophic landslide events.**

**Abstract:** Rock slopes involved in extensive landslide processes are often characterized by complex morphodynamics acting at different scales of space and time, responsible for different evolutionary scenarios. Mass Rock Creep (MRC) is a critical process for long-term geomorphological evolution of slopes and can likewise characterize actively retreating coastal cliffs where, in addition, landslides of different typologies and size superimpose in space and time to marine processes. The rocky coast at the Conero promontory (central Adriatic Sea, Italy) offers a rare opportunity for better understanding the predisposing role of the morphostructural setting on coastal slope instability on a long-time scale. In fact, the area presents several landslides of different typologies and size and state of activity, together with a wide set of landforms and structural features effective for better comprehending the evolution mechanisms of slope instability processes. Different investigation methods were implemented; in particular, traditional geomorphological and structural field surveys were combined with land surface quantitative analysis based on a Digital Elevation Model (DEM) with ground-resolution of 2 m. The results obtained demonstrate that MRC involves the entire coastal slope, which can be zoned in two distinct sectors as a function of a different morphostructural setting responsible for highly differentiated landslide processes. Therefore, at the long-time scale, two different morphodynamic styles can be depicted along the coastal slopes that correspond to specific evolutionary scenarios. The first scenario is characterized by MRC-driven, time-dependent slope processes involving the entire slope, whereas the second one includes force-driven slope processes acting at smaller space–time scales. The Conero promontory case study highlights that the relationships between slope shape and structural setting of the deforming areas are crucial for reaching critical volumes to induce generalized slope collapse as the final stage of the MRC process. The results from this study stress the importance of understanding the role of morphostructures as predisposing conditions for generalized slope failures along rocky coasts involved in MRC. The findings discussed here suggest the importance of the assessment of the slope instability at the long time scale for a better comprehension of the present-day slope dynamics and its major implications for landslide monitoring strategies and the hazard mitigation strategies.

**Keywords:** coastal landslides; mass rock creep; coastal cliffs; land surface analysis; data analysis; Conero promontory

#### **1. Introduction**

Landslides of different typologies and size can occur along active retreating coastal cliffs, where the gravity-induced processes often superimpose, in space and time, to marine processes in shaping the coastal slopes. Rocky coasts represent about 80% of coasts worldwide [1], and they represent a great part of coasts around the Mediterranean basin [2]. Rocky coasts have been broadly studied so far, representing zones with high economic, social, cultural, and touristic values. In particular, many studies focused on the factors controlling erosional processes along both hard and soft rock coasts [3], whereas other studies focused on the present morphodynamics, with a particular emphasis on the role of active landslides on the present evolution of the coastline and the associated hazards [4,5]. Nonetheless, the understanding of the role of the slope shape and tectonic structures on the long-term morphoevolution of coastal slopes is crucial to enhance the awareness of slope instability on rocky coasts.

Complex morphodynamics, implicating multiscale morphoevolutionary scenarios, often characterize rock slopes affected by extensive landslide processes. Mass Rock Creep (MRC) *sensu* [6], has been demonstrated to represent a critical process for the long-term morphoevolution of rock slopes in different geological contexts [7–9]. Large disruptive landslides may represent the ultimate stage of the MRC process that can culminate in paroxysmal events (like rock avalanches); nevertheless, these ultimate events are anticipated by gravity-driven deformations which represent the strain effect of the rock mass viscous rheology. When these deformations develop into abrupt and generalized slope failures, it is caused by an increasing strain rate over time [10]. In particular, threshold conditions can be reached when the stationary creep stage evolves into an accelerating creep stage [11,12]. Generalized slope failures are typically associated to mass strength reduction due to rock mass damage, the latter occurring over time as an effect of viscous deformation [13,14]. Catastrophic failures that occur in the accelerating stage of MRC process can progress rapidly into rock avalanches and slides that are generally due to the fragmentation of the involved rock masses throughout their propagation and emplacement [15].

Coastal slope processes can be prone to evolve in paroxysmal events [16–18] or can display a continuous slow activity associated with MRC process [19–22]. Along rocky coasts, proneness to landslides needs to consider the relations between the cliff rheology and lithostructural setting, as well as external stresses [23,24], such as rainfall, earthquakes, or storm sea waves. Furthermore, considering the coastal slope instability at a longer time scale, changes in sea level, occurring over a wide range of space–time scales [25], should be taken into account.

The investigation of rocky coasts often represents a research challenge due to the accessibility of rock cliffs in all their portions. Remote sensing techniques from satellite imagery and aerial photos can be sometimes unhelpful due to cliff attitude and/or in the case of high-angle sloping coasts, as for example plunging cliffs. [26] proposes a methodological approach for the continuous acquisition of geological, geomorphological, hydrogeological, and biological datasets along rocky coasts based on a snorkeling approach. The latter is helpful for a detailed mapping of landforms and processes at the foot of the coastal slopes. [27] proposed a combined approach based on snorkeling-based and traditional fieldwork in the inland for a detailed mapping of the rocky coasts affected by gravity-induced slope processes.

The rocky coast at the Conero promontory (Central Adriatic Sea, Italy) presents several features that account for the predisposing role of morphostructural setting on coastal landslide processes involving extensive rock slopes [28,29]. Moreover, this area provides a wide set of geomorphic and tectonic elements that can be regarded as conditioning factors for the MRC process evolution [27].

The main goal of this study consists in unraveling the role of morphostructures as a controlling factor for the coastal morphodynamics and the implication of the latter on the gravity-induced slope processes that act at different spatial and time scales. Different long-term morphoevolutive scenarios are considered based on different geomorphological/structural settings and the possible triggering mechanisms, then their implications on the present-day slope instability are discussed.

#### **2. Dataset and Methods**

A combination of laboratory- and field-based techniques consisting of land surface quantitative analysis and geomorphological–structural surveys were performed. The topographic dataset, available for the study area in vector format at the scale of 1:10,000, was used as a basemap for the field surveys. These surveys were preliminarily supported by the interpretation of panchromatic orthophotos, available for the years 1988–1989 and 1994–1998 at www.pcn.minambiente.it as a web map service (WMS), together with panchromatic aerial photos available at the scale of 1:33,000 (year 1954) and 1:75,000 (year 1989). Land surface quantitative analysis was performed, starting from a Digital Elevation Model (DEM) with a ground resolution of 2 m and derived from the LiDAR (Light Detection and Ranging) dataset available for the Italian coastline at www.pcn.minambiente.it. A Geographic Information System (GIS) using the ESRI® ArcGIS 10.6 desktop software platform was used for storing geomorphological and geostructural data derived from the surveys, other than for managing the whole dataset and elaborating the quantitative analyses.

#### *2.1. Geomorphological and Structural Survey*

A series of geomorphological surveys, organized between the spring and summer of 2019, were carried out for corroborating previous surveys [27,30], which mainly focused on coastal processes and landforms, with particular attention on the detection of tidal notches and their morphoevolutive significance. The surveys performed here aimed at refining the recognition and mapping of coastal and gravity-induced processes and landforms, as well as at enlarging this detection along the entire coastal slope between the Portonovo and the Due Sorelle sites (Figure 1).

**Figure 1.** Geomorphological map of the Conero promontory showing the main geological and structural elements (map coordinates are in WGS84 cartographic reference system). Near-coast bathymetry (0–10 m below the sea level) is also indicated. 1: study area; 2: lake; 3: normal fault; 4: strike-slip fault; 5: Maiolica Fm.; 6: Marne a Fucoidi Fm.; 7: Scaglia Rossa Fm.; 8: gully erosion; 9: tidal notch; 10: plunging cliff; 11: landslide-toe compressional ridges; 12: rocky scarp; 13: landslide trench; 14: swamp; 15: alluvial fan; 16: stratified slope deposit; 17: sea stack; 18: beach; 19: shore platform; 20: slope debris and talus deposits; 21: debris flow; 22: debris avalanche; 23: rock avalanche; 24: rock fall; 25: rock slide; 26: diffuse gravity-driven slope deformations.

The detection and mapping of the processes and landforms were based on the approach proposed by the new Italian Geomorphological Mapping guidelines [31], which suggest mapping each landform based on its genesis, evolution, and present dynamics. The peculiar amount of landslide phenomena affecting the coastal slopes imposed special attention on the detection and mapping of the gravity-induced processes and landforms, and their complex evolutionary scenarios. Geomorphological fieldwork was preceded by the interpretations of the panchromatic orthophotos visualized in the GIS platform and by the interpretation of the historical aerial stereopair images, which were helpful for a preliminary delimitation of both large landslide bodies and debris deposits. Results of these preliminary tasks were also helpful for tracing the main structural lineaments. The optical satellite images available in Google Earth, together with the ones available in ArcGIS 10.6 (i.e., TerraColor dataset) were consulted for planning fieldwork and for supporting the mapping of the largest landslide features, other than for delimiting the main cliff escarpments. The multitemporal remote sensing analysis was also useful for defining the recent activity of the main landslides. After the geomorphological fieldwork, final refining of the geomorphological mapping was completed through the 3D visualization of the topography at high-resolution, thanks to a hillshade map derived from the 2 m-cell-sized DEM in the GIS platform.

Geostructural surveys, conducted together with the geomorphological ones, were carried out to detect and characterize the main tectonic structures that strongly control the rock mass behavior, as well as for refining the mapping of the lithological units composing the bedrock. The final aim was to unravel the connection between the local bedrock discontinuity frequency and distribution and the rock slope deformation and instability.

The fieldwork was conducted using digital survey techniques, based on mobile devices that allowed for collecting field data in digital form, and its direct storage through remote cloud-uploading and real-time updating.

#### *2.2. Land Surface Quantitative Analysis*

The objective classifications of terrain types, at the base of the quantitative geomorphic analysis, require the use of both local and statistical terrain parameters helpful for defining the topographic signature of different surface processes and landforms [32]. Nonetheless, given the strong variability of landforms in size and tridimensional shape, due to different morphogenetic factors, many terrain parameters have been proposed to provide a quantitative description of topography [32,33]. In this research, two main geomorphometric parameters were selected among numerous DEM-based terrain parameters, considering the specific research aim as well as the spatial scale of the analyses performed. The selected parameters were calculated starting from the available 2 m-cell-sized DEM.

The first parameter is the slope angle, derived by means of the surface slope algorithm running in the Spatial Analyst extension of the GIS platform. For each DEM cell, this algorithm computes the first derivative of the elevation surface, or the maximum change in elevation between a cell and its eight neighbors (i.e., the steepest downhill descent from the cell) [34].

The second parameter is the surface roughness that can be successfully used to define different landslide features and to examine the landslide activity [35,36]. Different methods have been proposed to compute the surface roughness [37–39], although one of the most used methods is simply based on the statistical dispersion of heights or slopes. In this work, the surface roughness was calculated as the standard deviation of the elevation values within a 4 × 4 moving window. This method, among many calculation methods revised by [36], was demonstrated as one of those useful for evidencing the topographic signature of large landslides. In order to obtain the surface roughness values, the focal statistics tool running in the ArcGIS 10.6 platform was adopted.

For supporting the spatial interpretation of both parameters, color-coded maps were produced and superimposed with a hillshade map, thus favoring the visualization of both terrain morphology and the spatial variation of parameters. Descriptive statistics for both parameters were calculated and reported as an inset for each map. Finally, for both parameters the number of pixels per bin

(i.e., pixels in each parameter class) was reported using separate histograms for the southern and northern coastal sectors.

A composite altimetric profile parallel to the coastline was created, starting from the 2 m-cell-sized DEM. This composite profile reports the elevation distribution for progressive distance along three different tracks, thus allowing for a more accurate description of coastal morphology. The first track permits delineation of the cliff escarpment morphology (i.e., upslope zone); the second track favors the description of the morphology along the intermediate slope portion (i.e., midslope zone); and the third track portrays the base of the coastal slope, at about 50 m in from the present shoreline (i.e., footslope zone). The extraction and visualization of both slope angle and surface roughness values along two distinct transects support the study of the distribution of these parameters along different portions of the slope. This visualization increases the evidence of parameter value variation due to the presence of landslides of different type, size, and state of activity. The parameters were extracted along the midslope zone (i.e., at the +200 m a.s.l. contour line) and along the footslope (i.e., at the +20 m a.s.l. contour line).

#### **3. Study Area**

The Conero coastal sector stretches for about 5 km along the Adriatic Sea in the Central Italy, 10 km SE of the Ancona city, and represents its inflection point. Mt. Conero (43◦33 04 N, 13◦36 18 E) constitutes the highest relief along the Adriatic coast with an altitude of 572 m a.s.l. (Figure 1). The climate is of the Mediterranean type with mean annual precipitation of 780–790 mm and mean annual temperature of 14–15 ◦C. Summers are generally hot and dry, though a marked variability characterizes winters, mostly depending on the Atlantic cyclogenesis [40]. Prevailing winds approach from the SSE–SSW sectors [30] and account for the dominant SE main wave direction. During the winter months, sea-storm waves up to 3–3.5 m and with NE secondary wave direction can strike the coastline, during Bora wind conditions. The average tidal range in the Conero area is about 0.5 m and the action of tidal current is extremely low [27].

Concerning the geological setting, the Conero promontory consists of an asymmetric, NE–E verging anticline with NNW–SSE axial directions [29,41] and represents the easternmost outcropping part of the external stack of the Umbria–Marche foreland fold and thrust belt.

Limestones and marly calcareous terrains compose the bedrock of the study sector of the Conero promontory and are ascribable to the Umbria–Marche stratigraphic succession, from Cretaceous to Oligocene [29]. The complete stratigraphic sequence consists of the following formations, from oldest to youngest: Maiolica Formation (Fm.) (MA; Lower Cretaceous; prevalently composed of decimetric well-stratified limestones with chert and black shales), Marne a Fucoidi Fm. (MF; Lower Cretaceous; prevalently composed by centimetric marls and marly calcareous layers with black shales), and Scaglia Rossa Fm. (SR; Upper Cretaceous; composed of alternating and well-stratified centimetric to metric limestones, calcarenites, biocalcarenites and marly calcareous layers). In the area, this stratigraphic sequence presents some peculiarities with respect to the outcropping in the internal area of the Apennine chain. In particular, the thickness of the MF is about 15 m, whereas the regional thickness is generally up to 80 m. Furthermore, the paraconformity erosive contact between the MF and SR formations is likely caused by syn-sedimentary slumps [42]; in fact, in the area, the Scaglia Bianca Fm. together with a part of the basal portion of the Scaglia Rossa Fm. is totally absent [27]. The western side of the Mt. Conero anticline is gently dipping (around 25◦) and affected in the inner part by N–S right-lateral transtensive strike-slip faults with an offset of hundreds of meters. The eastern, coastward flank of the Conero asymmetric anticline is strongly inclined and is intersected by several transversal left-transtensive faults with an average E–W direction and offset a few km (Figure 1). The shear zones are characterized by cataclastic deformation, often with calcite veins and a thickness up to a few meters. The joint pattern is well developed in the rock masses with three prevalent directions NE–SW, NW–SE, and E–W. A pervasive cleavage with subvertical SW dipping planes is observable especially in the Maiolica and in the calcareous Scaglia Rossa bedding. Subhorizontal diagenetic

stylolites are quite common in the Maiolica Fm. Subvertical joints are developed mainly in the Scaglia Rossa Fm., with systems striking NW–SE and ENE–WSW, associated with breccias in related shear zones and systematic calcite veins.

The morphogenesis along the coastal slope is mainly conditioned by gravity-induced and marine processes, with the corresponding landforms and deposits (Figures 1 and 2). Landslides of different typologies, size, and state of activity are widespread in the area and principally characterize the entire coastal slope morphodynamics (Figure 1). Among the main mass movements, the landslides occurring in the northern coastal sector, at Portonovo and La Vela, are the largest ones (Figures 1 and 3A) [29]. Along the coastline, landslide deposits and plunging cliffs interrupt the poorly developed beaches and pocket beaches. The latter mainly consist of rocky blocks and gravels. Flat shore platforms, slightly inclined seaward, are also observable. Typical sea stacks at La Vela and Due Sorelle are located at short distances from the coastline (Figure 2). The current tidal notch occurs around these sea stacks [27]. The absence of either hanging or uplifted tidal notches, together with the absence of Tyrrhenian marine deposits and terrace remnants, testify of an active and intense morphodynamic in the entire coastal slope [27], mainly due to gravity-induced processes and the interrelated factors, such as geological factors or climate variations [43]. Several anthropic structures and man-made landforms, such as a local harbor, some touristic infrastructure, and parking areas, appear at the Portonovo landslide toe, close to the coastline.

**Figure 2.** Main landforms due to coastal processes at the foot of the cliff. (**A**) Shore platform along the central sector of the study area. Calcareous boulders due to recent rock falls are visible in the emerged portion of the platform. (**B**) Due Sorelle sea stacks and the adjacent pocket beach along the southern coastal sector.

**Figure 3.** Landforms and deposits due to gravity along the northern portion of the Conero rocky coast. (**A**) Panoramic view, towards the south, of the Portonovo and La Vela rocky coast (northern coastal sector). **(B**) Toe area of the Vela rockslide. In this sector, a late-glacial stratified slope deposit (i.e., èboulis ordonneè) overlays the landslide deposit. (**C**) Rock avalanche deposit at the toe of the Portonovo landslide exposed along a wave-cut scarp. The inset shows chaotic debris prevalently belonging to the Scaglia Rossa Fm. (SR) and incorporating a marly calcareous boulder belonging to the Marne a Fucoidi Fm. (MF). (**D**) Ridge-top depression in the Portonovo landslide backscarp area.

#### **4. Results**

The Conero promontory provided a set of data useful for better defining the morphoevolutive implications of the gravity-induced processes and landforms along the coastal zone. The entire coastal slope is affected by landslides of different typologies, which account for erosional landforms at the upslope zone and abundant landslide deposits along the mid- and footslope sites (Figure 1, Figure 2A, Figure 3, and Figure 4). Landslides are predisposed by low persistency of rock mass joints and by strata attitude whose dip is very close to the slope face direction.

The largest mass movements surveyed in the area occur in the northern coastal sector, between the Portonovo village and the La Vela sea stack (Figure 1), where two adjacent, large landslides masses can be observed and ascribable to rock-avalanche/slide mechanism [44], as they are associated with wide detachment scars that typify the sharp edge of the coastal cliff and delimit seaward the calcareous ridge top (Figures 1 and 3A). The northern landslide (Portonovo landslide) shows a runout of about 500 m that flowed seaward, accounting for the shaping of a bay where Portonovo village developed. A series of typical compressional ridges characterize the distal zone of this mass movement. The southern landslide (La Vela landslide), on the contrary, shows a truncated toe, though observing the trend and geometry of the contour lines down to −10 m b.s.l., the nearshore morphology of the seafloor indicates that a terminal landslide lobe is still detectable (Figure 1).

**Figure 4.** Landforms and deposits due to gravity along the southern portion of the Conero rocky coast. (**A**) Panoramic view of the Due Sorelle cliff. Falls and slides of rock blocks occur in the uppermost cliff portions, characterized by a sharp subvertical escarpment. Bedding attitude coincides with topography along the intermediate slope portion. At the cliff base, thick talus deposits occur, often interested by small debris flows and falls. (**B**) Evidence of a rock-fall detachment zone along the main cliff escarpment.

Patches of a thin layer of stratified slope-waste deposits form outcrops along the cliff base (Figures 1 and 2B) and unconformably overlay the chaotic landslide deposits. Stratified slope deposits are widespread within the Umbria–Marche Apennines, even along the coastline [45], principally where bedrock particularly prone to frost shattering, as Maiolica and Scaglia Rossa formations, is exposed. These deposits typify cold, periglacial environments [46] and represent in the Apennine area the termination of the last glacial slope–fluvial sedimentation cycle [47], thus their chronological collocation is at the end of Late Pleistocene [48].

A characteristic wave-cut scarp parallel to the present shoreline delimits the toe of both the landslide masses and provides suitable exposures of the deposits that here are those typical of rock avalanche/slide phenomena: disrupted rock masses, including rock blocks; chaotic calcareous debris with sandy clay; and a subordinately clay matrix (Figure 3B,C). Rock blocks and calcareous debris belong to the Scaglia Rossa Fm. and subordinately to the Marne a Fucoidi Fm., confirming that both bedrock formations are involved in the mass movement. At the toe of the Portonovo landslide, a deposit belonging to a shallow reactivation of the main landslide body uncomfortably overlays the rock avalanche/slide deposit (Figure 3C). Along the backscarp of both landslides, a series of trenches and ridge-top depressions mark the releasing zone (Figure 3D).

A plunging cliff characterizes the southern coastal sector, whose upslope and midslope zones are generally affected by rock-block slides and falls that strongly contribute to the cliff regression (Figure 4). A wide active talus deposit occurs at the footslope, often affected by secondary debris flows and fall phenomena. The main instability of the coastal cliff is strictly controlled by lithostructural factors, and just a few and limited zones along the cliff are subjected to the prevailing action of undercutting due to the sea actions.

Location and scale of gravity-induced instability processes appear related to the structural settings of the Mt. Conero anticline; in the southern sector, NE-dipping, high-angle bedding conditions the rock mass stability, providing the kinematic arrangement for planar and wedge-sliding, which seasonally affect the coastal slopes by mobilizing block volumes from tens up to hundreds of cubic meters. The repeated rock-block instabilities feed the slope system with a large amount of detritic material, which is mobilized by channelized flows after intense rainfalls.

The northern sector is instead featured by an E–ENE dip direction of bedding, resulting from fold-axis plunging and periclinal folding, as identifiable in the east escarpment of the Portonovo landslide. In this area, high-angle preorogenic transtensional faults [41] dissected the anticline forelimb, providing the inherited passive framework where slope-scale deformation evolved [49,50].

Evidence of slope-scale, gravity-driven processes occur along the entire coast where mesoscale features of MRC processes were surveyed (Figure 5).

**Figure 5.** Evidence of ongoing gravity-induced deformations along the coastal slope due to the MRC process. (**A**) Buckling folds affecting the calcareous layers at the Vela cliff. **(B**) Flexural toppling south of the Vela cliff. (**C**) Rock fall due to a buckling fold breakout. (**D**) Flexural block-toppling at the Due Sorelle cliff.

The structural observations indicate macroscopically ductile features of slope-scale flexural toppling along a steeply dipping and narrow-spaced primary discontinuity set, despite the presence of horizontal cross joints. The latter act as a breakout surface, driving the relative rotation of blocks with respect to a defined hinge surface [51]. These observations lead to a challenging distinction between flexural-toppling and brittle block-toppling [52], thus suggesting a transitional behavior strongly dependent on local lithological and joint density properties of the rock masses. Buckling folds also affect the calcareous bedding at the Vela cliff, causing tension and shearing fractures that can enhance the occurrence of rockfalls.

Along the analyzed coastal slope, the slope angle parameter shows a wide range of values (i.e., 0−88◦) and its spatial distribution is good evidence for the subvertical geometry of the upslope zone of the cliff (Figure 6A). In particular, the slope map reveals a very good match between the maximum values and the main rock-block detachment zones and the large landslide scar areas, respectively, in the southern and northern coastal sector. In the latter sector, the lowest slope angle values also appear at the foot of the cliff, delimiting the Portonovo landslide toe. Statistics show the prevalence of the highest slope values in the southern sector (Figure 6B), while within the northern sector, the values are distributed heterogeneously (Figure 6C) with a prevalence of the lowest values at the foot of the cliff (Figure 6A,C).

**Figure 6.** Spatial distribution of the slope values in the study area, based on a LiDAR-derived elevation dataset with ground resolution of 2 m. (**A**) Color-coded slope map overlain with a hillshade map, including a table with descriptive statistics. (**B**) Histogram with the distribution of slope values within the southern and northern sectors of the study area (y-axis reports the number of pixels per bin). (**C**) Histogram with the distribution of slope values within the footslope and upslope portions of the northern sectors of the study area (y-axis reports the number of pixels per bin).

The surface roughness parameter shows values up to 33 m (Figure 7). The highest values occur along the main rocky scarp, whereas the lowest ones occur in more homogeneous, flat or smoothed zones. The zones with the highest roughness coincide well with the ones with the highest slope values, along the main detachment areas for the rock-block falls. As for the slope parameter, even the surface roughness permits delineation of two almost homogenous zones along the coastal slope. The southern sector shows the largest values, while average small values are recorded in the northern sector (Figure 7B). More in detail, within the northern sector two distinct zones can be distinguished, based on surface roughness (Figure 7C). The first zone, involving the entire body of the La Vela landslide and the upper zone of the Portonovo landslide, rests within the base of the scar area downhill to the upper termination of the toe. The second zone coincides with the Portonovo landslide toe. This particular configuration is easily recognized along the transect at +200 m a.s.l. (Figure 8). This latter zone, along its northern portion shows similar values that range between 3 and 6 m; more in

detail, the values are the same where the transect crosses both the La Vela and Portonovo landslide bodies (Figures 1, 3A and 9). Along the transect at +20 m a.s.l., on the contrary, the roughness values are dissimilar for the toe areas of the two adjacent landslides, with the La Vela foot zone showing the same values as the uppermost transect and the Portonovo footslope zone showing lesser values. The highest roughness values occur along the transect at +20 m a.s.l. along the southern coastal sector, marking well the topographic setting of the base of the active retreating cliff.

**Figure 7.** Spatial distribution of the surface roughness values in the study area. (**A**) Color-coded surface roughness map overlain with a hillshade map, including a table with descriptive statistics. (**B**) Histogram with the distribution of surface roughness values within the southern and northern sectors of the study area (y-axis reports the number of pixels per bin). (**C**) Histogram with the distribution of surface roughness values within the footslope and upslope portions of the northern sectors of the study area (y-axis reports the number of pixels per bin).

The quantitative analyses performed by means of altimetric profiles in the coastline direction improved the interpretation of the outcomes described above, as well as allowed the definition of the different topographic and geometric characteristics of the analyzed southern and northern sectors, and the ones of the upslope and midslope zones (at +200 m a.s.l.), and the footslope (at 50 m from the present shoreline).

**Figure 8.** Composite altimetric profile of the entire coastal slope. The three superimposed profiles illustrate the morphology of the upslope, midslope, and footslope zones (the latter at 50 m in from the present shoreline). For the principal geomechanics measurement positions, the stereo plots illustrate the bedding attitude with respect to the slope geometry (topographic steepest descent). Below the profiles, the variations of both slope angle and surface roughness along the two distinct transects at the elevation of +20 and +200 m a.s.l. are reported.

**Figure 9.** Conceptual landscape model. The 3D sketch illustrates the different morphodynamic styles characterizing the southern (left sketch) and northern (right sketch) coastal sectors. Morphostructural settings strongly influence the types and development of different gravity-induced processes, as well as diverse landform associations and their spatial organization and evolution at different scales. Slope-scale, MRC-driven processes dominate the long-term morphoevolution of the northern coastal sector, whereas the southern sector is generally characterized by cliff-scale, force-driven processes. See Figure 10 for further details.

**Figure 10.** Conceptual space–time evolutive model for the active retreating coastal slope along the Conero promontory. Starting from the same process affecting the entire coastal slope, different long-term morphoevolutive scenarios are expected, depending on two different morphodynamic styles. The first scenario refers to slope-scale, time-dependent processes, while the second one refers to cliff-scale, time-independent (i.e., force-driven) processes. Both scenarios coexist along the analyzed coastal slope. The first one is typical of the northern coastal sector, where it favors active coastal slope downward erosion; the second scenario prevails along the southern coastal sector, where it accounts for active cliff regression.

The distribution of the slope parameter along the transect at +200 m a.s.l. (Figure 8) displays similar values, generally not exceeding 40–45◦ with two major exceptions at the Due Sorelle cliff, in the southern sector, and at the right flank of the Portonovo landslide, in the northern sector. Except for the Portonovo landslide toe zone (Figure 3A), where the average slope values are below 20◦, along the transect at +20 m a.s.l. the average slope values are >40–45◦ (Figure 8), highlighting the presence of a high-angle or subvertical active retreating cliff (Figures 2A and 9). The +20 m transect crosses either the base of the rocky slope or the talus deposits at the footslope in those sectors where it occurs (Figure 4).

#### **5. Discussion**

Two distinct sectors can be recognized along the Conero slope, typified by convex, dip-slopes [53]. The first one includes the northern sector embracing the Portonovo and La Vela zones, whereas the second sector encompasses the whole southern portion of the coastal slope at the Due Sorelle sea stacks. The two slope portions, featured by the same geological formations, show many differences in both topographic (Figure 8) and structural settings (Figures 8 and 9). Consequently, these sectors experienced different morpho-evolutive styles due to gravity-induced processes acting at different space- and time scales (Figure 10). In the northern portion of the area, at the slope scale, massive rock slope failures and their related scar areas mainly characterize the landscape of the coastal area, whereas cliff-scale, rock fall, and rock slide typify the landforms along the southern sector.

It is reliable to assume that a first-time massive rock slope failure, representing the generalized slope collapse as the final stage of the MRC process, characterizes the entire slope encompassing both Portonovo and La Vela localities with an estimated volume of about 40 <sup>×</sup> 10<sup>6</sup> m3. The latter likely developed in at least two steps (i.e., multistep first-time event): the first one involving the whole coastal sector and the second one only the Portonovo coastal slope. Both events likely occurred during the last glacial time. The lack of remnants of the marine terrace representing the Middle Pleistocene highstand [27] and the presence of the periglacial deposits sealing the main landslide bodies both support this assumption. Moreover, the relative freshness of the landslide-related erosional and depositional landforms allows exclusion of the collocation of the mass movements during earlier glacial cycles. The spatial distribution of surface roughness values along the northern coastal sector seems to confirm that the two landslide bodies belong to the same first-time event, and that only the Portonovo area experienced a second massive rock slope failure event. In fact, the La Vela landslide body shows nearly the same roughness values as the upslope portion of the Portonovo landslide, indicating a similar postevent remodeling, hence the same timing. Both landslide bodies were partially reactivated during historical time (mainly debris avalanche/slide phenomena) and field evidence confirms this observation, especially at the Portonovo locality. At this location, a catastrophic large landslide event occurred during the Middle Ages, destroying the local harbor and an ancient Benedictine Abbey [29].

From an evolutive point of view, the Conero coastal slope can be considered a typical palimpsest landscape (*sensu* [54]), where landforms belonging to the same process that make up the landscape (the gravitational one in the present case) developed at different spatial scales, have different age and, at present, display dissimilar dynamics. In fact, the nested landslide system that characterize the Portonovo and La Vela coastal slope likely started to develop during the beginning of the last glacial time, under cold environmental conditions and relatively low sea level (only the surficial landslide reactivations, though catastrophic, developed under the present environmental conditions). This morphodynamic style accounts for the down-wearing evolution of the rocky coast, with the upslope portion eroding downward as the footslope receives the supply of upslope sediments and advances seaward (Figures 8 and 9). On the contrary, the southern sector is a typical active plunging cliff, whose morphodynamic style, though dominated by the gravity-induced processes, is in part actively related to the coastal processes due to continuous sea wave action, tidal range, and storm waves, in particular along the footslope portions. Rock-block slides and rockfalls dominate the upslope portion accounting for the upper escarpment retreat, whereas at the base of the coastal slope, in addition to rock-block falls and slides, debris flows and falls also account for the cliff retreat. Several earthquake-induced landslides were reported for this portion of the Conero promontory [55], often triggered by the low-magnitude, high-frequency seismic events that typically characterize the seismicity of the external portions of the Apennine chain.

It is assumed that the dominant role of the different morphoevolutive styles along the Conero coastal slope is strictly related to the different structural settings (Figures 8 and 9) that account for the rather different development of MRC processes.

A double scale morphodynamic conceptual scheme is proposed in Figure 10, where time is reported on the x-axis and space on the y-axis. According to the time scale, the processes considered evolve from gravity- to force-driven during recent times, as these processes involve more limited portions of slopes, causing smaller size landslide events. These last events can be triggered by time-independent local and transient forcing, such as sea storms and earthquakes, while the larger size rock mass failures (i.e., rock avalanches and slides), which can be observed in the northern sector of the Conero promontory, are substantially triggered by inertial forces as they result from MRC, causing viscous deformations over time (red arrows in the conceptual scheme, Figure 10). For a higher dip of strata attitude, the MRC deformations (occurring at the entire coastal slope scale) do not envelop critical volumes that can become prone to failure over a longer time [7]; in this last condition, the morphoevolution results in several small-size failures (occurring at the cliff scale), which are represented principally by the rock-falls and slides in the southern sector of the Conero promontory (green arrows in the conceptual scheme, Figure 10). Reactivations are also possible, after initial mass rock failures, due to debris avalanche/slide or debris flow occurrence. Concerning the

climatic context and referring to the last glacial–interglacial cycle, the low-stand phase corresponded to the gravity-driven deformations evolving toward failures, while the high-stand phase mainly corresponds to the cliff-scale processes. The implications of this conceptual model in terms of risk are related to a higher probability of occurrence for the cliff-scale processes, acting at an entire slope scale, which represent the most hazardous ones with respect to the gravity-driven MRC deformations and failure.

The coexistence of large-scale MRC deformational processes and small-scale, time-independent failures is clearly visible in the Due Sorelle sector (southern portion of the Conero promontory) where observation of flexural toppling and buckling demonstrate that progressive rock mass deformations rapidly evolve toward failure (i.e., either falls or slides), thus avoiding larger volumes involving creep processes that reach the critical size for a generalized slope collapse.

Traditional field-based surveying and mapping with quantitative geomorphic analyses based on a high-resolution digital elevation dataset were integrated to represent a more suitable tool for better defining the morphostructural setting of the area under investigation, even in those sites where the field surveys were difficult to perform. DEM-based quantitative land surface analysis, thanks to different visualization techniques and the computation of two geomorphometric parameters, revealed a high-performance mapping methodology that allows for rapid completion of the time-consuming traditional techniques used in refining the survey and mapping slope processes and landforms due to gravity.

Results of this research can be valuable for coastal landslide analysis and propose a conceptual scheme for assessing and monitoring the current dynamics of slopes affected by mass movements of different typologies and size. In particular, this research reveals the importance of the evaluating the long-term morphoevolution of the entire slope for better comprehending the present activity and future development of the coastal landslides.

#### **6. Conclusions**

Based on field evidence, MRC processes involve the entire coastal slope along the Conero promontory. The coastal slope can be zoned into two distinct sectors, depending on different morphostructural settings and consequent rock slope instability processes. For the long timescale, two separate morphodynamic styles can be delineated along the coastal slope, which are each associated with specific evolutionary scenarios. The first scenario represents MRC-driven, time-dependent slope processes acting at the slope scale, whereas the second one includes force-driven, time-independent slope processes acting at smaller space–time scales. First-time rock slope failures (i.e., rock avalanche/slide mechanism) characterize the northern area of the promontory, which is interpreted as the result of extensive, multistep catastrophic failure that occurred at the final stage of the MRC process. Shallower landslides (i.e., debris avalanche/slide) can locally reactivate the main landslide masses. First-time rock-block slides and falls, occasionally triggered by earthquakes and/or storm sea waves, occur mainly along the southern sector of the coastal slope where they account for the active cliff retreat.

The present study highlights that slope shape and geostructural conditions of deforming slopes are concurrent in allowing the critical volumes for failures to be reached, thus inducing generalized slope collapse that can be regarded as the ultimate stage of MRC processes.

The findings discussed here also confirm the relevance in understanding the role of geomorphological vs. structural features as predisposing factors for slope collapse along rocky coasts. Moreover, the assessment of long-term slope instabilities allows for a better comprehension of the present-day, gravity-induced slope dynamics, suitably supporting landslide monitoring strategies and appropriate mitigation strategies even based on stress-strain time-dependent numerical analyses.

**Author Contributions:** Conceptualization, F.T. and S.M.; methodology, F.T. and S.M.; validation, M.M.; laboratory analysis, D.T. and G.M.M.; field investigation, F.T., D.P., M.M. and D.T.; data curation, G.I.; writing—original draft preparation, F.T., S.M., G.M.M. and D.P.; visualization, D.P. and G.I. All authors have read and agreed to the published version of the manuscript.

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

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Physical Model Experiments on Water Infiltration and Failure Modes in Multi-Layered Slopes under Heavy Rainfall**

#### **Junfeng Tang 1, Uchimura Taro 1, Dong Huang 2, Jiren Xie <sup>3</sup> and Shangning Tao 1,\***


Received: 30 March 2020; Accepted: 11 May 2020; Published: 17 May 2020

#### **Featured Application: Authors are encouraged to provide a concise description of the specific application or a potential application of the work. This section is not mandatory.**

**Abstract:** To assess the influence of an intermediate coarse layer on the slope stability during heavy rainfall, knowledge about water movement and how slope failure occurs is important. To clarify the characteristics of water infiltration in a multi-layered slope and assess its influence on the slope failure modes, eight groups of physical slope models were investigated. It was found that the unsaturated hydraulic conductivity in the coarse layer (5.54 <sup>×</sup> 10−<sup>6</sup> cm/s) was much lower than that of the fine layer (1.08 <sup>×</sup> 10−<sup>4</sup> cm/s), which resulted in the capillary barrier working at a lower water content. Intermediate coarse layers embedded between finer ones may initially confine the infiltration within the overlying finer layers, delaying the infiltration and eventually inducing a lateral flow diversion in the inclined slope. Two different failure modes occurred in the model experiments: surface sliding occurred at the toe in the single-layer slope group and piping occurred at the toe in the multi-layered slope as the rainfall water accumulated, was diverted along the interface, and then broke through in the downslope direction of the intermediate coarse layer. The lateral flow diversion caused by the capillary barrier and the tilt angle may be the major factors influencing the difference of the failure modes. The result also revealed that the coarser layers may have negative effects on the slope stability.

**Keywords:** unsaturated soil; capillary barrier; multi-layer slope; slope failure

#### **1. Introduction**

Rainfall-induced slope failure is one of the most destructive natural disasters that occur in shallow natural slopes. The impacts of such catastrophic events are known, but these recurring natural hazards still result in many significant casualties and economic losses [1,2].

This study deals with a slope consisting of a fine layer and an intermediate coarse layer. The presence of soil layers with different unsaturated hydraulic conductivities in a shallow depth of soil affects the process of water infiltration, distribution, and pore water pressure in the slope, which results in different failure modes [3]. Rainfall infiltration water is due to the build-up of capillary barriers [4], which accumulate at the interface between fine and coarser soil layers. The capillary barrier effect has been widely studied in terms of its use in cover systems in waste disposal sites [5,6], whilst more recent research has focused on slope stabilization [7].

Capillary barriers can maintain a high degree of saturation in the soil above them, which results in different failure parts in a multi-layer slope [8]. These phenomena are related to the capillary tension, which limits the downward movement of the wetting front from a finer soil into underlying coarser soil. In an inclined interface, under continuous water infiltration from the slope surface, the accumulated water above the interface between the fine and the coarse material leads to the formation of a gravity-driven sub-surface water flow along the interface. This depends on the geometrical or boundary variations, at a certain distance downslope, often reported as the "diversion length", which can be estimated using a model proposed by Ross [9]. In natural slopes, characterized by a length in the order of hundreds of meters and an irregular layer geometry, the conditions leading to the penetration of infiltrating water into the underlying coarse layer are mainly governed by the materials, slope angle, and infiltration rate [10].

In some field investigations, a multi-layered slope with different hydraulic conductivities is a common situation in layered hillside slopes. Such a slope consists of fine sand deposits and medium-coarse sand deposits that are a few meters thick [11]. Additionally, the inclined angles range from 10 to 35 degrees. The coarse sand deposits involved have been found to have large pores and a higher hydraulic conductivity in a saturated condition, which have an important influence on the subsurface water flow and water content distribution during the infiltration, steady percolation, and drainage [12,13]. A schematic of a three-layer distribution in a slope is presented in Figure 1.

**Figure 1.** (**a**) Schematic profile of a multi-layered slope; (**b**) a description of the geometry of a multi-layered slope.

Many previous field investigations have been conducted and have found that layered soil influences the process of rainfall water movement and the distribution of water content, which determines the slope failure. However, clarifying the influence of an intermediate coarse layer on slope failure is still a complex issue, since many uncontrollable factors exist in a natural slope, such as the slope angle, rainfall intensity, etc. Hence, slope model experiments replacing field monitoring have been conducted, as these allow the same slope in different conditions to be investigated, reduce the cost by reducing the duration, and consider variable slope types [14].

To investigate the infiltration process in a layered soil profile in simplified and known geometrical and boundary conditions, physical model experiments have been performed by many researchers. Infiltration experiments have been conducted and have indicated that the wetting front temporarily stops above the interface of the fine layer and coarse layer, and the infiltration rate slows down [15,16]. Additionally, the pore water pressure head of the finer layer above the soil interface could not increase when the water started to infiltrate into the coarse layer [17]. The presence of a coarser lower-most layer may confine the infiltration within the upper finer layer up to a high degree of saturation. This capillary barrier effect occurred in a slope and was considered to be the cause of landslide initiation [7]. A capillary barrier at the upper interface of a coarse layer could have developed, favoring the accumulation and a lateral distribution of infiltrating rainfall and a possible diversion of flow down the slope, thus leading to a localized increase in the water content and loss of strength [18]. A diversion of the flow from the vertical direction towards the slope direction occurs and part of the water crosses the coarse layer infiltrate into the deepest layer. The settlement of soil is accompanied by progressive saturation of the soil, but general slope failure does not take place in the whole process [19].

Although some model experiments have been conducted to analyze the water movements and post-failure evolution due to the capillary barrier effect, few model experiments have focused on the slope failure time and modes in single- and multi-layer slope conditions or analyzed how an intermediate coarser layer may affect the slope stability.

Therefore, the aims of this study were as follows: (1) to experimentally clarify the influence of intermediate layers on water movement in a horizontal and inclined slope, respectively, especially after the breakthrough of the capillary barrier and (2) to investigate its influence on slope failure modes and the failure time in a single-layer and multi-layer slope by monitoring the water content, pore water pressure, and failure modes at different tilt angles. Soil-water characteristic curves (SWCC) [20] of fine sand and coarse sand were also measured to explain how the capillary barrier works. The proposed physical model experiments can provide a perspective on the failure time and modes of a single- or multi-layered slope during a rainfall event.

#### **2. Materials and Methods**

#### *2.1. Basic Material Properties*

Silica sands S1 and S7 were chosen for the physical model experiments as the fine material and coarse material, respectively, by the difference of D50 and size distribution. Figure 2 shows the size distribution of S1 (fine material) and S7(coarse material) used in model experiment. Table 1 summarizes the specific properties of S1 and S7 used in the model experiments. The sieve tests were conducted using the JGS Geotechnical Society standard test methods (JGS0131-2009).

**Figure 2.** Grain size distribution of S1 and S7 in physical model experiments.

**Table 1.** Basic properties of sand materials.


The unsaturated hydraulic properties were measured by the variable head method (ASTM 2006) D2434-68. Table 2 and Figure 3a show the fitting parameters of the VG (van Genuchten–Mualem) model [21] and soil-water characteristic curves of silica No 1 and 7 in both the drying process and wetting process. The unsaturated hydraulic conductivities of sands were obtained by a pressure plate apparatus in the lab. As shown in Figure 3a, the drying (desorption) process and the wetting (absorption) process of the SWCCs of sand cause hysteretic behavior [15] for the same suction value, and the sand can retain more water in the drying process than in the wetting process. In the wetting process, the air entry value (AEV) of the fine layer (silica No 7) is about 1.5 kPa. The AEV of sand and gravel is about 0.45 kPa, being lower than the AEV of the fine layer. The capillary barrier effect occurs at the fine-coarse sand interface during rainfall infiltration. In Figure 3b, the unsaturated hydraulic conductivities of two types of sand are reported. The obtained values range between 10−<sup>2</sup> and 10−<sup>7</sup> cm/s under different suction conditions. The results show that the hydraulic conductivity of the coarse sand is higher than that of the fine sand at almost saturation, while it is significantly lower when the soil is unsaturated.

The soil-water characteristic curves of the soil were modeled with the van Genuchten–Mualem model [5], as follows:

$$\mathcal{S}\_{\varepsilon} = \{1 + (-ah)^{\eta}\}^{-m}, \mathcal{S}\_{\varepsilon} = \frac{\partial - \partial\_{r}}{\partial\_{s} - \partial\_{r}} \mathcal{m} = 1 - \frac{1}{n}.\tag{1}$$

In the above equation, the water retention curve is expressed in terms of the effective degree of saturation. θ is the volumetric water content; θ*<sup>r</sup>* and θ*<sup>s</sup>* indicate the residual and saturated values of the water content, respectively; *a*, *m*, and *n* are the fitting parameters; *h* is the matric suction; *a* is a scaling parameter (units of m<sup>−</sup>1); and the exponents n and m are parameters that determine the shape of the retention curve. The hydraulic parameters are given in Table 2.


**Table 2.** Hydraulic properties of sand materials.

**Figure 3.** (**a**) Soil-water characteristic curves (SWCCs) of silica No 1 and No 7 in the drying and wetting process; (**b**) unsaturated hydraulic conductivity of silica No 1 and No 7 in the wetting process.

#### *2.2. Water Flow and Slope Failure in Multi-Layer Slope Models*

The multi-layer models were set up to evaluate the effect of the capillary barrier acting on a slope. This model system was built to show the advancement of the wetting front, and monitor the volumetric water content in the soil, pore water pressure, and slope failure process. In addition, the failure process was directly recorded by cameras to obtain a better understanding of failure modes and the effect of the capillary barrier. Table 3 summarizes the physical experimental conditions.


**Table 3.** General information for the slope model experiment.

#### 2.2.1. Flume Model System

Figure 4a shows the apparatus used for the physical model experiments, which consisted of an inclined steel box, a rainfall simulation system, and a set of pore water pressure and VWC (volumetric water content) sensors. The details of each subsystem are as follows: (i) the inclined steel box had dimensions of 1.0 m (length) × 0.3 m (width) × 0.5 m (height); (ii) the sidewalls of the box were made of an acrylic plate to observe the advancement of the wetting front and failure process during rainfall. The gap between the steel plates and the acrylic plate was sealed with epoxy adhesive.

**Figure 4.** (**a**) Schematic diagram of the experimental apparatus used for a multi-layered slope under rainfall: side-view of the multi-layered slope; (**b**) side view of the experimental apparatus used for the slope. The yellow dashed line shows the fine–coarse interface.

Furthermore, the inclined angle of the model box could be lifted by a crane to simulate the inclined cover from 0 to 60 degrees. Experiment pictures were recorded by cameras in different locations around the model. Table 3 summarizes the physical model experimental conditions.

#### 2.2.2. Rainfall Simulation System

A rainfall nozzle was placed 60 cm above the model box to simulate rainfall with a constant intensity. The intensity and duration of rainfall were controlled by a control value and air pressure gauge. The rainfall intensity was kept at a constant intensity by air pressure, and ranged from 35 to 100 mm/h. The sensors ECH2O EC-5 were used to determine the volumetric water content. Soil-specific calibration is recommended for obtaining the best possible accuracy in volumetric water content measurements [22,23]. Calibration of the EC-5 sensors has been shown to result in an increased accuracy of 1–2% for all soils with soil-specific calibration [24,25]. The calibration of EC-5 sensors employed in silica No 1 and No 7 is shown in Figure 5b.

**Figure 5.** (**a**) Nozzle body and spray pattern and distribution for rainfall simulation; (**b**) calibration of the volumetric water content sensor in silica No 1 and No 7.

#### 2.2.3. Theory of Measurement Devices

The instruments used in the model experiments were calibrated before installation, including the pore water pressure, tilt sensors (Figure 6a), VWC sensors (Figure 5b), and rainfall simulators (Figure 5a). The intensity and uniformity of artificial rainfall created by the simulator were calculated based on the weight of the sample at a certain time.

**Figure 6.** (**a**) MEMS (Micro-Electro-Mechanical System) tilt sensor and its working theory; (**b**) typical pore water pressure response in soil and the structure and size of the transducer.

#### 2.2.4. Testing Procedure

#### Soil Preparation

S1 and S7, which were used to make the slope, were dried in an oven for 48 h. Then, an amount of water was added to the soil to achieve an initial water content of around 6%.

#### Compaction of Soils

The prepared soil was compacted and placed in a series of horizontal layers. Silica No 7 and No 1 were placed in the model box in layers and compacted to achieve a dry density of 1.33 and 1.43 g/cm3, respectively. Each layer was tamped equally with a rod to a thickness of 5 cm and these procedures were repeated until the height of the slope was achieved.

#### Positions of Sensors and Cameras

During the soil placement, VWC sensors, tensiometers, and tilt sensors were placed at specific locations in the three different layers and the time of recording the quantity of water content was 10 s.

Table 3 summarizes all of the experimental conditions.

#### **3. Results**

#### *3.1. Failure Situations in All Cases*

Figure 7 shows the failure situations for different cases during rainfall. No runoff on the surface of the slope was observed in all cases, as the infiltration capacity of the soil layer was higher than the rainfall intensities used in this study. Two failure modes were observed in this study. One was soil piping, which occurred in the multi-layered slope. Slight soil piping occurred at the toe of the slope, finer materials mixed with amounts of infiltrate water flowed out from the piping as the seepage surface, and the seepage surface grew gradually and increased until cracks appeared (Figure 7c). As Figure 7a shows, the piping phenomenon occurred slightly, and as the tilt angle increased, the soil pipe size developed rapidly and more seriously (see Figure 7b,c).

**Figure 7.** Comparison of failure situations of the different model experiment cases (each case was performed three times to ensure the repeatability of the results, and typical results of individual experiments were chosen to show the failure situation in each case). (**a**) case I (tilt angle α = 0◦); (**b**) case II (tilt angle α = 7◦); (**c**) case III (tilt angle α = 15◦); (**d**) case V (tilt angle α = 0◦); (**e**) case VI (tilt angle α = 7◦); (**f**) case VII (tilt angle α = 15◦).

Another failure mode was surface sliding from the toe of the slope in a single-layer slope. Small slide failure occurred at the toe at first, and then relatively larger slide failure followed. This type of failure mode was clearly observed for cases IV, V, and VI, as shown in Figure 7d.

The failure modes and initial failure time *t* for all eight cases are summarized in Table 4. The initial slope failure occurred at the toe of the slope in case I and case IV at about 1.18 and 0.95 h after the rainfall has been applied, respectively. This means that failure in a multi-layer slope occurs later than that in a single layer under the same condition in the horizontal group. In Figure 7b,e, in case II and case V, failure occurred at around 0.84 and 0.89 h, respectively. In case III and case VI, more rapid movement of slope failure occurred in an almost fully saturated condition and they thus had a lower strength. The failure times of the 15 degrees inclined group were 0.78 and 0.85 h. It should be noted that failure occurred earlier in the multi-layer slope than in the single-layer slope in the inclined group, which is different from the result obtained for the horizontal group.


**Table 4.** Failure conditions of physical model experiments.

#### *3.2. Profile of the Volumetric Water Content in a Slope*

During the tests, the hydrological response was monitored by means of VWC (volumetric water content) sensors located at different locations within the slope and crossing the entire soil thickness, allowing the retrieval of volumetric water content profiles at different depths. The VWC sensors in section I in different cases were shown to explain how the capillary barrier effects influence the water infiltration and distribution during the rainfall and drying process, as is shown in Figure 8 (the sensor locations are shown in Figure 4b). Throughout the experiment, the VWC of the toe of slopes increased slowly with time toward a saturated value in response to the saturation process, until failures occurred. It should be noted that the capillary barrier clearly controlled the rate of changes of the VWC in a multi-layer slope, which determined the failure time of slope.

**Figure 8.** *Cont.*

**Figure 8.** Time histories of VWC during rainfall and the drying process in case I, III, V, and VII (the left figures show the time from 0 to 8000 s, while the right figures show the whole wetting and dry process in the rainfall event). (**a**) Case I. VWC trends in a multi-layer slope in the flat group (tilt angle α = 0◦); (**b**) Case V. VWC trends in a single-layer slope in the flat group (tilt angle α = 0◦); (**c**) Case III. VWC trends in a multi-layer slope in the inclined group (tilt angle α = 15◦); (**d**) Case VII. VWC trends in a multi-layer slope in the inclined group (tilt angle α = 15◦).

In case I and case V (see Figure 8a,b), the advancement of the wetting front is evident from the time history of the VWC sensors. It shows that the rise of the VWC was rapid, and then became gradual, as the toe of the slope approached a fully saturated condition. When the VWC at point L reached around 0.4, failure occurred at the toe of the slope. In addition, compared with the VWC at point L in the flat group, point L reached the saturated condition later in the multi-layer model, since the capillary barrier prevented the water from infiltrating into the bottom. This made the slope more stable and caused a delay in the failure time, which was 1.2 and 0.92 h in case I and case III, respectively. This capillary barrier effect can also be explained by the VWC histories obtained during the experiments. For example, it took around 26 min for the wetting front to pass through the interface from point B to point C [11].

Figure 8c,d show that when rainfall was applied, the VWC increased quickly above the interface (point B), and when rainfall was stopped, the soil above the coarse layer (point B) remained wetter than the same location in a single-layer slope. Additionally, for the VWC at the toe of the slope (point L) in the inclined group, the failure occurred earlier (0.78 h), while it was 0.84 h in a single-layer slope, which was contrary to the flat group experimental results.

#### 3.2.1. Beginning of Rainfall (t = 0 h)

In case I and V, the VWC in the bottom of the upper fine layer (point B) increased from 0.05 to 0.33 with the depth in the upper fine layers (height was from 25 to 45 cm), while that in the coarse layer was about 0.03 before capillary barrier breakthrough at the beginning of rainfall in case I (t = 0.5 h). This indicates that the wetting front above the interface was stopped for a while and rainfall water was stored in the upper fine layer due to the influence of the capillary barrier. In the lower section I (height was from 0 to 20 cm), the VWC maintained a constant condition while the capillary barrier was present.

#### 3.2.2. The Breakthrough of the Capillary Barrier (Case I, II, and III)

During the first 0.8 h of rainfall, as shown in case I and case II in Figure 9 a,b, the presence of the intermediate coarse layer caused a significant time delay in the infiltration process: the VWC in section I above the interface increased to 0.34 in case I, whereas the VWC in the coarse layer changed from 0.03 to 0.05 in case I at the depth of 22 cm, while the VWC was from 0.04 to 0.22 at the same depth in case II (t = 0.6 h). In this stage, the interface between the fine and coarse layer acted as a capillary barrier and gradually broke when the bottom of the upper fine layer was almost at saturation. In case II, the VWC in section I above the interface increased from 0.06 to 0.32, while that of the coarse layer was from 0.03 to 0.04. A lateral diversion flow may have occurred along the inclined interface, resulting in capillary barrier breakthrough occurring later in section I compared with case I.

**Figure 9.** *Cont.*

**Figure 9.** Profile of the VWC in section I in different cases (the VWC distribution after slope failure was verified by obtaining soil samples from various locations at the end of the experiments). (**a**) Case I; (**b**) Case II; (**c**) Case III; (**d**) Case V; (**e**) Case VI; (**f**) Case VII.

#### 3.2.3. Post-Breakthrough of the Capillary Barrier (t = 0.6–0.78 h in Case I, II, and III)

After 0.6 h of rainfall, the VWC of the upper fine layer maintained a constant value in section I in case I and III, since it reached steady state infiltration in the soil. Additionally, rainfall water infiltrates into the next layer following capillary barrier breakthrough. For example, the VWC of the coarse slayer increased from 0.05 to 0.12 at 0.83 h and 0.13 at 1.2 h (see Figure 9a) in case I and the VWC of the coarse layer increased from 0.05 to 0.07 at 0.78 h in case III. It should be noted that the VWC of the bottom of the lower finer layer increased to 0.36 at 0.8 h in case I and to 0.08 at 0.78 h in case III, respectively. This means that rainfall water could not arrive at the lower fine layer along section I because of the lateral water flow along the inclined interface in case III.

#### 3.2.4. Failure Occurred (t = 0.78–1.2 h)

In this part, the failure process is analyzed using the response of the tilt angle, pore water pressure, and VWC measurements in different locations in the slope. In case III and VII, the experimental flume was tilted to 15◦ and subjected to the same rainfall intensity. This section aims to evaluate the effects of an inclined angle on the water movement and slope stability.

The main results of this part are reported in Figures 10 and 11. The initial pore water pressure values were negative at the locations of S-I and S-II. The pore water pressure and VWC increased during the infiltration and then stabilized when the infiltration reached a steady state condition. During the rainfall, the tilt angle of point G and point A almost maintained a constant value, while the soil in these locations maintained a negative value. A sudden change occurred when the pore water pressure approached 0 kPa at point L in all cases.

Comparing the difference of the pore water pressure measured at S-II in case I and case V, presented in Figure 10a,b), the pressure grew from −4 kPa at the beginning of the rainfall to around −1 kPa after around 1.2 h above the coarse interface, while it reached −3.5 kPa in case V. This suggests that there was a higher water content and pressure head above the coarse layer than at the same location in a single-layer slope, as the capillary barrier effect led to the storage of water and high water pressure head. In this respect, it is worth noticing that in case I, the pore water pressure at point L (S-I) increased suddenly after capillary barrier breakthrough at the interface.

**Figure 10.** Volumetric water content, pore water pressure, and tilt angle change trends at different locations during infiltration. (**a**) case I (tilt angle α = 0◦, multi-layer), where the pore water pressure increased after capillary barrier breakthrough; (**b**) case V (tilt angle α = 0◦, single layer). For the positions of sensors, see Figure 3b.

The same behavior, in terms of the pore water pressure and VWC, was observed in tests with the tilt angle α = 15◦ (Figure 11). In these two cases, the greater slope angle induced the most sudden variations of water content and pore water pressure, and also resulted in an earlier response of tilt sensors located at the toe of the slope. For case III with the tilt angle α = 15◦, failure occurred in about 0.78 h, and the increase of VWC was more rapid than that of case VII, which was contrary to case VII. In the inclined group, the later failure occurred at 0.84 h in a single-layer slope, which proved that an inclined multi-layer slope is more dangerous than a single-layer one under this condition, which is contrary to the findings of the flat group.

It is evident in Figures 10a and 11a that the increase of the pore water pressure can be divided into two periods. A negative pore pressure was measured at the beginning of rainfall (S-I), and then, around 30 min later, a sudden increase of the pore water pressure value Δu could be measured after capillary barrier breakthrough at point L (Figure 3b). Moreover, a few minutes after breakthrough, piping occurred for a pore pressure of 0 kPa at the toe, which seems to indicate a negative influence for the slope instability. Similar to case I, pore pressure sensor S-I in case III also showed the same trends after the breakthrough of the capillary barrier. It is clear from both data sets that the pore water pressure value increased from −1.5 to −0.3 kPa in case I and from −1.1 to −0.21 kPa in case III, respectively.

**Figure 11.** Volumetric water content, pore water pressure, and tilt angle change trends at different locations during infiltration. (**a**) Case III (tilt angle α = 15◦, multi-layer), where the pore water pressure increased after capillary barrier breakthrough; (**b**) case VII (tilt angle α = 15◦, single layer). For the positions of sensors, see Figure 3b.

#### 3.2.5. Capillary Barrier Restoration in the Drying Process (t = 3.5–36 h)

The VWC in soil decreased once the rainfall stopped in case I (Figure 9a) and V ((Figure 9c). However, the bottom of the upper fine layer maintained a higher water content condition at 0.33 compared with the single-layer slope, which was only 0.19 in case I and V, respectively. The capillary barrier broken through could be restored to its pre-breakthrough condition once rainfall had stopped. The VWC of the finer soil above the interface decreased as the rainfall water continued to drain out. As the VWC at the interface decreased, the unsaturated hydraulic conductivity of the coarser lower layer also decreased, and eventually approached zero. The capillary barrier was completely restored when the intermediate coarser layer could not accept any more rainfall water from the overlying fine layer.

#### *3.3. Influence of the Tilt Angle*

#### 3.3.1. Influence of the Tilt Angle on Water Movement

To show the water movement and its lateral diversion length in a multi-layer slope, the slope surfaces were marked with blue colors before the rainfall. With the aid of the backlighting and the dye traces, it was simple to visualize the dyed streamlines with the water movement once a steady state was achieved. The lateral diversion length and breakthrough zone could also be measured directly. Figure 11 shows photographs of the dyed trace streamlines from single-layer and multi-layer

slopes with different tilt angles. The coarse layer is shown in deep gray in the pictures. For all of the inclined experimental cases, the dye traces were diverted downslope (referred to as capillary diversion, see Figure 12) and, in most cases, penetrated the coarse layer at different points (breakthrough). The slope at 7◦ with a breakthrough zone maintained a constant value when infiltration reached a steady state. In these three cases, a clear lateral flow region without infiltrate water passing through the coarse layer formed near the toe of the coarse layer. An amount of water entered the bottom layer through the breakthrough zone. Additionally, the breakthrough region was measured as the total length along the interface through which breakthrough was observed. This will be referred to as the piping in the subsequent discussion.

**Figure 12.** Side view of model experiments with blue dye traces during rainfall for the slopes in the inclined group. (**a**) case II (tilt angle α = 7◦); (**b**) case III (tilt angle α = 15◦); (**c**) case VI (tilt angle α = 21◦); (**d**) case VI (tilt angle α = 21◦); (**e**) case VII (tilt angle α = 21◦); (**f**) case VIII (tilt angle α = 21◦). The capillary barrier diversion in the multi-layer slope is shown by the movement of blue dye traces. The white arrow shows the approximate flow direction (case II and case VI).

Figure 13 shows the lengths of the three cases observed, and the capillary diversion length occurring upslope on the fine-coarse interface. The length of capillary diversion was measured from the initial point of the blue dye trace closest to the upper interface to the point where the dye first penetrated into the coarse layer. Figure 14 shows the VWC contour maps at 0.5 h, when failure occurred, rainfall stopped (t = 3.5 h), and the drying process stopped (t = 36 h) for case I, II, III, and VII. The initial VWC values were about 0.06 in the slopes.

In Figure 14a,b, comparing the VWC in the flat group at t = 0.5 h, the area above the interface reached a higher degree of saturation in the multi-layer slope than in the single-layer slope, since the capillary barrier prevented the water from infiltrating into the coarse layer, which made the slope more stable and caused a delay in the failure time. Comparing Figure 13, it was proven that lateral diversion occurred along the interface in the inclined slope, which resulted in a higher water content near the toe of the slope, and water could not infiltrate into the lower finer layer in Figure 13. This result may be associated with the sloping of the cover system [24]. The lateral diversion in the coarse layer was around 0.4 m, possibly because of the inclined angle, material properties, and rainfall intensity used in the test. Afterward, water could infiltrate into the deeper finer layer at the end of lateral diversion [26].

**Figure 13.** Effect of the tile angle on the capillary barrier diversion length in a multi-layer slope in case I, II, III, and IV (each case was performed three times for the same experimental conditions).

**Figure 14.** The VWC changes with elapsed time in case I, II, III, and IV (the VWC distribution after slope failure was verified by obtaining soil samples from various locations at the end of the test). (**a**) Case I. Multi-layer slope, 0 deg, 75 mm/h; (**b**) Case II. Multi-layer slope, 7 deg, 75 mm/h; (**c**) Case III. Multi-layer slope, 15 deg, 75 mm/h; (**d**) Case VII. Multi-layer slope, 21 deg, 75 mm/h.

3.3.2. Influence of the Tilt Angle on the Pore Pressure and VWC

To investigate the effects of the tilt angle, the pore pressure sensor (S-I) and VWC sensor (W-L) at the bottom were considered for comparison in different cases, as these are most crucial to the

knowledge on slope failure. Detailed in Figure 14 is the progressive build-up of VWC and pore water pressure throughout the experiments I, III, and IV, until failure occurred.

For Figure 15a,b, it is clear that an increased tilt angle has a drastic effect on the build-up of volumetric water content and pore water at the bottom of the slope, progressively resulting in quicker failure times as the tilt angle increases.

**Figure 15.** Time series of the VWC and pore pressure at the bottom of the slope for case I, III, and IV, with α ranging from 0◦ to 21◦ (multi-layer slope). (**a**) Time series of the VWC from initiation to failure; (**b**) time series of the pore pressure from initiation to failure.

Once the wetting front arrived, sensor spikes in measurements occurred. It took around 1.18, 0.78, and 0.62 h for the pressure head and VWC to approach the maximum value in case I, III, and IV, respectively. The peak value of VWC and pore pressure values were similar for each case ranging between 0.348 and 0.351 and −0.021 and −0.016 kPa, respectively, at times of failure. Due to the rapid progression of the wetting front at an increased tilt angle, the pore water pressure and VWC increased at a faster rate, as is evident in results showing that the α = 21◦ (case V) failed 47% sooner than α = 0◦ (case I) and 34% sooner than α = 15◦ (case III). These similar trends were exhibited in all experiments where the tilt angle was increased [27].

#### **4. Discussion**

#### *4.1. Mechanism of the Capillary Barrier*

A capillary barrier (Figure 16a) forms in unsaturated conditions when the hydraulic conductivity of finer soil is higher than that of coarse soil, which limits the downward infiltration of water due to the difference of the hydraulic conductivity. Figure 15b shows the relations between the hydraulic conductivity and suction of the two sands measured in the lab. The intersect of two hydraulic conductivity curves is 0.4 kPa. From the intersect to the higher suction (suction > 0.4 kPa), the hydraulic conductivity in the coarse layer is much lower than that in the fine layer, which caused the capillary barrier to form at the fine-coarse interface. Otherwise, suction at the interface decreased gradually after the wetting front arrived and was located to the left of the intersection, and the hydraulic conductivity in the coarse layer was higher than that of the fine one, which allowed the rainfall water to infiltrate into the next layer after capillary barrier breakthrough.

**Figure 16.** (**a**) Capillary barrier at the interface between different layers with different conductivities. (**b**) The relationship between the hydraulic conductivity and suction of silica No 1 and No 7.

#### *4.2. D Flow and Multi-Layer Slope Stability*

As is shown in Figure 17a, the blue color was used to stain the surface of the slope and to show the movement of rainfall water, and it was applied under the coarse layer to check that the capillary barrier still worked. At the beginning of case III, the water flow showed a vertical movement, mainly controlled by gravity, before the wetting front arrived at the interface. Secondly, the water flow could not go across the interface directly upward of the slope and showed a significant velocity component parallel to the inclined interface after the wetting front arrived [28].

**Figure 17.** (**a**) Side view of the diversion length and breakthrough area of the capillary barrier (the movement and the breakthrough zone are shown by the blue dye traces); (**b**) calculation of the capillary diversion length *L*.

Lateral diversion is essentially gravity-driven unsaturated drainage within the finer layer of a sloped capillary barrier [29,30]. The unsaturated hydraulic conductivity is the maximum as the VWC of the upper finer soil increases with depth, where lateral diversion is concentrated at the fine-coarse interface. Besides, the diverting water will increase the water content in the downslope direction with the influence of the inclined angle, which may result in failure of the barrier and then infiltrate into the next layer.

In Figure 17b, the maximum length of the capillary diversion can be calculated by [31]:

$$K(\psi) = \frac{\left\{1 - (a\psi)^{mn} \left[1 + (a\psi)^n\right]^{-m}\right\}^2}{\left[1 + (a\psi)^n\right]^{m/2}} \ast \text{K}\_{\mathfrak{s}}\,\,\,\,\tag{2}$$

$$L \le \tan(\beta) \frac{K(\psi)}{q} \left[ \frac{1}{\alpha} + \left( \left| \psi\_a \right| - \left| \psi\_w^\* \right| \right) \right] \tag{3}$$

where, ψ is the matric suction in the soil; ψ∗ *<sup>w</sup>* is the water entry value of the coarse layer; ψ*<sup>a</sup>* is the air entry value of the fine layer; β is the tilt angle of the slope (◦); *r* is the infiltration rate (mm/h); *K*(ψ) is the hydraulic conductivity of the fine layer; and *a*, *m*, and *n* are the fitting parameters.

Physical model experiments with a tilt angle α = 0◦, 7◦, 15◦, and 21◦ were performed three times in order to ensure the repeatability of the results. An important point to be noted is that every physical model was assumed to be identical, but there were slight differences in the model construction and preparation process, which could have resulted in some dissimilarity between experiments. Taking this fact into account, the result of the capillary diversion length at different tilt angles (Figure 18a), which produced similar results, demonstrated the repeatability of the experiment. The pictures of failure situations in case III (tilt angle α = 7◦) are described in detail in Figure 18b. Although the details of soil piping of each experiment were unique, the results of different experiments were consistent.

**Figure 18.** (**a**) Repeatability of the capillary diversion length (tilt angle α = 0◦, 7◦, 15◦, and 21◦, where each case was performed three times to ensure the repeatability of the results); (**b**) experiments of failure situations were conducted in case III under the same experimental conditions three times.

The capillary barrier effects caused by the inter-mediate coarser layer initially confines the rainwater infiltration in the uppermost soil. Depending on the characteristics of rainfall and the inclined angle condition, top finer soils are almost permanently unsaturated [32]. In some slope steeper than 40◦ at higher elevations, this may lead to instability of the top fine soil before capillary barrier breakthrough, when the coarse layer and bottom layer are still far from the point of saturation [33]. Depending on the slope angle and shear strength of soil, failure can occur at the bottom before complete saturation, while the pore water pressure is still negative. These phenomenon in natural slope seem to disprove the possibility of failure mechanism due to piping in this study.

In our experiment, the piping failure occurred above the impermeable layer after the loss of the capillary barrier. This is consistent with the field investigation of the development of soil pipe in the base of the fine soil [34]. However, the actual mechanism of soil piping in these slope that still remains unclear. Some mechanisms have been proposed by previous researches to explain this phenomenon such as internal erosion, flowslides and grain coarsening [34,35]. Internal erosion of the finer soil fraction driven by the seepage forces is thought to have played a significant influence on the slope failure [36,37]. It seems reasonable that the rainfall water caused the build-up of excess pore water pressure inside slope and the water inflow occurred during the rainfall condition, as the sharp increase of pore water pressure Δu measured at base of multi-layered slope (see Figure 10a). Further measurement such as the pore water pressure inside slope, internal displacements are still needed to improve the understanding of mechanism of failure of the slope [38].

#### **5. Conclusions**

Four groups of laboratory model experiments were performed to investigate the water movement, failure time, and modes in multi-layer and single-layer slopes caused by rainfall infiltration. In addition, the unsaturated hydraulic conductivity and SWCCs were also measured in the lab to clarify how the capillary barrier works under different conditions.

The results of hydraulic conductivity show that the unsaturated hydraulic conductivity *K* in the coarse layer is lower than that in the fine layer in a lower suction condition (0.4 kPa), which results in the development of the capillary barrier at the interface of the fine-coarse layer.

Different failure modes occurred in an inclined multi-layer slope and single-layer slope: sliding and piping failure. In the flat group, the capillary barrier was presented, which prevented the rainwater from infiltrating into the coarse layer for a while and caused a delay of the failure time. However, in the inclined group, the inclined intermediate coarser layer formed a capillary barrier, resulting in a significant amount of water being diverted to the downward slope and causing piping failure at the toe of the slope that resulted in earlier failure, which has a negative influence on the slope stability.

An increased tilt angle has a drastic effect on the capillary diversion length, in which more infiltrate water will be diverted to the downslope side and then infiltrate into the bottom of slope, resulting in quicker failure times.

The present study does not provide a model to be used in a specific site problem. Instead, the model is suitable for studies on hypothetical multi-layer hillsides to assess water movement patterns and general failure mechanisms. The results from such studies can prove useful in the development of an appropriate strategy for resolving problems in individual, site-specific multi-layer slopes.

In the present manuscript, we have mainly focused on the influence of the tilt angle and the effect of the capillary barrier on water movement and slope failure modes compared with a single-layer slope. It should be noted that in multi-layer slopes, different geometric characteristics, such as different heights, slope ratios, and fronts, also have an influence on the water infiltration and water content distribution in the slopes, significantly affecting slope failure initiation. Model experiments with different geometries of the shallow soil cover will be conducted in further research [39].

**Author Contributions:** Conceptualization, J.T. and U.T.; investigation, D.H.; methodology, S.T.; resources, J.X.; supervision, U.T.; Writing—original draft, J.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by a Japanese Government Scholarship for PhD studies of the first author, the Sichuan Science and Technology Program (2018JY0613), Grants-in-Aid for Scientific Research of Japan Society for the Promotion of Science (JSPS), and Research Fellowships of Japan Society for the Promotion of Science for Young Scientists.

**Acknowledgments:** The authors would like to thank the reviewers for constructive comments that improved the paper.

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

#### **References**

1. Santi, P.M. Debris-flow hazards and related phenomena: (Edited by matthias jakob and oldrich hungr). *Environ. Eng. Geosci.* **2007**, *13*, 75. [CrossRef]


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Exploring the Impact of Multitemporal DEM Data on the Susceptibility Mapping of Landslides**

### **Jiaying Li 1,2,3, Weidong Wang 1,2, Zheng Han 1,\*, Yange Li <sup>1</sup> and Guangqi Chen <sup>3</sup>**


Received: 5 March 2020; Accepted: 2 April 2020; Published: 6 April 2020

#### **Featured Application: Landslide susceptibility assessment and other geological disaster assessment issues.**

**Abstract:** Digital elevation models (DEMs) are fundamental data models used for susceptibility assessment of landslides. Due to landscape change and reshaping processes, a DEM can show obvious temporal variation and has a significant influence on assessment results. To explore the impact of DEM temporal variation on hazard susceptibility, the southern area of Sichuan province in China is selected as a study area. Multitemporal DEM data spanning over 17 years are collected and the topographic variation of the landscape in this area is investigated. Multitemporal susceptibility maps of landslides are subsequently generated using the widely accepted logistic regression model (LRM). A positive correlation between the topographic variation and landslide susceptibility that was supported by previous studies is quantitatively verified. The ratio of the number of landslides to the susceptibility level areas (RNA) in which the hazards occur is introduced. The RNA demonstrates a general decrease in the susceptibility level from 2000 to 2009, while the ratio of the decreased level is more than fifteen times greater than that of the ratio of the increased level. The impact of the multitemporal DEM on susceptibility mapping is demonstrated to be significant. As such, susceptibility assessments should use DEM data at the time of study.

**Keywords:** multitemporal DEM; control factors; susceptibility assessment; LRM; historical landslide events

#### **1. Introduction**

Geohazards are some of the most uncontrolled impacts on local and global economies, as well as on people's livelihoods. Almost nowhere on the planet is free from the damage of geohazards [1]. An essential component of predicting possible geohazard zones is the identification of an area that is vulnerable to future landslides [2,3]. Landslides are common geohazards that destroy local resources and environments. Therefore, susceptibility assessments of landslides have been widely investigated to improve their capability for use with these hazards [4,5].

Current studies commonly focus on two major issues regarding the susceptibility assessments of landslides—selection of evaluation indices and establishment of a rational assessment model. Relevant reviews [6–9] demonstrate that exiting studies have proposed various remarkable models for susceptibility assessment of landslides, such as logistic regression models (LRMs) [10]. Meanwhile, in recent decades, the rapid development of machine learning algorithms in susceptibility assessments has been seen over time. LRMs are some of the most classic and commonly used methods, which

have the advantages of minimal computation, high detection speed, and good adaptability [11]. LRMs have extensive applications in forecasting [12,13] and susceptibility assessment [14,15]. These studies provide a solid foundation for the susceptibility assessment of landslides, and advance the knowledge of susceptibility assessment with machine learning algorithms.

In addition to an assessment model, the rational selection of indices is quite important for susceptibility assessment of landslides. Generally, indices such as lithology, elevation, and topography are introduced as the controlling factors, which are strongly related to the susceptibility assessment of landslides [16–20]. The rational selection of these indices remains a subject of scientific debate in many studies because the importance of different indices may vary from case to case. A review of previous studies [21] provides a guideline for addressing this issue. However, the generalizability of these results is subject to certain limitations, while another essential problem regarding the temporal variation of indices and the influence of this on susceptibility is still confusing and requires more attention. The susceptibility result will be questionable if the data for one index is out of date owing to the absence of a current data source, which is the case in digital elevation models (DEMs). So far, however, few researchers have highlighted the impact of outdated data sources on the reliability of susceptibility results, and few previous studies have investigated the mechanisms involved.

A DEM is the digital representation of a terrain surface. These models have been widely applied in geohazard planning, terrain surface analysis, and other fields [22]. Generating a DEM generally involves data from different sources [23], including global multiresolution topography (GMRT) models, shuttle radar topography mission digital elevation models, and advanced space-borne thermal emission and reflection radiometer global digital elevation models (ASTER GDEMs).

Among the abovementioned indices, elevation, topography, and slope are closely related to the DEM data [24]. In the majority of the previous studies, it has been suggested that DEM data should be consistent with the study date [25–29]. However, because of landscape changes and reshaping processes (e.g., earthquakes and engineering construction), the DEM data usually show obvious temporal variation. An example of significant temporal variation of DEM data was shown by Cucchiaro et al., who substantiated that the DEM difference of the eastern Italian Alps within one month is up to 0.14 m [30]. Pineux et al. demonstrated that the changes experienced by a DEM over time are obvious and unpredictable [31]. DEM data have significant influence on susceptibility assessments of landslides, mainly through influencing factors such as elevation, topography, and slope. Elevation is one of the most influential factors controlling landslide occurrences in a study area, while topography and slope are also important in susceptibility assessments [18,29,32]. However, studies have not considered the obvious temporal variation in DEMs.

In the present paper, multitemporal DEMs area are obtained for the study and the elevation changes demonstrated in the DEMs are analyzed. The DEM data for three years (1992, 2000, and 2009) are used to analyze factors affected by a DEM. The susceptibility results for the three years are evaluated using the different DEM-dependent factors and the same DEM-independent factors. Quantifying the index, the susceptibility assessments are then obtained using a generally accepted LRM model. Additionally, the maps of susceptibility levels in the study area are obtained within the geographic information system (GIS). The influence of DEM data at different times on a susceptibility assessment and the influence of the specific values of the elevation differences on the assessment level are obtained. Meanwhile, the impact of historical landslide events on the susceptibility assessment is analyzed. The study aims to achieve a better evaluation of the impact of multitemporal factors on the susceptibility assessment of landslides.

#### **2. Study Area**

#### *2.1. Study Area and Landslide Data*

The southern area of Sichuan province in China is a region with frequently occurring landslides. The study area covers the three districts of Panzhihua, Liangshan, and Ya'an (Figure 1). The area is

located at a latitude of 26◦03 –30◦56 north and a longitude of 100◦03 –103◦52 east. The area also belongs to the Sichuan Basin, adjacent to the Qinghai–Tibet Plateau in the west. It has a complex topography mainly composed of plains, hills, mountains, and plateaus; furthermore, the relative elevation is more than 3200 m.

**Figure 1.** Location of the study area: (**a**) Sichuan province in China; (**b**) the study area; (**c**) an aerial photograph of the study area.

The climate in the study area is a typically humid subtropical monsoonal climate that is hot in summer and mild and humid in winter, with plentiful rainfall and mist. The vertical climate differs greatly, and the annual average temperature distribution gradually increases from northwest to southeast. The average temperature ranges from 5.3 to 15.7 ◦C, with annual average rainfall ranging from 1500 to 1800 mm. The study area is, therefore, extremely rainy. The rainfall is mainly concentrated from June to October, which accounts for approximately 70%–75% of the annual rainfall. The maximum daily precipitation reaches 300–500 mm. Meanwhile, part of the study area is located in the Himalaya earthquake zone, which suffers from permafrost, avalanches, and landslides. The stratum of the study area is complete (from Archean to Quaternary). The dominant lithology is sedimentary rock, majorly consisting of dolomite, limestone, siliceous rock, shale, and sandstone [33]. The surface water system is developed, and there are many tributaries of the Yangtze River. The groundwater is distributed widely and buried shallowly, which is affected significantly by the rainfall and landform [34].

The Chinese Geological Environment Monitoring Institute collected the landslides information in this region through the China Geological Survey, including locations and occurrence times, and released it to public in its Bulletin of National Geological Hazards [35]. The annual total number of historical landslide events is shown in Figure 2. More than 85% of the historical landslide events were induced by natural environmental factors, such as rainfall and earthquakes, while only a few events were induced by human factors, such as mining and slope cutting. As the source data in the Bulletin of National Geological Hazards are not distinguished by landslide type, the general landslides analyzed in the study are of various types.

**Figure 2.** Number of historical landslide events in various years.

In this study, the total number of historical landslide events from 2000 to 2016 and the number of those hazards in various years after 2004 are obtained. However, we can only obtain the total number of historical landslide events for the period between 2000 and 2004. Thus, to better analyze the impact of a DEM on susceptibility assessment, we divide the historical landslide events into two broad types: the landslide events during the period ranging from 2000 to 2009, and the period ranging from 2010 to 2016. The historical landslide events in the study area from 2000 to 2009 and from 2010 to 2016 (Figure 3) can be obtained from the China Geological Environment Information database (http://www.cigem.cn). A total of 341 historical landslide events from 2000 to 2009 and 218 historical landslide events from 2010 to 2016 were recorded.

**Figure 3.** Locations of historical landslide events.

#### *2.2. Influencing Factors*

Selecting the influencing factors is a key step in a landslide susceptibility assessment [36,37]; however, the causes of landslides are complex. To date, hundreds of influencing factors have been identified that are potentially important to susceptibility assessments [38]. For example, areas with high elevation and steep slopes are highly prone to landslides [39]. Lithology is also related to landslide occurrences, because different lithologies can withstand different levels of triggering factors. The distances to structure lines, rivers, and roads have important impacts on the spread and size of landslides in the study area [40]. Pourghasemi et al. investigated global susceptibility during the period of 2005–2016 [28]. Nearly 100 factors were summarized to reveal the commonly used factors. In accordance with the relevant studies [29,41,42] and the study area, we select eight influencing factors: elevation, topography, slope, lithology, distance to a structure line, distance to a river, average annual rainfall, and distances to roads.

Table 1 lists all the DEM data for the study area and their sources. The Open Topography Facility provided the GMRT data, which is hosted at the San Diego Supercomputer Center, University of California San Diego. This facility has built a strong cyberinfrastructure framework for managing and processing high-resolution topography data from light detection and ranging (LiDAR) (http://opentopo.sdsc.edu/datasets?listAll=true). Meanwhile, the Shuttle Radar Topography Mission (SRTM) DEM and Global Digital Elevation Model (GDEM) DEM were provided by the Geospatial Data Cloud site of the Computer Network Information Center in the Chinese Academy of Sciences (http://www.gscloud.cn/). The DEM data at the same resolution for 1992, 2000, and 2009 were then obtained using a resampling tool in the GIS environment (Figure 4). There are three resampling methods in the resampling tool, namely nearest neighbor, bilinear interpolation, and cubic convolution methods. Nearest neighbor is selected in the present study because the method is simple, fast, and applicable.

**Table 1.** Digital elevation model (DEM) data used in this study and their sources.


\* GMRT is short for Global Multi-Resolution Topography. \*\* SRTM is short for Shuttle Radar Topography Mission. \*\*\* GDEM is short for Global Digital Elevation Model.

**Figure 4.** DEM data for the study area in (**a**) 1992, (**b**) 2000, and (**c**) 2009.

The DEM data for 1992, 2000, and 2009 were used to analyze factors affected by the DEM. Meanwhile, the data of other factors, namely lithology, distance to the structure line, distance to the river, and average annual rainfall, remained the same and were provided by China Railway Number 4 Engineering Group Co., Ltd. (745 Heping Road, Wuhan, Hubei, China), and the Roads and Traffic Authority of China. Because of the diversity of lithological layers in the study area, the lithology is divided into four groups in Table 2 [43]. Thus, the zoning maps of the four factors were obtained (Figure 5).


**Table 2.** Lithology and structure in the study area.

**Figure 5.** Zoning map of (**a**) lithology, (**b**) distance to a structure line, (**c**) distance to a river, (**d**) and average annual rainfall.

#### **3. Methodology**

#### *3.1. Flowchart*

The research methodologies applied in the present study are as shown in Figure 6. The flowchart consists of four major steps, as follows:


**Figure 6.** Methodology of research applied in this study. LRM—logistic regression model.

#### *3.2. Confidence Interval*

In the present study, because of the noise in the DEM data reducing the accuracy of data analysis significantly, not all DEM data are reliable. Therefore, the noise needs to be eliminated to determine the reliability of the DEM data. The confidence level of the DEM data is assessed to improve the accuracy of data analysis.

The confidence interval is a commonly-used interval estimation method used for sample statistics, which shows the confidence probability of the measured parameter value [42,44]. The elevation difference is considered as the sample data in the present study. The average value of elevation difference is μ, and the standard deviation is σ. The confidence probability is obtained as follows:

$$\Pr(c\_1 \le \mu < c\_2) = 1 - \alpha \tag{1}$$

where α is the significance level and the interval (*c*1, *c*2) is the confidence interval. Therefore, the confidence interval of the average value is (μ − σ*Z*α/2, μ + σ*Z*α/2), where = *Z*α/2 is the corresponding standard score. In general, the confidence probability in the literature [42,45,46] is 90% or 95%, and *Z*α/2 is 1.645 or 1.96, respectively.

#### *3.3. LRM*

The LRM is a statistical model used to predict the probability of a categorical occurrence using one or more independent variables [47]. The purpose of LRM is to obtain the relationship between a dependent variable and multiple independent variables that have been identified.

In the present study, the independent variables are the influencing indices of landslides, and the dependent variable is the probability of the landslide occurring. By transforming the universal formula of LRM, *Y* is obtained:

$$Y = \mathbb{C}\_0 + \sum\_{i=1}^{n} \mathbb{C}\_i I\_i \tag{2}$$

where *C*<sup>0</sup> is the LRM constant coefficient, *Ci* is the LRM coefficient, *Ii* is the landslide index, Y = ln(P/(1 − P)), and *P* is the probability of the landslide occurring. The LRM coefficient *Ci* and the LRM formula can be obtained using Statistical Product and Service Solutions (SPSS) software.

#### *3.4. Ratio of Number of Landslides to Area (RNA)*

*RNA* denotes the ratio of the number of historical landslide events to the area of the susceptibility levels at which historical landslide events are located. The RNA of the level *i* is obtained with the following equation:

$$RNA\_i = n\_i / A\_i \tag{3}$$

where *RNAi* is the number of landslides per area at a particular level *i*; *ni* is the number of historical landslide events at level *i*; and *Ai* is the area of level *i* in the assessment.

#### **4. Results**

#### *4.1. Elimination of Noise*

As mentioned above, the DEM has a significant effect on susceptibility assessment through three indices, namely elevation, topography, and slope. The changes in the DEM from 1992 to 2000 and from 2000 to 2009 are obtained to explore the impact of the temporal variation of DEM data on susceptibility assessment. Meanwhile, the average values and standard deviation of the elevation differences are obtained from the topography changes revealed by the DEM data. The confidence probability is considered as 95% in this study. The corresponding standard scores and confidence intervals are, therefore, calculated (Table 3). After the noise is eliminated, the optimized map of topography changes from 1992 to 2000 and from 2000 to 2009 are obtained (Figure 7).

**Figure 7.** Topography changes revealed by the DEM data: (**a**) from 1992 to 2000; (**b**) from 2000 to 2009.


**Table 3.** Confidence interval of elevation differences for the periods ranging from 1992 to 2000 and from 2000 to 2009.

#### *4.2. Susceptibility Assessments of Landslides*

The susceptibility assessment maps of the study area in 1992, 2000, and 2009 are obtained in the GIS environment using LRM (Figure 8 and Table 4). The susceptibility is classified into four levels based on the natural breaks classification (NBC) method, which is based on natural groupings inherent in data. The NBC identifies groups of similar values and maximizes the differences between classes. The features are divided into classes whose boundaries are set based on relatively big differences in the data values. Here, level I, level II, level III, and level IV denote low, moderate, high, and very high susceptibility, respectively. The area percentages of various susceptibility levels in the study are shown in Table 4.

**Figure 8.** Susceptibility assessment maps for (**a**) 1992, (**b**) 2000, and (**c**) 2009.


**Table 4.** Susceptibility assessment of landslides using LRM scores.

As shown in Figure 8 and Table 4, the results of the susceptibility assessment vary from year to year. The areas with low and moderate susceptibility (levels I and II, respectively) in 1992 are the largest, while the area with very high susceptibility (level IV) is the largest in 2000. The rapid decrease in the level IV area from 2000 to 2009 is noticeable. One possible explanation for the decrease is that the slopes where landslides occurred became stable in the short term, owing to the excessive landslides that occurred from 2000 to 2009; thus, the susceptibility level of the study area decreased after the hazards.

The susceptibility maps are obtained using different DEM-dependent factors, while the DEM-independent factors are constant in the assessment. Comparing these results illustrates the significant influence of the DEM data from different years on the susceptibility assessment.

#### *4.3. Assessment Levels of Historical Landslide Events*

Figure 8 and Table 5 show the susceptibility assessment levels at which historical landslide events occurred from 2010 to 2016.


**Table 5.** Historical landslide event levels.

\* RNA denotes the ratio of number of landslides to area.

It can be clearly seen in Table 5 that the combined number of historical landslide events designated as levels III and IV is 106, 140, and 154, respectively, in 1992, 2000, and 2009. The RNA for level IV in 2009 is 5.44 <sup>×</sup> 10−3/km2, which is the largest among the RNA values over the three years. The most ideal assessment of landslides would be most landslides occurring in level IV areas. There is a significant positive correlation between the RNA and the accuracy of the assessment. Therefore, the results of the susceptibility assessment using the DEM for 2009 are more accurate than the results using the DEM for 1992 and 2000.

It is apparent from Table 5 that the small change in elevation over time has a great influence on the susceptibility assessment of landslides. However, in previous research, the survey data are not explained or the data are significantly different from the assessment time.

#### *4.4. Di*ff*erences between the Assessment Levels*

Comparisons are performed using the susceptibility assessments of landslides in the study area from different periods. Figure 9 and Table 6 compare the differences between the assessment results for 1992 and 2000 and between those for 2000 and 2009, showing that the area where the level difference is zero is the largest area. Meanwhile, the area with increased susceptibility is much larger than the area with decreased susceptibility in the period ranging from 1992 to 2000, which is contrary to the results for the period ranging from 2000 to 2009.



**Figure 9.** Maps of the differences between the assessment levels: (**a**) difference between 1992 and 2000; (**b**) difference between 2000 and 2009.

#### **5. Discussion**

#### *5.1. Comprehensive Comparison of Assessment Results for 1992, 2000, and 2009*

We compare the areas and RNA values of various assessment levels for the years 1992, 2000, and 2009 (Figure 10a,b), and the level differences from 1992 to 2000 and from 2000 to 2009 (Figure 10c). Figure 10 is quite revealing in several ways. First, the RNA for level IV in 2009 is larger than that in 1992 or 2000. Because of the positive correlation between the RNA and the accuracy of the assessment, the accuracy of the assessment is further exemplified in studies using the DEM for 2009. As Figure 10c shows, there is a significant difference between the two results. The area with increased susceptibility in the period ranging from 1992 to 2000 is larger than that in in the period ranging from 2000 to 2009. The observed increase in the area with increased susceptibility could be attributed to the frequent landslides from 2000 to 2009 (e.g., the various landslides caused by the 2008 earthquake).

**Figure 10.** Comprehensive comparison of the percentages of (**a**) the assessment levels, (**b**) the RNA values of various assessment levels, and (**c**) the level differences.

#### *5.2. Impact of Elevation Di*ff*erence on the Susceptibility Assessment*

In addition to the above analysis, we analyze the impact of elevation differences on the susceptibility assessment and the relationships between the elevation differences and the differences in assessment levels. Therefore, the areas with varying level and elevation differences are presented (Table 7).

**Table 7.** Numbers of historical landslide events at various assessment levels from 1992 to 2000 and from 2000 to 2009.


Table 7 demonstrates that the decrease in elevation results in a decrease in the susceptibility assessment level. There is a positive correlation between the elevation difference and the difference in the assessment level. Regardless of the level difference, the elevation difference is much larger in the area ranging from −20 to 20 m than in other areas.

#### *5.3. Impact of Historical Landslide Events on Susceptibility Assessment*

Most of the factors influencing susceptibility assessments have been explored in several studies [48–50]. However, much of the historical research overlooks the impact of historical landslide events on the assessment [51]. In contrast, because the slopes where landslides occur become stable in the short term following these events, the susceptibility of the area will consequently decrease. Therefore, historical landslide events have a significant effect on susceptibility assessments.

The numbers and percentages of historical landslide events at various assessment levels from 2000 to 2009 and the RNA values of the various assessment levels are investigated, allowing the change in the assessment levels of the historical landslide events to be obtained (Table 8 and Figure 11).

In Figure 11, the circles represent the percentages of historical landslide events or the RNA values at various assessment levels. It can be seen from Figure 11 that over one-third of the historical landslide events (38.71%) occur at a decreasing assessment level. A total of 54.24% of the historical landslide events are within the same assessment level. If the assessment level area is larger, then there may be more historical landslide events within that assessment level. To rule out the possible influence of area size, the RNA is considered to be the most important factor. The RNA at decreased levels (81.31 <sup>×</sup> <sup>10</sup>−3/km2) is more than fifteen times greater than the RNA at increased levels (5.25 <sup>×</sup> <sup>10</sup>−3/km2). Previous studies [52,53] do not take into account the impact of historical landslide events; therefore, the above results can establish the importance of historical landslide events for susceptibility assessments of landslides.


**Table 8.** Numbers and percentages of historical landslide events at various assessment levels and RNA values of various assessment levels.

**Figure 11.** Changes in the assessment levels of historical landslide events from 2000 to 2009.

#### *5.4. Advantages and Limitations*

The elevation, topography, and slope—depending on the DEM—are usually the main factors affecting susceptibility [23]. The existing studies assume that the influencing factors for susceptibility are constant over time [15,54]. However, the DEM data always change significantly over time owing to landscape reshaping, such as landslides and engineering construction [30,31,55]. A susceptibility assessment is, therefore, closely related to the data source [56], and as such the susceptibility maps of the study area are not consistent.

Despite this, the existing studies have analyzed the susceptibility of landslides using a single time-related data source without considering whether the date of the data source is consistent with the study date [57,58]. Therefore, it is necessary to understand the temporal evolution of influencing factors (e.g., topography) and the influencing mechanisms of the factors on the susceptibility assessment [59,60]. The changes in influencing factors are significant in earthquake- and landslide-prone areas [61], such as the study area. Therefore, more research regarding the impacts of factors on the susceptibility of landslides needs to be undertaken. It is better not to evaluate susceptibility simply using constant factors.

The present study performs a preliminarily analysis of the impact of multitemporal DEM on susceptibility and shows that DEM variation has an impact on the susceptibility assessment using DEM-dependent factors. The analysis undertaken here can extend the knowledge regarding the impacts of multitemporal data on susceptibility and can aid preliminary studies on the mechanism causing the influence of factor variation.

DEM variation closely relates to local geotechnical properties, which change the susceptibility assessment of landslides. However, the reasons for DEM variation are debated. The mechanism explaining the consequent influences of geotechnical properties on DEM variation remains a scientific challenge. In addition, considering the large spatial extent and period (about 77,000 km2, 17 years) in the present study, changes between the three susceptibility maps could not be completely explained by DEM variation. There are also data limitations and difficulties in analyzing the interactions of various factors.

#### **6. Conclusion**

In the present paper, DEM data for 1992, 2000, and 2009 are obtained to evaluate the susceptibility of landslides. Based on the DEM-dependent factors and the same DEM-independent factors, the susceptibility results are evaluated using LRM.

From the assessment, we find that the results are different by using DEM data from different times. The rapid decrease in the area of level IV from 2000 to 2009 is noticeable. The RNA for level IV in 2009 is larger than in 1992 and 2000. The area of increased susceptibility is much larger than that of decreased susceptibility based on the results of the assessment level differences from 1992 to 2000, which differ from the results from 2000 to 2009.

From the assessment results, we conclude that the DEM data have an impact on the susceptibility assessment in the study area. It is also worth noting that the influences of the specific values of the elevation differences on the assessment results are obtained. Meanwhile, the influence of historical landslide events on the susceptibility assessment is obtained by analyzing the hazard data and the differences between the assessment results. The most obvious finding to emerge from the analysis is that the assessment level of the area with historical landslide events decreases. The present study goes some way towards enhancing our understanding of the impacts and the mechanism of multitemporal factors on susceptibility.

**Author Contributions:** W.W., Z.H., and G.C. designed the study. J.L. and Y.L. performed the in situ investigation. Z.H. and J.L. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the National Key R&D Program of China (Grant No. 2018YFC1505401), the National Natural Science Foundation of China (Grant No. 51478483, W.D. Wang; and Grant No. 41702310, Z. Han), and the Natural Science Foundation of Hunan (Grant No. 2018JJ3638) These financial supports are gratefully acknowledged. Thanks to the Scholarship of Central South University for Undergraduate Students Studying Abroad.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Spatial Proximity-Based Geographically Weighted Regression Model for Landslide Susceptibility Assessment: A Case Study of Qingchuan Area, China**

#### **Yange Li 1, Xintong Liu 1, Zheng Han 1,\* and Jie Dou <sup>2</sup>**


Received: 27 December 2019; Accepted: 4 February 2020; Published: 7 February 2020

#### **Featured Application: Landslide susceptibility assessment and other geological disaster assessment issues.**

**Abstract:** Landslides pose a serious threat to the safety of human life and property in mountainous regions. Susceptibility assessment for landslides is critical in landslide management strategy. Recent studies indicate that the traditional assessment models in many previous studies commonly assume a fixed relationship between influencing factors and landslide occurrence within an area, resulting in an inadequate evaluation for the local landslides susceptibility. To address this issue, in this paper we propose a spatial proximity-based geographically weighted regression (S-GWR) model considering spatial non-stationarity of landslide data for assessing the landslide susceptibility. Spatial proximity is the basic input condition for the proposed S-GWR model. The challenge lies in defining the spatial proximity expression that shows the geographical features of landslides and therefore affects the model ability of S-GWR. Our solution chooses the slope unit as spatial adjacency, rather than the grid unit in DTM. The multicollinearity between landslide influencing factors is then eliminated through variance inflation factor (VIF) method and principal component analysis (PCA). The proposed model is subsequently validated by using data in Qingchuan County, southwestern China. Spatial non-stationary is identified for landslide data. A comparison with grid unit and four traditional evaluation models is conducted. Validation results using the area under the ROC (receiver operating characteristic) curve and success rate curve indicate that the spatial proximity-based GWR model with slope unit has the highest predictive accuracy (0.859 and 0.850 respectively).

**Keywords:** landslide susceptibility assessment; geographically weighted regression; spatial non-stationary; spatial proximity; slope unit

#### **1. Introduction**

Landslides are catastrophic natural hazards frequently posing risks to the major societal, economic, and environmental on an international scale [1]. According to the report from EM-DAT [2], 21,412 landslides occurred worldwide between 1900 and 2014, resulting in 38,521,499 fatalities, with 7,229,487,068 people affected and total direct economic losses exceeding \$2.7 trillion.

Landslide susceptibility assessment has long been recognized as a useful tool for landslide hazard management through land use planning and better decision making in landslide prone areas [3]. It is generally based on heuristic, statistical, or deterministic models [4–8]. Heuristic models are subjective and much susceptible to the expectation of the results [9,10]. Deterministic models have

been reported with higher accuracy, but are limited by the difficulty of obtaining detailed landslide database [11,12]. Statistical models are the most widely used models due to their simplicity and high efficiency [13,14]. Many remarkable studies on the above aspects have been made, laying a solid foundation for landslide susceptibility mapping. However, in general, most of the previous studies consider the relationship between triggering factors and landslide occurrence as a fixed effect within an area, whereas different degrees of parameter influence may occur, such that, with the change of location, the effect of parameters can be consequently changed. The uncertainties due to this varied relationship remain a scientific challenge.

The second law of geography suggests that there exists variability over space of a given relationship between variables widely in spatial data [15], which is the so-called spatial non-stationarity. In view of that landslide susceptibility assessment is heavily based on spatial data, the relationship between influencing factors and landslide susceptibility may also have the characteristics of spatial non-stationarity. Previous studies, e.g., [16], have also suggested that the effective parameters in the occurrence of a natural disaster phenomenon do not have the same importance in different parts of an area. The existence of spatial non-stationarity indicates that average relationships fitted to the whole study area of traditional models might be inappropriate since they do not accurately fit local conditions [17]. This spatial non-stationarity characteristics in the data pose difficulties in landslide susceptibility assessment based on the traditional models.

Geographically weighted regression (GWR), the most popular local regression format, shows great capability in dealing with spatial non-stationary relationships [18]. It allows the relationship between dependent and independent variables to vary over space, as well as that regression coefficients in the model are calculated for each spatial zone [19]. This method has been applied in various fields of study such as social economics, geography, and meteorology [20–22]. However, previous studies applying GWR model for the assessment of geological hazard susceptibility have not yet been reported.

One difficulty limiting the application of GWR model in landslide susceptibility assessment is spatial proximity. It is the basic input condition and core problem for GWR model, and the issue regarding an adequate expression for the spatial proximity at different locations directly affects the modeling ability of GWR model [23,24]. Spatial proximity is the distance relationship between two units in space, and the closer the distance, the greater the impact. The key to determining the spatial proximity is segmenting the study area into map units to effectively express the spatial adjacency relationship between landslide data. The relationship should satisfy the requirements of GWR for good internal homogeneity and between-units heterogeneity. The commonly used segmentation methods in previous studies relating to GWR model can be categorized into two major kinds, i.e., administrative units [25] and grid units [26]. Administrative unit is mostly used for social and economic issues, and its segmentation does not accord with the neighborhood characteristics of landslide data. As such, administrative unit is rarely used in geological hazard assessment. Grid unit is a popular mapping unit for susceptibility assessment since it is easily accessible, but it is not associated with geological environments. Slope unit is a relative new mapping unit for evaluating landslide susceptibility, which is generated according to hydrology theory and is the watershed area defined by drainage lines (valley lines) and water divide lines (ridge lines) [27]. It is the basic topographical unit of geological hazard occurrence. Slope unit has higher internal homogeneity and between-unit heterogeneity than grid unit. It is closely related to geological environment conditions. In this sense, slope unit provides an alternative solution for spatial proximity expression of the GWR model for landslide susceptibility assessment.

Two other key issues of GWR model are the multicollinearity elimination and the kernel function establishment. Previous studies, e.g., [28], indicated that GWR is highly susceptible to the effects of multicollinearity between explanatory variables, and collinearity among pairs of explanatory variables or multicollinearity among more than two variables often lead to problems such as parameter estimate instability and unintuitive parameter signs. These problems remain significant owing to the complicated conditions of landslide posing a high possibility of correlations between explanatory variables. Kernel function is based on the distances between observations and calibration units to place emphasis on observations that are closer in space [28]. The selection of kernel function type and the determination of its bandwidth are crucial to the spatial proximity modeling of GWR.

In this study, we attempt to propose a spatial proximity based on geographically weighted regression (S-GWR) model for landslide susceptibility assessment. The presented model resolves the spatial non-stationarity of landslide susceptibility assessment with GWR model. Firstly, we generate slope units to establish spatial adjacency. Then, variance inflation factor (VIF) method [29] and principal component analysis (PCA) method [30] were adopted to eliminate multi-collinearity, and kernel function was determined according to the characteristics of landslide data. Finally, we chose Qingchuan County, Sichuan Province, China, as the study area to validate the applicability of the model, and further compared the established model with the grid-unit GWR model and other evaluation models.

#### **2. Study Area**

The study area is the Qingchuan County in the transitional region between the Sichuan Basin and the Western Sichuan Plateau. This area has long been recognized as one of the most landslide-prone areas of China [31]. It locates between 32◦12 ~32◦56 N in latitude and 104◦36 ~105◦38 E in longitude, covering a total area of 3217 km2. The minimum elevation of the Qingchuan County is approximately 500 m and the maximum is 3820 m, characterized by northwestern part with higher elevation than the southeastern. Slope gradient reaches a maximum of about 80◦, with a mean value of 38◦.

The tectonics and geological settings in the area are complex. Because of the neotectonics, soft-lithology and hard-lithology usually appears alternately. There are about eight types of lithological outcrops throughout the study region (as shown in Figure 1), including the sedimentary rock (limestones, sandstone, and conglomerate) from Cambrian to Jurassic age, magmatic (granite), metamorphic rock (shales, schists, gneiss) from Cambrian to Jurassic age and Quaternary loess unconsolidated sedimentary. Two main active faults cross the area: the Pingwu–Qingchuan fault located in the north and crossing the whole territory, and the Yingxiu–Beichuan fracture which belongs to the Longmenshan fault belt, is a thrust fault 60◦–70◦ NW dipping. Bailong river, Qingzhu river and Qiaozhuang river are distributed in the area. The discharges of the three rivers are measured approximately 525, 30, and 40 m3s−1, respectively, serving as the main channel for atmospheric precipitation and groundwater drainage.

**Figure 1.** Lithological map of Qingchuan area and locations of landslide points.

#### **3. Methodology**

#### *3.1. Flowchart of Research*

The methodologies used in this study are as shown in Figure 2. The flowchart consists of three major steps.

**Figure 2.** Methodology of research applied in this study.

The first step is dataset preparation, including the preparation of basic data and the determination and extraction of landslide-related causal factors (LRCFs);

The second step is the establishment of the landslide susceptibility assessment model, which is the most critical step. Firstly, slope units are generated to establish spatial adjacency of landslide data. After obtaining the data in the slope units, the multicollinearity between data is eliminated by variance inflation factor method (VIF) and principal component analysis (PCA). Then, the spatial kernel function of GWR model is determined, and finally the GWR model based on slope unit (S-GWR) is established for landslide susceptibility evaluation.

The last step is model validation, including applying the S-GWR model to the actual area and generating the landslide susceptibility. The results are then used for validation and compared with the grid-unit GWR model (G-GWR) and other evaluation models including ordinary least squares regression (OLS), artificial neural network (ANN), and information models (I).

#### *3.2. Dataset Preparation*

The basic data utilized in the study include the digital elevation model (DEM) of the study area with a spatial resolution of 10 m; a geological map with the scale of 1:100,000; aerial photographs taken by Ministry of Land and Resources. After field identification and interpretation of the color aerial photographs, 973 landslide points are identified, and the distribution is as shown in Figure 1.

Owing to that a variety of factors make a contribution to the landslide occurrence, an adequate selection of the effective parameters is essential to identify areas prone to landslide. We choose these factors in accordance with the well-illustrated criterion in many previous studies (as listed in Table 1). In detail, firstly, the actual conditions of research area are considered. The study area has a complex topography, so topographic factors, such as slope and slope height, provide terrain conditions for landslides. The alternating appearance of soft and hard lithology provides a material basis for the development of landslides. Besides, there are two large faults across the study area, and the landslide frequency may increase as the distance to a major or minor fault decreased [26]. There are also many river systems in the study area. In the dense-river-areas, with the increase of catchment area and reduction of erosion basis, the shearing force of river flow is enhanced, which gives rise to the steep gorges in the downstream and provides premise conditions for the growth of geological disasters.



Based on the analysis of the research area, the factors commonly used in previous studies on landslide susceptibility are further referenced (Table 1), and 7 variables (elevation, slope, slope height, aspect, distance-to-stream, distance-to-fault, and lithology) are finally identified for modeling. All these factors are processed within ArcGIS. Slope, elevation, slope height, and aspect are calculated using DEM. All data are directly obtained by ArcGIS without data transformation. Lithology map is digitized from the existing geological map. Lithology map is digitized from the existing geological map. Distance-to-fault is obtained by calculating the distance to the nearest fault, and distance-to-stream is the distance to the nearest stream.

#### *3.3. Establishment of Spatial Proximity*

Spatial proximity is the distance relationship between two units in space. GWR model establishes a model for each unit according to the spatial proximity, and the areas within a unit are considered homogeneous. Therefore, it is essential to determine the spatial adjacency relationship between units, i.e., the boundary of map units. Administrative boundaries [25] and grid boundaries [26] were commonly used as the spatial proximity expressions of GWR in previous studies, but they are inconsistent with the neighborhood characteristics of landslides. These boundaries could not perform well to express the heterogeneity between units and the homogeneity within units, affecting the modeling ability of GWR model.

Therefore, we incorporate the methodology of slope units to express spatial proximity. Slope unit is the basic unit of geological disasters. It divides the terrain into mapping units with similar hydrological and geomorphological conditions, and is shaped by similar processes occurring in the natural landscape under the same geo-environmental conditions [37]. A slope unit division map is formed by the GIS-based hydrologic analysis tool [38]. Firstly, reverse DEM is generated by subtracting the elevation value from the highest elevation value in each unit. Secondly, fill the DEM and reverse DEM, and the flow direction can be obtained by these filled DEMs. Then, by setting the minimum number of cells that flow to the calculating point, the watershed can be calculated. Eventually, by combining the watershed by DEM and the watershed by reverse DEM, slope units can be obtained.

#### *3.4. Elimination of Multicollinearity*

Collinearity among pairs of explanatory variables or multicollinearity among more than two variables in regression analysis is known to cause problems such as parameter estimate instability, unintuitive parameter signs, high coefficient of determination (R2) diagnostics despite few or no significant parameters, and others [28]. Due to the huge and complex landslide data, the complicated relationship between geological and topographical factors should be considered. Therefore, the multicollinearity elimination in landslide susceptibility assessment is very important. Correlations in GWR parameters, both within a set of local parameter estimates for all locations (global multicollinearity) and among different parameter estimates at each location (local multicollinearity), are the symptom of multicollinearity among explanatory variables [28].

Global multicollinearity of the entire area can be easily distinguished by the variance inflation factor (VIF) method [29], which measures the effect of collinearity on the estimated variance of a regression coefficient. Local multicollinearity is the linear dependencies in the design matrix of local regression model. Principal component analysis (PCA) [30] is adopted to eliminate local multicollinearity since the diagnosis of local multicollinearity is very complicated. Through data transformation and processing, the influencing factors of landslide susceptibility can be grouped into less integrated factors, which not only maintains the main information of original factors, but also weakens the correlation among them. First obtained initial factor loading matrix. There are two main principles for selecting the principal component: (1) the principal component eigenvalue is greater than 1; (2) the cumulative contribution rate of principal components reaches 80%. Then, the factor rotation is performed so that the obtained factors have clear professional interpretation significance. Quartimax method [39], a common factor rotation method, was used for factor rotation.

#### *3.5. GWR Modeling*

GWR model extends the ordinary least squares regression (OLS) [40] by weighting the spatial dependence [41]. Based on the established spatial proximity, the coefficients of the model are estimated for each unit, and the value and symbol of the coefficients vary at different units [16]. This model is in the form of Equation (1)

$$y\_i = \beta 0(\mu\_i, v\_i) + \sum\_{j=1}^p \beta\_{ij} (\mu\_i, v\_i) x\_{ij} + \varepsilon\_i \tag{1}$$

where (μ*i*, *vi*) represents the coordinates of an *i*th unit in space and *p* is the number of independent factors. *xij* is the *j*th independent variable of the *i*th unit. β<sup>0</sup> (μ*i*, *vi*) is the intercept parameter in position *i* and ε*<sup>i</sup>* is the random error. β*ij* (μ*i*, *vi*) is the local regression coefficient for the *j*th explanatory variable in position *i*, which varies with the change of spatial position and is a very important parameter for the embodiment of spatial non-stationary.

The establishment of the weight kernel function is an important step of GWR, which is used to determine the scope and degree of spatial dependence. The establishment process includes the selection of the type of kernel function and the determination of its bandwidth, and previous studies have found that the latter has a greater impact on the result of GWR than the former [42]. This paper used a common kernel function, Gauss kernel function [43], as the type of kernel function. Its function form is

$$w\_{ij} = \exp\left(-\left(d\_{ij}/b\right)^2\right) \tag{2}$$

where *wij* is the weight for unit *j* in the neighborhood of unit *i*, and *dij* is the distance between the center point of the unit *i* and *j* as the measurement of spatial proximity degree. *b* is the bandwidth of the Guass kernel function. Many approaches are available for determining the bandwidth, including cross-validation (CV) method [42], Akaike information criterion (AIC) method [44]. Compared with CV method, AIC method is easier to avoid over-fitting problems, and the selected optimal model is often more effective. Therefore, the AIC method was adopted, and the bandwidth corresponding to the weight function with the minimum AIC value is the optimal bandwidth. Equation (3) is used for the calculation of AIC,

$$AIC = 2n\ln(\vartheta) + n\ln(2\pi) + n\left[\frac{n + tr(S)}{n - 2 - tr(S)}\right] \tag{3}$$

where *tr*(*S*) is the function of *b*, and σˆ is the maximum likelihood estimation for GWR model.

#### *3.6. Validation Processes*

The GEZM validation process is important; without this, the study lacks scientific credibility. In this research, we first compare the prediction results with actual units free of or containing landslides [45]. Sample points with 30% of the number of units, accounting for totally 291 landslide points, are randomly selected. The values of the predicted result belonging to the low-susceptibility and very-low-susceptibility area are regarded as the stable points and the others are as the unstable points. Then, calculate the missing rate, which is the ratio of actual unstable slopes classified as stable slopes.

Although the missing rate is one indicator that evaluates the model results, the results are subject to the threshold value, which limits the accuracy of results [46]. Therefore, we introduced ROC curve [47] and success rate curve [34] for further evaluation. The two methods are independent of the specific decision threshold and are able to further verify the accuracy of the result. ROC curve method is a measure of the ability to discriminate landslides from non-landslide locations. Success rate curve, different from the ROC curve, only considers the prediction of the landslide samples. The areas under the curve (AUC) for ROC curve and success rate curve (0.5 to 1.0) are used to assess the accuracy of the models.

#### **4. Results**

#### *4.1. Landslide Susceptibility Map Using the Proposed Model*

Based on the hydrologic analysis tool in ArcGIS environment, we segment the topography of the study area into 55,899 slope units in total. The schematic diagram of slope units is shown in Figure 3.

**Figure 3.** Schematic diagram of spatial proximity expression in slope units.

Firstly, in order to test the global multicollinearity, the variance inflation factor (VIF) value of each influencing factor is calculated (results shown in Table 2). A greater VIF value indicates a greater multicollinearity of the influencing factors. A value greater than 10 denotes that multicolinearity problems may exist. In this case, elevation factor for instant, the maximum VIF value is 1.88, which is apparently lower than 10. Table 2 indicates that all of the influencing factors used in the proposed S-GWR model pass the global multicollinearity test.


**Table 2.** Test results of variance inflation factor method of S-GWR.

Then, the local multicollinearity problem is subsequently processed with PCA method. According to the principles of selecting principal components, four principal components are selected in our study, and the components explained 81.10% of the total variances. After the factor rotation, four new principal components were obtained. As is shown in Table 3, the accumulative loadings of four principal components decrease systematically. The first component can be interpreted as geological factor, as it has high positive loadings of distance-to-fault and lithology. The second component has high positive loadings of slope and slope height, indicating that this component represents slope shape factor. In the same way, the third and the fourth component represent hydrographic factor and aspect factor, respectively.


**Table 3.** Rotated component loadings of the indicators on selected principal components.

The S-GWR model is constructed using the Gaussian kernel function and the AIC method. The neighbors of the model estimated by the AIC method is 1000. According to the prediction results, landslide susceptibility of Qingchuan County is mapped. The landslide susceptibility map is classified into five classes based on the natural breaks method [48], including very low, low, moderate, high, and very-high (shown in Figure 4). The landslide susceptibility map shows that the high-susceptibility area and very-high-susceptibility area are highly consistent with the actual landslide points, and very-low- to very-high-susceptibility classes occupy 46.37%, 25.67%, 17.62%, 8.69%, and 1.65% of study area.

**Figure 4.** Landslide susceptibility mapping (LSM) by the proposed S-GWR model.

#### *4.2. Validation of the S-GWR Models*

#### 4.2.1. Validation of the Spatial Non-Stationary

In order to validate the spatial non-stationary of the S-GWR model, we explore the spatial distribution of the regression coefficients between four principal components and landslide susceptibility. Results are shown in Figure 5. The coefficients estimated by S-GWR vary greatly with districts in space. This result implies a noticeable spatial non-stationarity of the relationship between four principal components and landslide susceptibility. The regression coefficients of each component have positive and negative values, which indicates that the relevance between each influencing factor and landslide susceptibility has varying direction and strength in space.

Spatial distribution of the regression coefficient of geological factors is approximately from northeast to southwest, similar to the distribution of actual lithology features and fault zones. The coefficient values of geological factors tend to decrease with the increase of the distance to fault zone and the large absolute values are mainly distributed in the metasandstones, phyllite, and sandy slates area. The distribution area with large absolute values of the coefficients of slope shape factor consist of the actual landslide points. In most of the eastern part of the study area, slope shape factor is positively correlated with landslide susceptibility, but the opposite effect exists in some central regions. The majority of the coefficient of hydrographic factor is negative, indicating that there is a significant negative correlation between the distance-to-stream and landslide susceptibility in most regions. This correlation is quite obvious in the central area, where most of the landslides occurred, such as Chaoba township, Courtyard Hui township, and the northeast area near the Bailong river. The aspect coefficient distribution map shows that the absolute value in the eastern region is larger than that in the western region. Compared with the first three principal components, the difference between the maximum and minimum values of aspect factor coefficients is the smallest.

**Figure 5.** Spatial distribution of the principal component coefficient function estimate. (**a**) Regression coefficient of geological factor (PC1); (**b**) Regression coefficient of slope shape factor (PC2); (**c**) Regression coefficient of hydrographic factor (PC3); (**d**) Regression coefficient of aspect factor (PC4)).

#### 4.2.2. Validation of the Accuracy

In order to validate the accuracy and effectiveness of the proposed S-GWR model, we compare the prediction results with actual units free of or containing landslides. Table 4 shows the prediction result of the proposed S-GWR model. Among the 291 actual landslide points, 36 points are classified as stable points with a missing rate of 12.37%. ROC curve and success rate curve are used for further evaluation. Figure 6 shows that the AUC value of ROC curve is 0.859, and the AUC value of success rate curve is 0.850.

**Table 4.** Sample evaluation contingency table of the proposed S-GWR model.

**Figure 6.** Validation of S-GWR model. (**a**) ROC curve; (**b**) success rate curve.

#### *4.3. Comparison with Other Models*

#### 4.3.1. Comparison with Grid-Unit-Based GWR (G-GWR)

In order to verify the applicability of slope units in expressing spatial proximity of the GWR, the commonly used grid unit is adopted as a comparison. For the validity of contrast and simplicity of calculation, grid size is selected as 200 × 200 m, approximately equal to the slope unit, generating 80,264 grid units.

The following procedure is the same as the process of the proposed S-GWR. The variance inflation factor (VIF) values of each influencing factor are calculated, with the maximum test result of 2.14 (as shown in Table 5). It indicates that the original G-GWR passes the global multicollinearity test. Based on the principles of selecting principal components, G-GWR also obtained four components, including geological factor, slope shape factor, hydrographic factor, and aspect factor. The establishment of kernel function of G-GWR is the same as S-GWR.

**Table 5.** Test results of variance inflation factor method of G-GWR.


The established G-GWR model is then applied to the landslide susceptibility assessment in Qingchuan County. The sample evaluation contingency table of the G-GWR model is shown in Table 6. It shows that 61 points of the 291 actual landslide points are classified as stable points, with a missing rate of 20.96%, which is higher than that of the S-GWR model with 12.37%. The ROC curves and success rate curves of the two models are shown in Figure 7. The AUC values of ROC curve and success rate curve of G-GWR are 0.837 and 0.827 respectively, both lower than those of S-GWR. The comparative verification of the S-GWR model and the G-GWR model shows that the S-GWR model has a better performance in classification precision, thus can better identify susceptible areas to landslide occurrence. This comparison indicates that slope unit is more suitable for GWR model in the expression of spatial relationship than traditional grid unit, in particular for the purpose of landslide susceptibility assessment.

**Table 6.** Sample evaluation contingency table of G-GWR model.

**Figure 7.** Validation of S-GWR and G-GWR models. (**a**) ROC curve; (**b**) success rate curve.

4.3.2. Comparison with OLS, ANN, and I Models

In order to further verify the applicability and accuracy of the landslide susceptibility assessment using the GWR model, we choose the ordinary least squares (OLS) model [40], artificial neural network (ANN) model [32], and information (I) model [49] as compared models, which are commonly used for the susceptibility assessment. The slope units and grid units are adopted as the sampling units of the three models respectively.

The ROC curve and success rate curve of eight models are shown in Figures 8 and 9. It can be seen from the figures that for the same unit type, the AUC values of ROC curve and success rate curve of GWR models are both higher than those of ANN, I, and OLS models. It shows that the GWR model performs better, and further implies the importance of spatial non-stationarity in landslide data sets. As to the results by the same model via different unit types, the AUC values of ROC curve and success rate curve of the models with slope unit are all higher than those of the models with grid unit. It indicates that adopting slope unit as mapping unit can express the landslide characteristics better than the traditional grid unit, highlighting that slope unit is more suitable for landslide assessment.

**Figure 8.** ROC curves and area under the curves (AUC) for the susceptibility maps produced by GWR, ANN, I, OLS. (**a**) Slope unit; (**b**) grid unit.

**Figure 9.** Success rate curves and area under the curves (AUC) for the susceptibility maps produced by GWR, ANN, I, OLS. (**a**) Slope unit; (**b**) grid unit.

#### **5. Conclusions**

In this paper, a spatial proximity-based geographically weighted regression (S-GWR) model is proposed for assessing the landslide susceptibility. The presented model solves the issues of the spatial non-stationarity between landslide factors and its occurrence that usually neglected in previous landslide susceptibility assessment studies. In order to express spatial proximity properly, slope units are adopted. The multicollinearity between the data is eliminated through VIF method and principal component analysis.

The Qingchuan area in China after the 2008 Ms. 8.0 Wenchuan earthquake is selected as a case study to illustrate the effect and validity of the proposed S-GWR model. The result shows that the four influencing factors (including geological factor, slope shape factor, hydrographic factor, and aspect factor) in the S-GWR model all have noticeable spatial non-stationary effects on the landslide susceptibility. By quantitatively testing, the missing rate of S-GWR is 12.37%, which is much lower than that of G-GWR of 20.96%, and the AUC values of ROC curve and success rate curve of S-GWR (0.859 and 0.850, respectively) are both higher than those of G-GWR (0.837 and 0.827, respectively). Besides, the AUC values of the ROC curve and success rate curve of the GWR models are higher than those of ANN, I, and OLS models, and the accuracy of each model using slope unit is higher than those of the grid unit with similar cell size.

Our study verifies the importance of considering the spatial non-stationary and the applicability of the GWR model in the landslide susceptibility assessment. It also suggests that slope unit can better express the spatial relationship between landslide data and make the evaluation results more accurate.

**Author Contributions:** Conceptualization and methodology, Y.L. and Z.H.; Writing—validation, X.L.; Review and editing, J.D.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the National Key R&D Program of China (grant no. 2018YFD1100401); the National Natural Science Foundation of China (grant no. 41702310); the Natural Science Foundation of Hunan (grant no. 2018JJ3638); and the Innovation Driven Program of Central South University (grant no. 2019CX011). These financial supports are gratefully acknowledged.

**Acknowledgments:** We would like to thank W.W. in Central South University, China, for his constructive suggestions. We also extend our gratitude to editor-in-chief and three reviewers for their insightful comments during the review stage while we submitted the manuscript.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-3788-7