Detailed Mapping of Lava and Ash Deposits at Indonesian Volcanoes by Means of VHR PlanetScope Change Detection

: Mapping of lava ﬂows in unvegetated areas of active volcanoes using optical satellite data is challenging due to spectral similarities of volcanic deposits and the surrounding background. Using very high-resolution PlanetScope data, this study introduces a novel object-oriented classiﬁcation approach for mapping lava ﬂows in both vegetated and unvegetated areas during several eruptive phases of three Indonesian volcanoes (Karangetang 2018/2019, Agung 2017, Krakatau 2018/2019). For this, change detection analysis based on PlanetScope imagery for mapping loss of vegetation due to volcanic activity (e.g., lava ﬂows) is combined with the analysis of changes in texture and brightness, with hydrological runoff modelling and with analysis of thermal anomalies derived from Sentinel-2 or Landsat-8. Qualitative comparison of the mapped lava ﬂows showed good agreement with multispectral false color time series (Sentinel-2 and Landsat-8). Reports of the Global Volcanism Program support the ﬁndings, indicating the developed lava mapping approach produces valuable results for monitoring volcanic hazards. Despite the lack of bands in infrared wavelengths, PlanetScope proves beneﬁcial for the assessment of risk and near-real-time monitoring of active volcanoes due to its high spatial (3 m) and temporal resolution (mapping of all subaerial volcanoes on a daily basis).


Introduction
Hazards triggered by volcanic activity have massive environmental, economic, and humanitarian impacts [1,2]. Lava flows can destroy entire settlements and infrastructure [3]. Ash and tephra deposits cause damage to critical infrastructure, buildings, and agricultural land [4] and pose a threat to human health [5]. Furthermore, pyroclastic flows [6,7] and lahars [8] are responsible for a large number of the fatalities due to volcanic hazards [1]. Secondary volcanic hazards, particularly volcano-induced tsunamis, also possess large destructive forces [9,10]. Due to the catastrophic impacts of volcanic hazards, monitoring these processes is an essential component for risk assessment and planning in disaster management [11][12][13]. For that, satellite-based data are increasingly being integrated into the monitoring process of volcanoes, as they allow large-scale views of areas that are difficult to reach and hazardous due to active eruptions. In addition to the frequently used infrared (IR) (e.g., [14][15][16]) and synthetic aperture radar (SAR) (e.g., [10,17]) data, optical satellite data are also useful for the monitoring of active volcanoes. Very high-resolution (VHR) optical systems have been used to map ash and tephra deposits [18], lava flows [19], pyroclastic flows, and lahars [20,21] at volcanoes. Lower-resolution optical imagery (e.g., Sentinel-2, Landsat-7/-8) (10-30 m) have also been successfully used to detect volcanic hazards (e.g., [22][23][24]). A commonly used method for mapping and quantifying volcanic hazard deposits using optical satellite data is the differencing of the Normalized Difference Vegetation Index (NDVI) [18,23]. The results of this methodology can provide information In this work, lava flow deposits were mapped during two different eruptive phases at Karangetang (Table 1). The first eruptive phase (Karangetang I) began in November 2018 and continued until March 2019. During this period, large lava effusions developed at the N crater, resulting in a lava flow towards NNW. A second lava flow also flowed WNW from the N crater in late November 2018 [31]. Effusive activity continued throughout December 2018 into January 2019 and was repeatedly accompanied by avalanches and pyroclastic flows. On 2 March 2019, the lava flow on the NNW flank crossed a road near the coast. Two days later, the lava flow reached the coastline in the NNW of the island and ended in the ocean. During the activity in February 2019, 190 people had to be evacuated. From 12 February 2019 onwards activity decreased and the eruption phase ended in March 2019 [31].  The second eruptive phase at Karangetang (Karangetang II) occurred between July 2019 and January 2020 (Table 1) and was characterized by a variety of incandescent block debris avalanches, lava flows, frequent gas and steam emissions, and ash clouds [32,33]. In addition to block debris avalanches with incandescent material, a WNW lava flow developed from the main crater in late July 2019. Lava effusion increased in the following weeks and continued to develop downslope in August and September 2019. Lava effusion was also observed at the N crater during this period. Following a large Strombolian eruption in late November 2019, steady lava effusion intensified, creating a large-scale lava field on the W flank of Karangetang [32]. During December 2019, lava effusion gradually Remote Sens. 2022, 14,1168 4 of 31 decreased and incandescent block debris avalanches became dominant on the W flank. This activity continued until mid-January 2020. Thereafter, volcanic activity abruptly stopped, and only isolated hotspots were sighted in both craters [33].
Mount Agung (8.343 • S; 115.508 • E) ( Figure 1) is an active stratovolcano (2997 m a. s. l.) located in the east of the island of Bali on the Sunda Island Arc [34]. The volcano is composed mainly of basaltic-andesitic rocks and has a partly explosive Plinian or Sub-Plinian eruptive form. This often results in pyroclastic flows and lahars [34]. Due to the high population density of the island of Bali, Agung is a particularly dangerous volcano [35].
In mid-August 2017, a new eruptive phase began (Agung I) after a prolonged period of inactivity [34,36]. This intensified in November 2017 (Table 1), with massive ash eruptions occurring on 21 November 2017 and the subsequent days. In the process, ash columns rose to heights ranging from 4.8 to 9.3 km a. s. l. between 25 and 26 November 2017 [37]. During this period, the ash cloud drifted mainly W and SW [36,37]. The eruptions resulted in the closure of the international airport in Bali (60 km SW of Agung) for several days between 26 and 29 November 2017. In addition, approximately 25,000 people were evacuated in the immediate vicinity of Agung. The largest eruptions continued until 30 November 2017. After that, only smaller ash and gas eruptions were observed [36].
Krakatau (6.102 • S; 105.423 • E) ( Figure 1) is a volcanic complex in the Sunda Strait that consists of a caldera with a stratovolcano inside. The volcanic complex is composed primarily of andesites or basaltic andesites [38]. It is characterized by frequent alternation between Strombolian or Vulcanian eruptions [10]. However, violent phreatomagmatic and Surtseyan eruptions are also common [9,10]. Mainly lava flows and pyroclastic flows occur during eruptions of Krakatau [10]. A caldera collapse during the catastrophic 1883 eruption destroyed the majority of the pre-1883 Krakatau Island. Volcanism continued after the 1883 events, eventually producing Anak Krakatau, "the son of Krakatau" [39]. The islands located in the immediate surroundings of Anak Krakatau are all uninhabited [38], which is why there is no direct threat to humans from the volcanic hazards of Anak Krakatau. However, a major secondary threat to the densely populated coastal regions of Java and Sumatra are volcano-induced tsunamis [9,40], which led to catastrophic consequences following major eruptions of Krakatau in 1883 [39,40] and in 2018 [9,10].
In this work, two volcanic events at Krakatau were analyzed ( Table 1). The first eruptive phase (Krakatau I) began in late June 2018, with lava flows developing on the S and SE flanks of Anak Krakatau accompanied with Strombolian activity. In addition, movement occurred on the SW flank due to the volcanic activity [10]. Lava effusion and Strombolian activity continued into September 2018, whereby an increase in eruptive activity was observed [41]. As a result, steady lava effusion flowed into the sea off the S coast of Anak Krakatau, creating 12.09 ha of new land [19]. In October and November 2018, lava flows also occurred more frequently in the SW of the island. Throughout the eruption period, volcanic activity was frequently accompanied by explosions with large ash plumes.
On 22 December 2018, the described eruptive phase (Krakatau I) subsequently resulted in a powerful explosion (Krakatau II), which brought significant morphological changes [42]. The explosion was accompanied by seismic activity that led to the collapse of the SW flank of Anak Krakatau [10]. Hence, a devastating tsunami developed which reached the coastal regions of Java and Sumatra and the surrounding islands of the Krakatau volcanic complex [9]. In total, 437 people were killed and over 31,000 were injured [42]. The eruption also resulted in large-scale tephra deposition on the surrounding islands [10,19].

Datasets
The basis of the developed change analyses for the detection of volcanic hazards is PlanetScope data. PlanetScope is a satellite constellation of over 150 CubeSats (Doves) (10 cm × 10 cm × 30 cm) capable of recording the entire Earth's landmass (approx. 1.49 × 108 km 2 ) on a daily basis [43,44]. The constellation has been continuously recording data since June 2016 and has been extended by additional Dove generations in subsequent years. In total, three different generations exist, all of which have slightly different sensor Remote Sens. 2022, 14, 1168 5 of 31 systems (1st generation: PS2, 2nd generation: PS2.SD, 3rd generation: PSB.SD). Depending on the generation, PlanetScope satellites are equipped with three (blue, wavelength λ = 455-515 nm; green λ = 500-590 nm; red λ = 590-670 nm) or four bands (additional near infrared (NIR) λ = 780-860 nm) [44]. PlanetScope has a geometric resolution of 3.7 m, which is resampled to 3 m within the provided data products [44]. Depending on the different Dove generations, the swath width thereby varies between 25-30 km × 8-10 km [18,44]. Furthermore, almost all satellites fly at an altitude of 475 km in a near-polar orbit with inclinations of 8 • and 98 • [43]. In addition to daily coverage of the entire Earth's landmass, overlap areas occur with successive Doves, which are recorded multiple times with a time offset of a few minutes and can thus be used to observe short-term, dynamic processes at the Earth's surface [45].
Complementary to the PlanetScope imagery, Sentinel-2 Multispectral Instrument (MSI) and Landsat-8 Operational Land Imager (OLI) data were included in the change analysis. First, these datasets were used to create multispectral false-color representations during eruption periods at the analyzed volcanoes. These were used for visual analysis and validation of the lava area derived from the PlanetScope. Second, thermal anomalies (so-called hotspot), detected by the Normalized Hotspot Indices (NHI) [46], were derived from the Sentinel-2 and Landsat-8 imagery during an eruption period and used as an input dataset for object-oriented lava flow mapping (Section 4.3). The Sentinel-2 satellites (2A and 2B) fly in a near-polar sun-synchronous orbit at an altitude of 786 km and have a temporal resolution of 5 days. The MSI sensor on board of the Sentinel-2 satellites has thirteen spectral bands in the wavelength range between visible light (VIS) (λ = 443 nm) and shortwave infrared (SWIR) (λ = 2190 nm) in three different geometric resolutions (10 m, 20 m, 60 m) [47].
Landsat-8 flies in a near-polar sun-synchronous orbit at an altitude of 705 km and has a repetition rate of 16 days. The OLI sensor has nine spectral bands with a geometric resolution of 15 m (panchromatic) or 30 m (VIS-SWIR) in the wavelength range from VIS (λ = 443 nm) to SWIR (λ = 2290 nm) [48]. In addition, Landsat-8 has the Thermal Infrared Sensor (TIRS), which acquires two bands in the thermal wavelength range (λ = 10.6-12.51 µm) with a geometric resolution of 100 m [48].
For the generation of false-color representations for monitoring thermal anomalies (Section 6.1), Sentinel-2 and Landsat-8 data were obtained via Google Earth Engine (GEE) for the Karangetang I, Karangetang II, and Krakatau I eruption events. Furthermore, two additional Sentinel-2 L2A scenes were downloaded for eruption phase Karangetang II. These were used to test the application and functionality of the developed object-oriented change detection method (Section 4.3) with data having a lower geometric resolution than PlanetScope (Section 6.3). For this purpose, a pre-(9 May 2019) and post-(4 December 2020) eruption scene were selected.
Due to the good geometric resolution and the availability of SWIR bands, Sentinel-2 and Landsat-8 are well suited for monitoring thermal activities at active volcanoes [15,49].
Using the Normalized Hotspot Indices (NHI) [46], volcanic hotspots can be mapped in detail based on the well-resolved IR bands of Sentinel-2 (20 m) and Landsat 8 (30 m). Based on the particularly strong emission of hot objects at the Earth's surface in the SWIR wavelength range, as well as the shift of the maximum emission towards shorter SWIR wavelengths as the temperature of the objects increases, volcanic hotspots can be identified and classified according to their intensity. The detailed methodology of the NHI is described in Marchese et al. [46]. The NHI hotspots for all active Holocene volcanoes, listed by the Global Volcanism Program (GVP), can be accessed online using the GEE tool (https://nicogenzano.users.earthengine.app/view/nhi-tool) [50]. For the analyzed events (Karangetang I, Karangetang II, Krakatau I), NHI hotspot pixels were derived and used as input data for the object-oriented change detection methodology for lava flow classification (Section 4.3).

Pre-Processing
To ensure the comparability of the pre-and post-scenes, PlanetScope data were downloaded as an "Ortho Scene" product [44] and projected onto a local UTM zone of the volcanoes. Individual PlanetScope scenes (of the same acquisition date) were then mosaiced for the respective study area. In the process, the overlap areas of the scenes were averaged. Subsequently, clouds or vapor and smoke plumes were manually masked out in all scenes. To minimize geometric redundancies between the pre-and post-scene, the PlanetScope scenes were co-registered. Histogram matching was then applied to the Plan-etScope image pairs. Finally, the NDVI (Equation (1)) [55] was calculated for each scene and appended to the respective dataset.
For monitoring of thermal anomalies using Sentinel-2 and Landsat-8 the RGB band combination SWIR2, SWIR1, NIR was used. Thus, for Sentinel-2 MSI data, bands 12, 11, and 8a were used. For Landsat-8 OLI, bands 7, 6, and 5 were applied. The two Sentinel-2 scenes (pre-and post-scenes Karangetang II) were reduced to the bands with a pixel size of 10 m (B02 (blue), B03 (green), B04 (red), B08 (NIR)). Subsequently, clouds and vapor or smoke plumes were manually masked out and image-to-image co-registration was performed. In the last step, the NDVI was calculated for both scenes and added to the dataset.
The NHI hotspot pixels were extracted with a corresponding spatial distance filter, so that hotspots were detected exclusively in the close surroundings of a volcano. The NHI vector data include the total hotspot area per observation from Sentinel-2 and Landsat-8, respectively, for which the three intensity classes (low-mid, high, extreme) [50] were combined. The output NHI pixels were aggregated over the entire eruption period so that a final dataset was available for each eruptive phase. This represents the total area of all hotspot pixels during the eruption and suggests that these areas have been affected by volcanic activity (e.g., lava flows). The total hotspot area aggregated over time was then implemented in the object-oriented change analysis (Section 4.3).
The DEMs (TanDEM-X 12 m, Copernicus GLO-30 m, SRTM 30 m, ALOS PALSAR RTC DEM 12.5 m) were projected onto the UTM zones of Karangetang, Agung, and Krakatau, respectively, using a bilinear resampling method. The DEM scenes were then reduced to the extents of the volcanoes. These data were then used to generate the respective volcanoes' drainage system.
The scheme of O'Callaghan and Mark (1984) was used to derive the runoff networks from DEM data. First, the sinks present in the DEMs were filled to avoid the derivation of discontinuous flows. The subsequently computed flow direction matrix includes a directional coding that channels the runoff flow in the direction of the adjacent pixel with the steepest slope [56,57]. From this, the flow accumulation matrix can be derived. This layer represents the cumulative weighting of all cells located upstream and drain out into a given pixel [57,58]. Flow accumulation was calculated using the D8 algorithm [56], with no weighting of individual cells [59]. Accordingly, cells with high flow accumulation values are areas where runoff is concentrated. These areas can be considered the main outlets of the drainage network [57,58]. Thus, by setting a threshold, the main flows can be identified and extracted. In this process, finding an appropriate threshold can be done manually based on geographic and site-specific factors [56]. Due to different environments at the volcanoes, as well as the different geometric resolutions of the DEMs, the threshold values for extracting the main flows were identified manually. For the Karangetang I and II events the TanDEM-X DEM (12 m) (flow accumulation threshold 50) was used. The Krakatau I event was analyzed using the Copernicus DEM (30 m) (flow accumulation threshold 10). Finally, the extracted flow networks were converted to a binary raster, which was used as input for the object-oriented mapping of lava flows (Section 4.3).

Volcanic Ash Deposits Mapping
For mapping volcanic hazards in vegetated areas, a dNDVI analysis was performed. The NDVI describes the normalized ratio of reflectance in the NIR and red (Equation (1)) and provides information about the photosynthetic activity or health of plants [55], where ρ NIR and ρ RED represent the corresponding reflections in the NIR and red wavelength ranges.
In this work, dNDVI analyses were performed for the Agung I and Krakatau II eruption events to identify and quantify ash and tephra deposits. In addition, dNDVI layers were created for the Karangetang I, Karangetang II, and Krakatau I eruption phases, which were subsequently incorporated into the object-oriented methodology for lava flow mapping (Section 4.3).
For the identification of volcanic induced vegetation loss, the NDVI layers of the pre-and post-scenes (exact dates of imagery are reported in Sections 5.1 and 5.2) of the respective eruption phase were differentiated using Equation (2): The direction of change can then be inferred from the resulting dNDVI layer [60], whereby areas of negative dNDVI values represent volcanic deposits that cover or destroyed vegetation. A threshold classification of the dNDVI layer was then used to classify the intensity of change. In doing so, the events used to analyze ash and tephra deposits (Agung I and Krakatau II) were classified with |0.2| intervals. For the analysis of eruption events with the lava flows (Karangetang I, Karangetang II, Krakatau I), the dNDVI change areas were classified using a −0.3 threshold and added to the object-oriented change detection method as an input dataset.

Lava Flow Mapping
The separation of cooled lava flows or other volcanic deposits and the background in unvegetated areas is hardly possible due to spectral similarities in the VIS and NIR [22], which is why a dNDVI analysis (Section 4.2) does not yield any valuable results [18]. Especially at volcanoes with high eruptive activity, more overlapping lava effusions take place. Due to these dynamics, the overlapping lava deposits exhibit mineralogical and morphological similarities that lead to identical spectral signatures [22,24,26]. Accordingly, at the end of an eruptive phase, it is difficult to distinguish which individual lava flows drained where and what its age structures are [22,24]. Vegetation overprinting of cooled lava or lahar deposits can also complicate the detection with optical remote sensing data [22,24,61].
To address this issue, a change detection method was developed to map lava deposits in unvegetated areas. The methodology follows an object-based image analysis (OBIA) approach [62]. This allows for the inclusion of spatial, structural, and contextual information into the classification, leading to the exploitation of the full data potential, especially in the case of VHR satellite imagery [62,63].
Input data of the OBIA change detection are the PlanetScope scenes before and after an eruption event (pre-and post-scenes), the classified dNDVI layer (threshold <−0.3) (Section 4.2), the temporally aggregated NHI layer, and the extracted drainage system (Section 4.1). The schematic workflow of the OBIA is shown in Figure 2. approach [62]. This allows for the inclusion of spatial, structural, and contextual infor-mation into the classification, leading to the exploitation of the full data potential, especially in the case of VHR satellite imagery [62,63].
Input data of the OBIA change detection are the PlanetScope scenes before and after an eruption event (pre-and post-scenes), the classified dNDVI layer (threshold <−0.3) (Section 4.2), the temporally aggregated NHI layer, and the extracted drainage system (Section 4.1). The schematic workflow of the OBIA is shown in Figure 2. pre-and post-PlanetScope scene, temporally aggregated NHI pixels, dNDVI layer, and the drainage system.
Step 1: extraction of the drainages in the overprint area (NHI or dNDVI area); step 2: classification of segments in the unvegetated areas using the combination of distance to drainage and changes in texture (GLCM) and brightness. Starting from the "lava flow area in vegetation or NHI area" objects, a pixel-based region growing into the "potential lava flow area with change" change areas is used to model the lava flow to its point of origin.
In the first step ( Figure 2), the dNDVI and NHI surfaces are intersected with the drainage system. In this way, the drains can be identified which lie within the dNDVI or NHI areas and therefore could potentially have carried a lava flow because of a volcanic or thermal overprint. Then, starting from the identified drainage sub-areas in the dNDVI or NHI surfaces, the entire coherent drainage flow is detected using a region-growing algorithm. In this process, starting from the runoff flow segments within the dNDVI or NHI areas, the object grows exclusively into areas that are unvegetated in the post-scene (NDVIpost layers <0.4). This step ensures that only runoff networks within unvegetated Figure 2. Workflow of OBIA change detection based on NHI pixels and the dNDVI area. Input data: pre-and post-PlanetScope scene, temporally aggregated NHI pixels, dNDVI layer, and the drainage system. Step 1: extraction of the drainages in the overprint area (NHI or dNDVI area); step 2: classification of segments in the unvegetated areas using the combination of distance to drainage and changes in texture (GLCM) and brightness. Starting from the "lava flow area in vegetation or NHI area" objects, a pixel-based region growing into the "potential lava flow area with change" change areas is used to model the lava flow to its point of origin.
In the first step ( Figure 2), the dNDVI and NHI surfaces are intersected with the drainage system. In this way, the drains can be identified which lie within the dNDVI or NHI areas and therefore could potentially have carried a lava flow because of a volcanic or thermal overprint. Then, starting from the identified drainage sub-areas in the dNDVI or NHI surfaces, the entire coherent drainage flow is detected using a region-growing algorithm. In this process, starting from the runoff flow segments within the dNDVI or NHI areas, the object grows exclusively into areas that are unvegetated in the post-scene (NDVIpost layers <0.4). This step ensures that only runoff networks within unvegetated areas are detected in the post-scene, as these are the only areas where lava deposits may be present. If no vegetation change occurred, only the NHI pixels can be used as the initial area for intersection with the drainage system. Similarly, if no NHI pixels are present during the eruption period, only the dNDVI layer can also be used as the initial change area. The outflows identified in this step are then used as the basis for OBIA lava flow mapping.
In the second step ( Figure 2), the extracted drains and the VHR PlanetScope data are used to identify areas where lava could potentially have flowed. These are inferred from the distance to the drainage system and texture and brightness changes derived from the pre-and post-event PlanetScope imagery. This combination can be used to ensure that differentiation between lava flows and background can be achieved even in the unvegetated area. For this purpose, the unvegetated area in the post-scene is first identified using a threshold value (NDVIpost < 0.4). Then, a multi-resolution segmentation [63] is performed within this area to create the segments in the highest hierarchical object level 1 ( Figure 2). Lava flows are not bound to a specific shape and can take on irregular shapes, especially for smaller objects that do not cover an entire lava flow. At the same time, despite spectral similarities between lava flows and background, differentiations can be made between objects based on spectral properties [24]. For this reason, the used multiresolution segmentation placed a stronger emphasis on spectral properties rather than geometric properties.
After segmentation, the distance of the segment center of each object to the nearest point of the flow line is calculated. Objects that are within a distance of ≤15 PlanetScope pixels (45 m) are assigned to the class "potential lava flow area" and form the object Remote Sens. 2022, 14, 1168 9 of 31 level 2. In addition, all objects of the class "non-vegetation area" are intersected with the dNDVI and NHI layers to identify objects with vegetation changes or thermal overprinting (class "lava flow area in vegetation or NHI area") ( Figure 2). In the next step, texture and brightness changes in the objects of the potential lava area (class "potential lava flow area") are calculated. Since almost no spectral differentiation is possible between objects with lava deposits and the background [22,24,26], information about texture and brightness are included in the classification. Texture was calculated using a grey level co-occurrence matrix (GLCM) based on Haralick et al. [64] for all objects in the potential lava flow area class. In this analysis, the GLCM was calculated based on the individual objects using the red band for the pre-and post-PlanetScope scenes. Following the recommendation of Haralick et al. [64], the average of all four possible viewing diagonals (0 • , 45 • , 90 • , 135 • ) was calculated. The differences of the secondary statistics homogeneity, entropy and contrast were used. For this purpose, each feature was calculated in the pre-and post-scene and then subtracted from each other. In addition, two post-scene GLCM statistics, contrast and spatial standard deviation, were integrated into the classification. Another differentiator is the overall brightness of the objects. Here, a ratio of the overall brightness from the pre-and post-scene was formed to detect changes in the brightness of volcanic deposits. The mentioned texture and brightness parameters were used with the help of a respective threshold value for the classification of the objects of the class "potential lava flow area". The thresholds are strongly dependent on the local factors of the volcanic environment and were therefore defined manually for each individual volcano. The exact threshold values for the individual analyses are presented in Table 2. If an object exceeds or falls below one of the defined limits, it is assigned to the class "potential lava flow area with changes". These objects are located at a plausible distance from the drainage system and contain texture or brightness changes in the PlanetScope scenes. Table 2. Optimal threshold values of GLCM statistics and overall brightness ratio for the classification of lava flows of the events Karangetang I and II and Krakatau I. If a segment meets one of the defined thresholds, it is assigned to the class "potential lava flow area with changes" (see Figure 2). In the final process of the OBIA change detection, a pixel-based region-growing procedure is performed, which has as a starting point the objects of the class "lava flow area in vegetation or NHI area" (Figure 2). From these objects, which have a certain volcanic overprint due to vegetation changes or thermal activities, the algorithm grows pixel-based into all neighboring objects belonging to the class "potential lava flow area with changes". A pixel is included in the objects of the initial class if it fulfils one of the defined conditions of GLCM differences, GLCM post-scene, or overall brightness ratio ( Table 2). This process is iterated until none of the adjacent pixels of the object class "potential lava flow area with changes" fulfil any of the defined threshold conditions. The resulting class "lava flow area" now contains all objects that represent the area of the lava flow. Finally, all objects of the class "lava flow area" were merged to one object. All objects that are <15 pixels and completely enclosed by the lava area were included in this object. Using a majority filter, the edges and borders of the lava flow object were refined. The final lava area was modelled upstream starting from the volcanically overprinted area (dNDVI < −0.3 or NHI hotspot), through the distance to potential runoff streams and changes in the VHR satellite data (texture and brightness), to the potential origin of the lava in the non-vegetated area.

Volcanic Ash Deposits Mapping
Differential NDVI analyses (Section 4.2) were used to detect ash and tephra deposits after the Agung I and Krakatau II events.
For the dNDVI analysis of ash and tephra deposits after the Agung eruptions, suitable PlanetScope scenes before (18 September 2017) and after (7 December 2017) the event were selected. The post-scene had partially heavy cloud cover, so only the southern region of Agung could be analyzed. Figure 3 visualizes the NDVI changes after the eruptions around 21 November 2017 on Agung. The corresponding dNDVI layer was classified into several change intervals to visualize the magnitude of change [60]. Here, the lower the dNDVI values are, the stronger was the vegetation loss from pre-to post-scene. Areas with low dNDVI values have therefore experienced a large vegetation change due to volcanic activity [18]. The total area of ash visible in the imagery (dNDVI < −0.4) was quantified to be 25.78 km 2 . A striking feature is the spatial distribution of ash deposits on the southern flank of Agung, which become more intense with increasing proximity to the crater. In contrast, the immediate surrounding of Agung again shows higher dNDVI values.
Using a dNDVI analysis, the detection of ash and tephra deposition on Anak Krakatau and surrounding islands (Krakatau II) was possible. For this purpose, a PlanetScope scene before (14 March 2018) and after (1 July 2019) the eruption on 22 December 2018 was selected. In Figure 4, the effects of ash and tephra deposition as well as vegetation destruction due to the tsunami and explosive force of the eruption are visible at Anak Krakatau and the neighboring island of Krakatau Kecil. Overall, vegetation changes (dNDVI < −0.2) are present on a total area of 3.4 km 2 . Particularly affected by ash and tephra deposition was Krakatau Kecil, located to the east [42]. Figure 4 shows the decrease in NDVI signal over nearly the entire island area, especially on the west coast. On Anak Krakatau, a sharp decline in vegetation signals can be seen in the east of the island. Using a dNDVI analysis, the detection of ash and tephra deposition on Anak K tau and surrounding islands (Krakatau II) was possible. For this purpose, a PlanetS scene before (14 March 2018) and after (1 July 2019) the eruption on 22 December 2018 selected. In Figure 4, the effects of ash and tephra deposition as well as vegetation des tion due to the tsunami and explosive force of the eruption are visible at Anak Krak

Lava Flow Mapping
Object-oriented change detection methods (Section 4.3) were used to analyze the flows of the different eruptive phases Karangetang I, Karangetang II and Krakat Based on GVP reports [31] and NHI hotspot time series, the Karangetang I eruption e was temporally constrained and thus PlanetScope scenes from 19 July 2018 (pre-s and 4 May 2019 (post-scene) were selected for change analysis. Lava flow mapping performed using the OBIA change detection based on vegetation change derived from classified dNDVI layer (Section 4.3, Figure 2). TanDEM-X DEM (12 m) was used fo derivation of the drainage system and subsequent identification of potential lava f (Section 4.1). The derived lava area and the unaltered pre-scene are shown in Figu Using the OBIA change detection methodology, the total area of the lava flow o NNW flank of Karangetang was determined to be 0.49 km². The total length of the flow is 3425 m. At the widest point near the summit area, a width of 250 m was meas In the region where the road overflow occurred, a width of 140 m was measured. F 6 shows the derived classes used to identify and model the potential lava outflow ar

Lava Flow Mapping
Object-oriented change detection methods (Section 4.3) were used to analyze the lava flows of the different eruptive phases Karangetang I, Karangetang II and Krakatau I. Based on GVP reports [31] and NHI hotspot time series, the Karangetang I eruption event was temporally constrained and thus PlanetScope scenes from 19 July 2018 (pre-scene) and 4 May 2019 (post-scene) were selected for change analysis. Lava flow mapping was performed using the OBIA change detection based on vegetation change derived from the classified dNDVI layer (Section 4.3, Figure 2). TanDEM-X DEM (12 m) was used for the derivation of the drainage system and subsequent identification of potential lava flows (Section 4.1). The derived lava area and the unaltered pre-scene are shown in Figure 5. Using the OBIA change detection methodology, the total area of the lava flow on the NNW flank of Karangetang was determined to be 0.49 km 2 . The total length of the lava flow is 3425 m. At the widest point near the summit area, a width of 250 m was measured. In the region where the road overflow occurred, a width of 140 m was measured. Figure 6 shows the derived classes used to identify and model the potential lava outflow area.
The Karangetang II eruption event was temporally constrained using the GVP reports [32,33] and NHI hotspot time series. OBIA change analysis was performed using PlanetScope scenes from 4 May 2019 (pre-scene) and 3 May 2020 (post-scene). Because the described volcanic activity was mainly confined to the vegetation-free area on the upper W flank of the Karangetang, the temporally aggregated NHI pixels were also included in the OBIA method in addition to the dNDVI. The combination of the NHI and dNDVI areas formed the initial surface for the pixel-based region growing for mapping the lava surface (Section 4.3, Figure 2).      In the northern part of the W flank, lava also flowed into the vegetation area. The length of these lava deposits could be determined to be 1650 m. Figure A1 shows the NHI and dNDVI areas, the potential lava areas (based on distance from the drainage system) and the change areas (based on texture and brightness changes). described volcanic activity was mainly confined to the vegetation-free area on the upper W flank of the Karangetang, the temporally aggregated NHI pixels were also included in the OBIA method in addition to the dNDVI. The combination of the NHI and dNDVI areas formed the initial surface for the pixel-based region growing for mapping the lava surface (Section 4.3, Figure 2). Figure 7 shows the pre-and post-scene of the W flank together with the derived lava area. During the Karangetang II eruptive phase, 0.66 km² was covered by lava. The lava surface occupies almost the entire W flank of Karangetang and measures 880 m at its widest point. The resulting lava field is mainly concentrated in the non-vegetated area of the W flank. In the northern part of the W flank, lava also flowed into the vegetation area. The length of these lava deposits could be determined to be 1650 m. Figure A1 shows the NHI and dNDVI areas, the potential lava areas (based on distance from the drainage system) and the change areas (based on texture and brightness changes). Based on temporal constraints using GVP reports [41,42] and NHI hotspot time series, PlanetScope scenes on 14 March 2018 (pre) and 17 December 2018 (post) were selected for the Krakatau I eruption event. Since hardly any vegetated areas exist on Anak Krakatau, no corresponding dNDVI changes could be included in the methodology. Therefore, Based on temporal constraints using GVP reports [41,42] and NHI hotspot time series, PlanetScope scenes on 14 March 2018 (pre) and 17 December 2018 (post) were selected for the Krakatau I eruption event. Since hardly any vegetated areas exist on Anak Krakatau, no corresponding dNDVI changes could be included in the methodology. Therefore, only the temporally aggregated NHI pixels were used to classify and model lava flows using OBIA change detection. Figure 8 shows the pre-and post-scene and the OBIA-derived lava area at Anak Krakatau. A total area of 0.94 km 2 was affected by the lava effusion, which was almost exclusively concentrated on the S flank of Anak Krakatau. The widest part of the lava surface on the S flank has an E-W extension of 1000 m, while the maximum distance from the crater to the newly formed coastline in the SW of the island is 1120 m. Figure A2 includes the derived potential lava areas that lie within 45 m of the extracted drains, as well as the thermally overprinted NHI area and the alteration area based on texture and brightness changes. Figure 8 shows the pre-and post-scene and the OBIA-derived lava area at Anak Krakatau. A total area of 0.94 km² was affected by the lava effusion, which was almost exclusively concentrated on the S flank of Anak Krakatau. The widest part of the lava surface on the S flank has an E-W extension of 1000 m, while the maximum distance from the crater to the newly formed coastline in the SW of the island is 1120 m. Figure A2 includes the derived potential lava areas that lie within 45 m of the extracted drains, as well as the thermally overprinted NHI area and the alteration area based on texture and brightness changes.

Discussion of the Results
The methodology presented in Section 4.2 for the classification of ash and tephra deposits obtained good results when applied to the case studies of Agung I ( Figure 3) and Krakatau II (Figure 4).
Due to cloud cover, only a partial area of the S flank of Agung could be analyzed. Despite this, the PlanetScope scenes used were suitable for analyzing ash and tephra deposition due to the temporal proximity to the eruptions of Agung I. It is reported that the

Discussion of the Results
The methodology presented in Section 4.2 for the classification of ash and tephra deposits obtained good results when applied to the case studies of Agung I ( Figure 3) and Krakatau II (Figure 4).
Due to cloud cover, only a partial area of the S flank of Agung could be analyzed. Despite this, the PlanetScope scenes used were suitable for analyzing ash and tephra deposition due to the temporal proximity to the eruptions of Agung I. It is reported that the ash deposition happened over a wide area outside of the PlanetScope imagery extent [36,37]. Therefore, the area derived from the dNDVI analysis (25.78 km 2 ) (Section 5.1) massively underestimates the total ash deposit area. However, the spatial distribution of ash deposits in the near surroundings of Agung is comparable to reports from the GVP [36], which shows increased deposition on the SW flank. Here, ash deposition increases with proximity to the crater. The comparatively higher dNDVI values in the close vicinity of the crater can be explained with the unvegetated crater environment, which does not reveal any changes in the NDVI signal [18].
Within the analyses of the Krakatau II eruption, much of the negative dNDVI values on neighboring Krakatau Kecil Island can be explained in terms of ash and tephra deposits created by the powerful eruption of Anak Krakatau. However, it can be assumed that a large area of change can be attributed to the destruction of vegetation. The surrounding islands of Rakata, Sertung, and Kecil were all massively impacted by the tsunami [9,42]. Modelling of the tsunami shows a wave run-up of 10-23 m on the S, W, and NW coasts of Kecil [9], indicating the destruction of vegetation there and explaining the significantly lower dNDVI values in this area. Additionally, in the E of Anak Krakatau, the sharp decrease in the NDVI signal can be explained by the complete destruction of vegetation due to explosive eruptive activity [19].
The presented OBIA change detection method for mapping lava flows (Section 4.3) using VHR PlanetScope data also achieved promising results. For the Karangetang I eruption event, based on the dNDVI areas, the lava flow could be well mapped in the vegetation area. Despite of the difficulty to spectrally differentiate newly deposited lava from the unvegetated background areas close to the crater, good results could be obtained using the OBIA method. During the Karangetang I eruptive phase, according to the GVP [31], the lava flowed over a length of up to 3.5 km and was about 160 m wide at the road crossing location. These results are in good agreement with the derived lengths (3.425 km) and widths (140 m at the road crossing) of the OBIA classification (Section 5.2). W and E of the lava flow, levees can be recognized (Figure 5b), which could be formed partly due to the channeling of the lava flow and the natural topography of the outflow stream [65]. The color deviation from the dark color, typical for basaltic-andesitic lava, in the upper part of the lava flow can possibly be explained by a renewed overprinting of brighter deposits from a pyroclastic flow or ash [21]. The small lava flow on the WNW flank of Karangetang was recognized almost exclusively based on vegetation changes (dNDVI) (Figure 6b). Similarly, no connection between the smaller lava flow in the WNW and the large lava flow in the NNW could be generated. Due to the thresholding used in generating the main flow networks from the DEM, this area was not included in the drainage system modelling. In Figure 9, the classified lava area was compared with a timeseries of multispectral false color imagery (Sentinel-2). Using this methodology, a qualitative assessment about the validity of the classification can be made. The false-color representations clearly show the evolution of the active lava flow on the NNW flank of Karangetang and are consistent with reports from the GVP [31]. Here, the derived lava area includes nearly all of the detected thermal anomalies.
In contrast to reports [32], no lava deposits were identified at the N crater during the Karangetang II eruption phase (Section 5.2). Only a small area within the crater, which can be attributed to a hotspot detection in the NHI layer, was classified as lava. The derived lava area was assigned to the main crater by the region-growing algorithm (Figure 7b). The combination of the NHI and dNDVI layers covers mainly the more northerly lava deposits on the W flank ( Figure A1). Based on these surfaces, the texture and brightness changes also identified areas located farther south. Brightness was an important factor in detecting these areas of change. Visually, the lava overprinting can be seen in the significantly darker areas in the post-scene, which also make the entire flank area appear much more homogeneous (Figure 7). Furthermore, incandescent debris avalanches are frequently mentioned in the reports of the GVP [32,33]. These deposits could not be identified both visually and by OBIA change detection. Furthermore, the overlap of individual lava flows throughout the eruptive phase makes detection of individual lava flows nearly impossible [24].
The analysis of the multispectral false-color imagery (Sentinel-2 and Landsat-8) allows the identification of individual lava flows during the eruption ( Figure A3). There, the development of the lava flow in WNW direction at the beginning of the eruption period is visible. In the following months, sporadic thermal anomalies occurred in the W and WSW direction, but the main activity remained in the WNW. No lava effusion from N crater can be observed in the false-color imagery, confirming the derived lava area. Interrupted and fragmented thermal anomalies with a tongue-like structure are common ( Figure A3c,f,h,j,m-p,u). These anomalies may represent the block debris avalanches described in the GVP reports [32,33], which contain incandescent pyroclastic material and usually form at the front of lava flows [66]. On the other hand, fragmented thermal anomalies can also represent differ- entially cooled regions of a coherent lava flow [11,65,66]. Because the derived lava area includes nearly all of the thermal anomalies in the multispectral false-color time series, a high-quality detection of lava areas can be assumed. In contrast to reports [32], no lava deposits were identified at the N crater during the Karangetang II eruption phase (Section 5.2). Only a small area within the crater, which can be attributed to a hotspot detection in the NHI layer, was classified as lava. The derived lava area was assigned to the main crater by the region-growing algorithm ( Figure  7b). The combination of the NHI and dNDVI layers covers mainly the more northerly lava deposits on the W flank ( Figure A1). Based on these surfaces, the texture and brightness changes also identified areas located farther south. Brightness was an important factor in detecting these areas of change. Visually, the lava overprinting can be seen in the significantly darker areas in the post-scene, which also make the entire flank area appear much more homogeneous (Figure 7). Furthermore, incandescent debris avalanches are frequently mentioned in the reports of the GVP [32,33]. These deposits could not be identified both visually and by OBIA change detection. Furthermore, the overlap of individual lava flows throughout the eruptive phase makes detection of individual lava flows nearly impossible [24].
The analysis of the multispectral false-color imagery (Sentinel-2 and Landsat-8) allows the identification of individual lava flows during the eruption ( Figure A3). There, the development of the lava flow in WNW direction at the beginning of the eruption period is visible. In the following months, sporadic thermal anomalies occurred in the W and WSW direction, but the main activity remained in the WNW. No lava effusion from The analyzed eruption phase Krakatau I (Section 5.2) can be seen as a precursor of the catastrophic event on 22 December 2018 (Krakatau II). Walter et al. [10] were able to detect movement of the SW and S flanks within this phase using interferometric SAR (InSAR) time series. This ranged from 4 to 10 mm per month, with higher rates of movement correlating with increased volcanic activity, especially in September and October 2018 [10]. The location of the lava deposits in the S-SW region of the island (Figure 8b) is consistent with the observations of Walter et al. [10]. The results of the Krakatau I event prove that the OBIA method also works using only the NHI pixels as the input ( Figure A2b). In addition, the newly formed coastline in the SW and S of Anak Krakatau due to the lava effusion was delineated and confirms the increment of the island quantified by Ginting et al. [19]. In the SE of the island, there is a section of a recent lava deposit, which was formed shortly before the pre-scene (Figure 8). This was partially overprinted by the lava flows of the eruption phase. Here, the part bordering the sea, which was not overprinted during the 2018 activities, could be well demarcated from the lava deposits. Figure A4 can be used for the qualitative validation of the lava area. The lava area contains nearly all detected thermal anomalies. Only on 14 November 2018 ( Figure A4o) thermal anomalies occur outside the detected lava area. Due to the proximity to the crater and the circular distribution of the anomalies, ejected incandescent tephra (e.g., blocks, bombs) can be assumed here [2,12]. Additionally, visible is the increase in thermal activity in September 2018 [10,41]. As the month begins, increasingly large thermal anomalies are detected ( Figure A4e-k), which exert visible smearing effects of the thermal signal on their neighboring pixels due to their high emission ( Figure A4i,j) [27]. It was during this phase of eruptive activity that there was increased movement of the SW flank, which subsequently collapsed due to seismic activity [10], triggering the catastrophic tsunami [9].

Advantages and Limitations of the Presented Lava Flow Mapping Approach
Using the dNDVI analysis (Section 4.2), an effective methodology was applied to detect vegetation changes with PlanetScope data. Similar to the results of Aldeghi et al. [18], ash and tephra deposits could be detected very well, which in turn is an important parameter for further analysis of possible secondary natural hazards, e.g., secondary lahars [25]. However, this methodology is applicable only in vegetated areas. The dNDVI methodology is not suitable for detecting volcanic hazards in unvegetated areas as there are no temporal changes in the NDVI signal [18,23]. In addition, the detected vegetation change cannot be automatically assigned to a volcanic hazard. Natural vegetation changes triggered by fire, landslides, or seasonal changes also cause a reduction in the NDVI signal [20,23]. Explosive eruptions also often result in complete destruction of vegetation in the impact zone (Section 5.1, Figure 4). To identify these areas and not confuse them with ash or tephra deposition, a preliminary visual analysis with expert knowledge of the pre-and postscene is necessary [18]. Ash deposits can be very easily eroded by wind or precipitation, which is why the selection of pre-and post-scenes that are close in time to the eruption is particularly important [18]. After deriving the dNDVI layer, a threshold must be used to extract the area of change [60]. The definition of this threshold value proves to be difficult. The method used for manually finding a threshold based on local-specific factors and expert knowledge achieved good results in this regard [18,21]. However, this introduces a degree of subjectivity into the classification. Alternatively, automatable adaptive methods of threshold finding exist (e.g., Otsu method [67]), which can speed up and potentially optimize the process [60].
The developed OBIA change analysis (Section 4.3) also allows, in contrast to the dNDVI analysis, the classification of volcanic hazards in areas without vegetation. Due to the use of PlanetScope data in combination with modelling based on topographic data (DEM), the methodology offers decisive advantages over classical numerical modelling of volcanic hazards (e.g., MAGFLOW [68], LAHARZ [69]). These modelling efforts rely on strong parameterization of input variables and are used only to predict volcanic natural hazards [58,70]. By linking such modelling with information from satellite imagery, posteruption or near-real-time analyses can also be performed [14]. The developed OBIA methodology models the potential area of lava flows via a simple distance analysis to the drainage system. Using the PlanetScope data, spectral and textual information from the real world can then in turn be transferred to the model, which is not possible in the aforementioned numerical modelling. In addition, the full potential of the VHR PlanetScope data is exploited via the object-oriented approach. Thus, on the one hand, the reduction in spectral variability in the scene can contribute to an improved classification of lava flows [21]. On the other hand, contextual information about individual image objects can be obtained [63]. Regarding the scarcely existing spectral differences between lava flows and background [22,24], this is crucial for the function of the methodology. In addition, the analysis of texture and brightness changes leads to homogeneous and continuous lava flow segments within the object-based approach. This would not be possible in a pixel-based classification [24].
However, there are some limitations and potential sources of error in the OBIA methodology. For the methodology to function, alteration or overprinting of pixels due to dNDVI changes or NHI hotspots is essential. If this is not the case, potential outflows, via intersection with these areas of change, cannot be identified and lava flow modelling fails. In addition, there is always the possibility that areas of volcanic overprinting were not detected by the dNDVI analysis or NHI detection.
Another source of error for the modelling of the lava flows is the creation of a runoff model based on the DEMs (Section 4.1). Here, a high degree of generalization is assumed, which has an influence on the result. The resolution and quality of the DEM is elementary, because only in this way morphological structures can be recognized. In Figure 10, the OBIA methodology for the event Karangetang I was tested using different DEMs. All DEMs achieved satisfactory results, with the TanDEM-X (12 m) drainage model representing the lava flow surface most accurately. The comparison confirms that for modelling runoff, coarser resolution does not necessarily lead to worse results [70]. Despite this, higher resolution DEMs achieve smaller-scale runoff models, which in turn allow for a more detailed modelling of lava flows [70]. The DEMs used for the methodology were all generated before the analyzed eruptive phases. Thus, they may not represent the current prevailing morphology of the volcanoes. In particular, explosive eruptions can result in drastic morphological changes [10,71] that are disregarded in this modelling. For applying this classification approach after an eruption with severe morphological changes, an updated DEM must be available. One possible solution for this is to derive DEMs from (tri-) stereoscopic images of VHR optical satellites. These not only provide an actual image of the relief of a volcano after an eruption, but also have high geometric resolution [43]. However, the analyzed events mainly resulted in the discharge of lava with no major changes in the volcanoes morphology (except the lava deposits), which is why the used DEMs provide sufficient quality for modelling the lava flow path.
Another way to optimize drainage network modelling, is to calculate flow accumulation with a multi-flow algorithm. In contrast to single-flow algorithms, multi-flow algorithms achieve better results when modelling diverging drainage networks and parallel flows [72]. The single-flow D8 algorithm used in this study [56] (Section 4.1), models runoff into one of the adjacent cells based on the slope gradient. In contrast, multi-flow algorithms, e.g., FD8 [73], DINF [74], MFD-md [72], assume that runoff can flow into more than one adjacent cell, which gives a more realistic representation of the flows [72].
A strong generalization, and thus potential source of error, arises when defining a threshold for distance analysis of objects to the drainage network. If a drain is identified as a potential lava flow, all objects that lie within a defined distance are assigned to the class "potential lava flow area" (Section 4.3, Figure 2). Thereby a spatial generalization takes place, which does not correctly reflect the morphological structures of a complex runoff system. A possible optimization to derive the threshold value in the classification of potential lava areas, can be achieved by using the HAND (height above the nearest drainage) layer [75]. The HAND layer is calculated using a topographic algorithm that identifies potential floodplains based on the relative vertical distance to the drainage network [75].
Using the NHI (Section 3), thermal anomalies during the eruptive phase can be detected and used as an initial area for modelling a lava flow. However, very weak thermal anomalies can often be difficult to detect [50]. If a thermal anomaly is detected by the NHI, the entire pixel is classified as a hotspot. There is no sub-pixel analysis to identify the hot object within a pixel, so the total hotspot area is overestimated [50]. In addition, cloud cover and the temporal resolution of Sentinel-2 and Landsat-8 are limiting factors for hotspot detection [46]. Furthermore, no atmospheric correction is applied to the data used in the NHI tool, which can sometimes lead to noise in the SWIR signal [50]. Figure 11 shows the results of the OBIA methodology for the Karangetang II event in a version with and without the inclusion of NHI hotspots. This comparison shows that for eruptions that produce volcanic deposits in unvegetated areas, significantly more accurate results can be obtained by incorporating the NHI pixels. Due to minor changes within the vegetation, the dNDVI area only intersects a small portion of the drainage network. Accordingly, based on the distance analysis, the "potential lava flow area" class ( Figure 2) does not represent the entire lava flow area. Thus, the NHI hotspot area contributes significantly to the improvement of the results, especially in the vegetation-free space. Another way to optimize drainage network modelling, is to calculate flow accumulation with a multi-flow algorithm. In contrast to single-flow algorithms, multi-flow algorithms achieve better results when modelling diverging drainage networks and parallel flows [72]. The single-flow D8 algorithm used in this study [56] (Section 4.1), models run- and without the inclusion of NHI hotspots. This comparison shows that for eruptions that produce volcanic deposits in unvegetated areas, significantly more accurate results can be obtained by incorporating the NHI pixels. Due to minor changes within the vegetation, the dNDVI area only intersects a small portion of the drainage network. Accordingly, based on the distance analysis, the "potential lava flow area" class ( Figure 2) does not represent the entire lava flow area. Thus, the NHI hotspot area contributes significantly to the improvement of the results, especially in the vegetation-free space. Figure 11. Results of OBIA change detection methodology with and without NHI hotspots. The lava area derived without the incorporation of NHI pixels (in yellow) is significantly smaller and concentrated on the northern part of the W flank of the Karangetang. This is because in step 1 of the Figure 11. Results of OBIA change detection methodology with and without NHI hotspots. The lava area derived without the incorporation of NHI pixels (in yellow) is significantly smaller and concentrated on the northern part of the W flank of the Karangetang. This is because in step 1 of the OBIA methodology, the areas of change based on the dNDVI were intersected with the drainage network. Due to the only minor changes in vegetation, which were limited to the WNW outcropping lava tongue, only the area to the N could be detected as a potential lava area. Subsequent analysis of texture and brightness changes was limited to this area. The lava area with NHI hotspots (in red) is located further south on the W flank. This is due to the detection of hotspots in this area by the NHI, which were subsequently intersected with the drainage network.
Based on the texture and brightness changes, objects that have undergone a change, but are located outside the volcanically overprinted dNDVI or NHI areas can be identified. However, for the used classification variables (Table 2), the threshold values for class separation have to be manually adjusted for each application. A common problem with this is that the GLCM changes between pre-and post-scene are very small and do not robustly represent volcanic deposits in different study areas. Therefore, the optical derivation and subsequent definition of the respective thresholds proves to be difficult. By adding SAR amplitude information, the object-oriented classification of the alteration surfaces could be optimized. The deposits of volcanic hazards (e.g., lava flows, pyroclastic flows) change the surface roughness, which in turn leads to a change of the backscatter intensity of the radar signal [27,76]. Changes at the surface due to volcanic activity also cause decorrelation of the interferometric coherence when comparing coherence data derived from pre-event and from co-event SAR image pairs [77]. In addition, the use of SAR data is highly recommended, especially in tropical regions which have a high probability of a cloudy sky [17].
A quantitative validation of the presented results is not an easy task. Reference data are scarce due to costly and sometimes hazardous surveys [20,21,24]. Accordingly, a classical accuracy assessment cannot be performed. One possibility is validation via visual derivation of the deposits from VHR satellite or aerial imagery. Usually, the deposits of volcanic hazards are very complex and difficult to interpret visually [21,25], especially in unvegetated areas where there is little contrast between deposits and background. Accordingly, expert knowledge has to be available for this type of validation. Despite all this, subjective assessments may occur [21,24,66]. In Figure 12, the lava flow of the Karangetang I eruption phase was visually interpreted and manually digitized. This was compared with the semi-automatically generated results of the OBIA method. The agreement between the manually and semi-automatically derived areas is 85.17%, whereby the OBIA-derived lava area showed an overestimation by 14.83% compared to the manually derived reference. Another way to qualitatively validate the lava areas is to produce multispectral false color time series (Sentinel-2 and Landsat-8) as performed in Section 6.1 (Figures 9, A3 and A4). Due to the sensitivity of the SWIR bands to thermal anomalies, greater contrast can thus be produced [66] and lava flows can be tracked during an eruption period. Other validation opportunities include comparison with current geologic maps, which are scarce at remote volcanoes [21,71], or comparison with numerical modelling results [50,58].
The OBIA methodology can be applied to other tropical volcanoes. It is important that the study area around the volcano has vegetation, as this allows potential dNDVI changes to be measured [23]. Application of the methodology in areas with little vegetation depends heavily on the availability of NHI pixels and is unlikely to provide satisfactory results ( Figure 11). In addition, the OBIA change detection methodology can theoretically be applied to map other volcanic hazards that exhibit flow behavior, e.g., pyroclastic flows and lahars, with the latter requiring vegetated volcano flanks, as lahars cause no thermal anomaly that could be detected by the NHI algorithm.

How Volcano Monitoring Benefits from PlanetScope Imagery
In this work, the potentials of PlanetScope data for monitoring volcanic hazards were investigated. Based on the results presented in Section 5, it can be confirmed that PlanetScope is a useful data source in volcano monitoring when using appropriate methods, especially the high geometric resolution (3 m), which brings great advantages over other sensor systems. With this, it is possible to identify small-scale volcanic deposits and activities. Harris et al. [23] encountered a significant underestimation of marginal areas due to partially covered pixels when mapping lahar deposits using Landsat data (30 m). This problem is resolved with the use of higher spatial resolution PlanetScope data. When mapping ash and tephra deposits, satisfactory results can also be obtained even with lower-resolution sensors (e.g., Sentinel-2 MSI (10 m)) [18]. Figure 13 compares the use of the OBIA change detection methodology using Sentinel-2 with the PlanetScope results. For the mapping of the lava area after the Karangetang II eruption phase, only the 10 m bands (B02 (blue), B03 (green), B04 (red), B08 (NIR)) of Sentinel-2 MSI were used to ensure comparability with PlanetScope. The methodology provides satisfactory results even with lower resolution Sentinel-2 (10 m) data. The agreement of the Sentinel-2 lava area with the PlanetScope lava area is 79.65%. At the same time, 13.4% of the Sentinel-2 lava area lies outside the boundaries of the PlanetScope lava area. This suggests that the methodology is robust to the geometric resolution of the input data. However, with lower geometric resolution, the detection of the boundary areas deteriorates, which was also observed in Harris et al. [23]. In addition, the advantages of the object-based approach become blurred [63]. Accordingly, high geometric resolution is an important factor for the application of object-based analysis methods in volcano monitoring [21,24]. Furthermore, coarser resolution makes texture analysis via GLCM more difficult because there is less contrast in the image [24]. However, when classifying small-scale volcanic hazards, e.g., fine lahar deposits within watercourses, even PlanetScope reaches its limits. Here, even better resolution systems (e.g., WorldView (0.30 m), Pléiades (0.50 m)) are useful. Remote Sens. 2022, 14, x FOR PEER REVIEW 23 of 34 Figure 12. Validation of lava area derived from OBIA methodology based on TanDEM-X drainage model (red) and manually digitized lava area (light blue). The manual digitization took place based on the PlanetScope post-scene (4 May 2019). In the lower area, the edge areas were manually captured in a slightly coarser way. In contrast, the lava delta in the sea was mapped better by manual digitization. In the upper area, the potential misclassifications in the E of the lava flow were omitted. In addition, there is the connection of the small lava flow in the WNW to the North Crater. In the sub-areas of the unvegetated area, manual digitization was difficult to perform because there was little contrast with the background. Accordingly, the optical validation of the lava flow from the OBIA method is given only qualitatively.
The OBIA methodology can be applied to other tropical volcanoes. It is important that the study area around the volcano has vegetation, as this allows potential dNDVI changes to be measured [23]. Application of the methodology in areas with little vegeta- Figure 12. Validation of lava area derived from OBIA methodology based on TanDEM-X drainage model (red) and manually digitized lava area (light blue). The manual digitization took place based on the PlanetScope post-scene (4 May 2019). In the lower area, the edge areas were manually captured in a slightly coarser way. In contrast, the lava delta in the sea was mapped better by manual digitization. In the upper area, the potential misclassifications in the E of the lava flow were omitted. In addition, there is the connection of the small lava flow in the WNW to the North Crater. In the sub-areas of the unvegetated area, manual digitization was difficult to perform because there was little contrast with the background. Accordingly, the optical validation of the lava flow from the OBIA method is given only qualitatively. Figure 13. Lava areas derived from PlanetScope and Sentinel-2 data for the Karangetang II eruption event. Both lava flows were classified under the same conditions using the OBIA change detection methodology based on dNDVI and NHI. A pre-(9 May 2019) and post-scene (4 December 2020) were used for classification and modelling using Sentinel-2. The results with the 10 m Sentinel-2 imagery are also satisfactory. However, marginal areas to the N of the surface are recorded much more coarsely. In the S of the lava field larger differences to the PlanetScope derived lava area were detected. Due to the lower resolution, the GLCM differences could possibly not detect detailed differences between the objects, which is why the lava areas appears somewhat coarser here.
The spectral resolution of PlanetScope is a limitation for monitoring active volcanoes. Due to the missing SWIR or TIR bands hotspots cannot be detected. Work is currently underway to improve the spectral resolution of PlanetScope. In the future, the Doves will be equipped with an 8-band sensor system, with the addition of a Coastal Blue, Green II, Yellow, and Red Edge band to the existing bands, but unfortunately no SWIR channels [44]. For other sensors with higher spectral resolution (e.g., Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS)), spatial resolution limits the exact localization of thermal anomalies [17]. Thus, via the combination of geometrically high-resolution PlanetScope data and spectrally high-resolution IR data, the full potential of both data sources can be combined. Sentinel-2 and Landsat-8 provide a hybrid solution here, combining IR bands with good geometric resolution, allowing hotspots to be located spatially (e.g., NHI) [16]. By integrating these data with PlanetScope, the spectral limitations of PlanetScope can be circumvented.
Another major advantage of PlanetScope data is the daily coverage of the entire earth's landmass. Due to the constellation with more than 150 Doves, this high temporal resolution can be achieved [44]. Especially in volcano monitoring, this is a drastic increase compared to the otherwise used optical sensor systems such as Sentinel-2 and Landsat. Due to the highly dynamic nature of volcanic hazards, a repetition rate of several days prevents the detection of fast-moving changes [26,50]. Individual lava flows of different eruptive phases separated in time can be distinguished [22,24,61,66]. However, differentiating lava flows within a single eruptive phase is difficult. Based on the high temporal resolution of PlanetScope data, it is possible to analyze individual lava flows and their dynamics during an eruptive phase, especially in conjunction with the developed OBIA Figure 13. Lava areas derived from PlanetScope and Sentinel-2 data for the Karangetang II erup-tion event. Both lava flows were classified under the same conditions using the OBIA change detection methodology based on dNDVI and NHI. A pre-(9 May 2019) and post-scene (4 December 2020) were used for classification and modelling using Sentinel-2. The results with the 10 m Sentinel-2 imagery are also satisfactory. However, marginal areas to the N of the surface are recorded much more coarsely. In the S of the lava field larger differences to the PlanetScope derived lava area were detected. Due to the lower resolution, the GLCM differences could possibly not detect detailed differences between the objects, which is why the lava areas appears somewhat coarser here.
The spectral resolution of PlanetScope is a limitation for monitoring active volcanoes. Due to the missing SWIR or TIR bands hotspots cannot be detected. Work is currently underway to improve the spectral resolution of PlanetScope. In the future, the Doves will be equipped with an 8-band sensor system, with the addition of a Coastal Blue, Green II, Yellow, and Red Edge band to the existing bands, but unfortunately no SWIR channels [44]. For other sensors with higher spectral resolution (e.g., Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS)), spatial resolution limits the exact localization of thermal anomalies [17]. Thus, via the combination of geometrically high-resolution PlanetScope data and spectrally high-resolution IR data, the full potential of both data sources can be combined. Sentinel-2 and Landsat-8 provide a hybrid solution here, combining IR bands with good geometric resolution, allowing hotspots to be located spatially (e.g., NHI) [16]. By integrating these data with PlanetScope, the spectral limitations of PlanetScope can be circumvented.
Another major advantage of PlanetScope data is the daily coverage of the entire earth's landmass. Due to the constellation with more than 150 Doves, this high temporal resolution can be achieved [44]. Especially in volcano monitoring, this is a drastic increase compared to the otherwise used optical sensor systems such as Sentinel-2 and Landsat. Due to the highly dynamic nature of volcanic hazards, a repetition rate of several days prevents the detection of fast-moving changes [26,50]. Individual lava flows of different eruptive phases separated in time can be distinguished [22,24,61,66]. However, differentiating lava flows within a single eruptive phase is difficult. Based on the high temporal resolution of PlanetScope data, it is possible to analyze individual lava flows and their dynamics during an eruptive phase, especially in conjunction with the developed OBIA change detection methodology (Section 4.3). The daily coverage of PlanetScope allows an almost daily (depending on cloud cover) mapping of the lava flows, which in turn provides information about the eruption activity. Hence, especially in near-real-time monitoring of volcanic activity, useful information can be generated [13].

Conclusions
In this work, a new change detection methodology for object-based mapping of lava flows was developed using PlanetScope data. In addition, ash and tephra deposits could be detected using difference analysis of the Normalized Difference Vegetation Index (dNDVI).
The developed change detection method combines the advantages of the very high resolution (VHR) PlanetScope data in an object-oriented approach. By combining the data with a drainage network model, Normalized Hotspot Indices (NHI) hotspot pixels (derived from Sentinel-2 and Landsat-8 data), and dNDVI analyses, lava flows could be mapped even in unvegetated areas where volcanic deposits and the background have high spectral similarities. This is a major added value compared to methods described in the literature for detecting volcanic deposits (e.g., lava flows) using PlanetScope data based solely on NDVI changes. The object-oriented change analysis obtained promising results, despite some generalizations in the modelling of potential lava areas.
The PlanetScope data turned out to be highly useful in the analysis of volcanic hazards. Due to the high geometric resolution, even small-scale change processes can be detected. In addition, the high temporal resolution allows for the monitoring of highly dynamic processes in volcanic environments with high eruptive activity. This represents a significant advantage over comparable remote sensing systems that are commonly used for volcano monitoring. Due to the lack of bands in the shortwave or thermal infrared, no thermal information about volcanic activity can be collected with PlanetScope data. Using a multisensor analysis, e.g., combining PlanetScope with infrared data, these spectral limitations of PlanetScope can be circumvented. The data sources can thus complement each other and exploit the full potential of Earth observation in volcano monitoring. The results of this work and the aforementioned advantages of the data show that PlanetScope is very well suited for monitoring volcanic hazards and active volcanoes. The derived lava areas or ash and tephra deposits can serve as a basis or additional data for further work and thus support natural disaster management in the risk analysis of volcanic hazards. In addition, the developed methodology based on PlanetScope data can make an important contribution to rapid damage mapping in disaster areas in the future, e.g., in the context of the International Charter 'Space and Major Disasters'.  (b) Areas based on thermal overprinting (NHI area). This area was considered to be cally overprinted and thus is certainly covered with lava deposits. It forms the basis for mo the lava flows in the areas with texture and brightness changes. Areas in the "change are were classified using the GLCM and brightness thresholds and are within the "potential lav class. Based on the NHI surface, region-growing modelling is performed into the "change a Figure A2. (a) Potential lava areas based on distance to the outflow model (Copernicus DEM). All areas in the "potential lava area" class are within 45 m of the nearest point to the extracted drainage flows. (b) Areas based on thermal overprinting (NHI area). This area was considered to be volcanically overprinted and thus is certainly covered with lava deposits. It forms the basis for modelling the lava flows in the areas with texture and brightness changes. Areas in the "change area" class were classified using the GLCM and brightness thresholds and are within the "potential lava area" class. Based on the NHI surface, region-growing modelling is performed into the "change area".