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

: Landslide mapping and analysis are essential aspects of hazard and risk analysis. Landslides can block rivers and create landslide-dammed lakes, which pose a signiﬁcant 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 ﬂight 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 inﬂuence the DEM generation from Sentinel-1 imagery.


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][5][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][10][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][13][14]. Optical imagery from different sources such as satellites, or unmanned aerial 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.

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 Snaefells 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]. 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 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, that approximately 7 million m 3 of material was released and that a total of 10-20 million m 3 of material was displaced. The landslide deposition covers an area of approximately 1.5 km 2 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.

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 preand post-event datasets, using the Google Earth Engine (GEE) platform.

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. 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 B perp ) 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 B perp 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 B perp 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 B perp between 96 and 170 m for landslide mapping with Sentinel-1. For this study, we downloaded the available Sentinel-1 images with a B perp 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). Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 22 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 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 B perp , the more accurate the altitude measurement [57,60]: 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 B perp , and a higher B perp 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 B perp 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.

Sentinel-1 Image Pairs
Orbit 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.

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

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

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][74][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.

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: where y i refers to the ith interpolated elevation, yt i 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.

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].
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]. (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.   (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. 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].

DEMs from Sentinel-1 and Landslide Volume Estimation
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]. (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.   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.

DEMs from Sentinel-1 and Landslide Volume Estimation
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.
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.  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.   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.  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. 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 m 3 of material was released from the landslide source area, and the volume of the debris in the depositions area is about 10-20 million m 3 . According to these statistics, the track 16 DEM shows an unrealistic value of 109 million m 3 , whereas the landslide volumes calculated from the other two DEMs are closer to the estimated reference value.

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 B perp 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 B perp [76]. The Sentinel-1 mission was mainly designed for the retrieval of surface deformations using differential InSAR (DInSAR) which requires a low B perp between consecutive SAR images, and DEM generation was not the primary goal of the mission [58]. The recommended B perp 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 B perp , whereas the common B perp 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 B perp 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 B perp 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 B perp , i.e., track number 155 with a B perp 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 B perp (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.

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