## 4 Methods  
### 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 PlanetScope 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 [PERSON] (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).