Eruption dynamics of Anak Krakatau volcano (Indonesia) estimated using photogrammetric methods

Analyzing video data from an uncrewed aerial vehicle (UAV) of two short-lived dome building events at Anak Krakatau volcano (Indonesia), we determine vertical and horizontal movements of the dome surface prior to explosions, as well as initial eruption velocities and mass eruption rates via automated feature tracking and other photogrammetric methods. Initial eruption velocities and mass eruption rates are estimated as a proxy for eruptive strength. Eruptive strength is found to correlate with deformation magnitude, i.e., larger pre-explosion surface displacements are followed by both higher initial eruption velocities and mass fluxes. In accord with other studies, our observations can be explained by an overpressure underneath the dome’s surface. We assume that the dome seals the underlying vent efficiently, meaning that pre-explosion pressure build-up controls both deformation magnitude and eruptive strength. We support this assumption by a simple numerical model indicating that pre-explosion pressure increases between 8 and 16 MPa. The model further reveals that the two events vary significantly with respect to the importance of lateral visco-elastic flow for pressurization and deformation. The video sequences also show considerable variations in the gas release and associated deformation characteristics. Both constant and accelerating deformation is observed. Our case study demonstrates that photogrammetric methods are suitable to provide quantitative constraints on both effusive and explosive activity. Future work can build on our or similar approaches to develop automated monitoring strategies that would enable the observation and analysis of volcanic activity in near real time during a volcanic crisis.


Introduction
Over 750 million people live in close proximity to active volcanoes (within 100 km), endangering their homes and lives (Cottrell, 2015). In the last 400 years, almost 260,000 lives were claimed by volcanic activity, partly due to the exceptionally high and constantly increasing population density in volcanic regions. To some extent, the latter results from the outstanding fertility of volcanic soils, which are thus often used as agricultural land. This leads to an increased risk for future volcano-related human casualties (Schmincke 2004).
To be able to evacuate endangered regions quickly and to save lives, the prediction of volcanic eruptions is crucial. However, this requires extensive monitoring of active and potentially dangerous volcanoes. Traditional monitoring techniques such as measuring seismicity (McNutt 2005), acoustic emissions (Johnson and Ripepe 2011;Fee and Editorial responsibility: M.R. James Matoza 2013), degassing (Burton et al. 2015), and geodesy (Lisowski et al. 2008) need the installation of expensive equipment such as seismometers, gas detectors, and tiltmeters close to often inaccessible active craters. The associated high costs, the considerable danger to personnel, and challenging logistics have prevented the development of extensive monitoring networks on many active volcanoes so far. This might change in the future due to innovative low-cost remote-sensing measurement techniques such as photogrammetry, enabling the safe and rapid acquisition of data even in inaccessible areas (Wakeford et al. 2019).
For over two decades, photogrammetric methods have been used to study a large variety of volcano-related phenomena. Pioneer work was even carried out much earlier by Chouet et al. (1974) who analyzed analogue movies to estimate the bulk and time-dependent properties of two eruptions at Stromboli volcano (Italy), including particle and gas velocities. More than 20 years later, Yamashina et al. (1999) measured the deformation at Unzen volcano (Japan) by visually comparing pairs of sequential photographs. Watts et al. (2002) combined photographs and video recordings with theodolite and electronic distance measurements to describe the chronological evolution of the lava dome at Soufrière Hills volcano (Montserrat) throughout the first episode of dome growth (November 1995to March 1998. More recent studies have used photogrammetry to constrain the geometry of an active vent (Dürig et al. 2015a) or to investigate lava flow extrusion , flank deformation and stability (Baldi et al. 2008;Saito et al. 2018), morphological changes in the crater area Schmid et al. 2021), as well as lava dome growth (Johnson et al. 2008;Major et al. 2008;Schilling et al. 2008;Diefenbach et al. 2012Diefenbach et al. , 2013Walter et al. 2013;Salzer et al. 2016;Zorn et al. 2020).
Lava dome growth has been examined by comparing sequential digital elevation models (DEMs) (Schilling et al. 2008;Diefenbach et al. 2012Diefenbach et al. , 2013 or by tracking features in the imagery (Johnson et al. 2008;Major et al. 2008;Walter et al. 2013;Salzer et al. 2016;Zorn et al. 2020). For instance, Major et al. (2008) manually tracked features of the lava dome at Mount St. Helens and showed an almost logarithmic decrease in the extrusion rate over 14 months, which correlated with trends in seismicity and geodetic deformation. Combining feature tracking with seismic data, Johnson et al. (2008) were able to attribute episodic long-period earthquakes at Santiaguito volcano (Guatemala) to abrupt sub-vertical surface displacements of the crater-filling dome. A similar study by Salzer et al. (2016) revealed a correlation between short-term meterscale downward displacements of the dome at Mount St.
Helens and seismic events, in both timing and amplitude. Walter et al. (2013) tracked temperature features in timelapse infrared images and found a temporal correlation between dome-deformation episodes and explosive activity at Volcán de Colima (Mexico). By combining optical feature tracking with thermal imagery, Zorn et al. (2020) were able to distinguish between the slow radial growth of the Caliente lava dome (Guatemala) and a narrow rapid lava flow extruded from the dome's summit. Instead of using ground-based data like the previous feature-tracking studies, Zorn et al. (2020) evaluated imagery recorded by an uncrewed aerial vehicle (UAV).
Traditionally, photogrammetric volcano monitoring has been based on imagery recorded by satellite-, crewed aircraft-, or ground-based instruments. Satellite data, for example, provide long-term synoptic observation of the investigated area usually at a comparatively low temporal resolution (James et al. 2020). The advantage of a more frequent data acquisition provided by crewed aircraft surveys is limited by significantly higher costs along with risk to the flight crew from geo-hazards (e.g., volcanic plumes; Melita et al. 2015). Ground-based data are characterized by both a high temporal and a high spatial resolution (Johnson et al. 2008;Carr et al. 2019). However, fixed measurement stations are likely to be destroyed by volcano-related phenomena, while handheld image recording depends on the accessibility of the study area and can pose a significant risk to the operating scientist. UAVs offer a safe alternative by allowing the rapid acquisition of imagery at a spatial and temporal resolution comparable to ground-based data (Saito et al. 2018;Wakeford et al. 2019;Schmid et al. 2021). Since 2010, UAVs can operate autonomously using a generalized GNSS-driven (Global Navigation Satellite System) autopilot system (Meier et al. 2011).
This study presents photogrammetric analyses of UAVbased video data recorded at Anak Krakatau volcano (Indonesia) in October 2018. The video sequences are used to examine both effusive and explosive activity, while the focus is on the quantitative analysis of pre-explosion dome deformation within the crater area. The study's main objective is to analyze the effusive phase with regard to both a spatial and a temporal deformation pattern. For this purpose, a 2D deformation vector field is derived from the imagery by applying the free research software Environmental Motion Tracking (EMT, Schwalbe and Maas 2017; https:// tu-dresd en. de/ geo/ emt). Additionally, the explosive phase is investigated by estimating both initial eruption velocity and mass eruption rate as a measure of eruptive strength. Finally, we discuss the strengths and limitations of the physical models tested and suggest potential improvements for the future under the light of our findings.

Study area
The Krakatau complex is located in the Sunda strait between Java and Sumatra, Indonesia (Fig. 1), and consists of four islands: Sertung, Panjang, Rakata, and Anak Krakatau ("child of Krakatau" in Indonesian). The latter represents the only currently active volcanic cone of the complex, which reached sea level in August 1927 (Harjono et al. 1989). Apart from Anak's recent emergence, at least two destructive historical eruptions contributed to the present island configuration (Effendi et al. 1986).
Anak Krakatau's frequent activity poses a significant threat to the coastal population of Java and Sumatra. Thus, a permanent multi-parameter monitoring system was installed (Hoffman-Rothe et al. 2006), consisting of five stations in the Krakatau archipelago and a volcano observatory on the western coast of Java. The monitoring network provides geophysical, environmental, and geochemical data including seismicity, GPS-based deformation, electromagnetics, soil temperature, weather, sea level, video, as well as measurements on fumarolic gas emission.
About 2 months following the acquisition of the video data used in this study, a major collapse event took place at Anak Krakatau volcano. On December 22, 2018, at 13:55:49 UTC (Walter et al. 2019), Anak's western flank collapsed into the sea, causing a tsunami with maximum run up heights of more than 14 m above mean sea level (update from January 2, 2019, Tsunami and Disaster Mitigation Research Center (TDMRC), http:// tdmrc. unsyi ah. ac. id/). As a result, 431 people were killed, 7200 injured, and 46,646 displaced, while the coastal infrastructure was significantly damaged (press release dated December 30, 2018, 12:48 Volcanoes are visualized by black triangles; SFZ labels the Suma-tra Fault Zone. The detailed view of the Sunda strait emphasizes the location of the Krakatau volcanic complex on a N20° volcanic line extending from Mt. Rajabasa to Panaitan island. The figure is modified from Deplus et al. (1995) UTC, Indonesian National Disaster Management Authority (BNPB), https:// bnpb. go. id/).
The instability of the volcanic edifice and thus the potential danger of a landslide-induced tsunami were already noted in numerous studies several years prior to the sector collapse event (Camus et al. 1987;Deplus et al. 1995). In hindsight, month-long precursors were identified by Walter et al. (2019), who evaluated multi-parametric ground-based and space-borne data recorded between January 2018 and January 2019. The authors identified an elevated stage of activity throughout 2018, while infrared data indicated the onset of a new intense eruptive phase on June 30, lasting until the December 22 sector collapse.
The video sequences analyzed in this study were recorded during the new intense eruptive phase and show persistent Strombolian to Vulcanian activity, which was characteristic for this period. The persistent explosive activity resulted in the deposition of eruptive products on the summit as well as on the western and southern flanks of the volcanic edifice. As a consequence, the island's surface area and the load acting on its summit (and especially on its southern flanks) increased progressively. InSAR time series revealed that the southwestern and southern flanks were already slowly subsiding and moving westward in January 2018 (see Walter et al. 2019). Over one third of the island was affected by precollapse deformation with a total displacement of 3.36 m. Considering satellite-borne and UAV-based imagery, Walter et al. (2019) conclude that only 45-60% of the previously deforming subaerial flank actually failed during the collapse event, reducing the island's height from about 320 to 120 m. The result was a steep amphitheater morphology enclosing a deep valley on the southwestern sector of the edifice, along with shifted coastlines due to the deposition of new volcanic material.

Field campaign
The video sequences analyzed in this work were recorded by Martin Rietze in October 2018, showing the active cone of the Krakatau complex prior to the December 22 flank collapse. Seven eruptions were filmed, the first six of them occurred between 08:28 and 09:17 a.m. on October 24, while the last one took place 1 day later at 11:45 a.m. local time (Rietze 2019, pers. comm.). Excerpts of the October 24 video sequences are available at http:// mriet ze. com/ web16/ Kraka tau18 air_ yout. htm.
A Mavic Pro drone with a sensor size of 6.3 × 4.7 mm was used for data acquisition. When utilized in video mode, the 4:3 format is cropped to 16:9, keeping the entire width and reducing the height of the image (Rietze 2019, pers. comm.). This results in a used sensor size of 6.3 × 3.5 mm and an effective pixel size of 1.64 µm. Twenty-five frames were recorded per second with an image resolution of 3840 × 2160 pixels.
To provide usable data for photogrammetric analysis, the individual video frames of an eruptive event must be as comparable as possible. For instance, camera perspective and position are preferably kept stable throughout the sequence, and it is crucial to keep the scale of frames and pixels constant. The former requires that each video is captured during a stationary hover of the UAV, which was not always the case (we note that none of the videos were recorded with a scientific goal in mind). Thus, only the videos that fulfill this requirement were considered for analysis. During recording, the drone-internal gimbal system was used for camera stabilization, providing a motor-based in situ compensation of small movements. Therefore, no remapping or post-processing was required to stabilize the data, and all frames maintained their size in terms of pixel quantity. Additionally, no changes in camera focus were carried out, neither manually nor automatically (Rietze 2019, pers. comm.). This allowed us to assume a constant focal distance value within our analyses.

Videos: deformation measurement
Due to position changes of the UAV during recording, adverse perspectives, and strong pre-explosion degassing, the photogrammetric evaluation for five of the seven recordings was found impossible. Both remaining videos show the first and the third eruption in the recording sequence, including their effusive and explosive stages. In the following, the effusive part of both recordings is described; the full videos are provided in the supplementary material (Online Resources 1 and 2).

Eruption 1
The effusive stage of the first eruption was filmed over a duration of almost 23 s from a distance of ~ 483 m, giving a spatial resolution of ~ 0.17 m per pixel. Since the video was recorded during a stationary hover of the UAV, it provides invariant object sizes along with a nearly constant image plane. The minor wind-related camera movements can be measured precisely in the imagery, allowing to remove their effect from the measured deformation. Due to the only slightly elevated, almost horizontal view of the crater area (12.30° deviation from horizontal), the image plane approximately equals a vertical plane in 3D space. This perspective is favorable since it allows us to project the measurements onto a vertical motion plane and translate pixel sizes into lengths of SI units. Additionally, low to moderate degassing allows for a clear view on the visibly deforming lava dome, reducing feature-tracking errors caused by rising volcanic gas.

Eruption 3
The effusive stage of the third eruption was recorded over almost 34 s from a distance of ~ 304 m, resulting in a spatial resolution of ~ 0.13 m per pixel. Similar to eruption 1, the UAV maintained a stationary position during recording and the camera was only affected by minor wind-induced movements. The camera perspective also provides a marginally elevated, nearly horizontal view of the crater area, deviating only 15.72° from the horizontal. As a consequence, visible deformation is restricted to an approximately vertical motion plane, allowing us to scale the measured trajectories. Degassing was slightly stronger than in case of eruption 1, increasing the probability of tracking errors due to rising volcanic gas. Yet, except for a few more distant sections, the deformation of the lava dome is visible. Since this data set was acquired from a shorter distance than the recording of eruption 1, the dome appears larger in size within the imagery. Therefore, deformation could be investigated at greater detail.

Set of photographs: DEM generation
In addition to the video sequences, a set of 27 individual oblique photographs was recorded during the October 2018 field campaign. From a subset of this imagery, we created the digital elevation model (DEM) of the volcanic edifice. DEM generation was performed using structure from motion (SfM) photogrammetry in the software Agisoft Metashape©.
SfM represents a photogrammetric method based on the simultaneous estimation of camera position and orientation as well as object geometry (structure) from a wide range of perspectives (Ullman 1979).
Initially, the resulting 3D model was not scaled or georeferenced. Therefore, it had to be transformed from its local coordinate system into the world coordinate system. This georeferencing was conducted by assigning three 3D control points that were derived from satellite images. To define the horizontal within the 3D model, we selected points as control points that were located on the coastline at sea level. This way, we could assume an identical height level for the control points, which we set to zero. The three control points, which are widely distributed around the volcano, allowed us to fix the position and orientation of the DEM in the world coordinate system. The distances between these control points were used to scale the whole DEM to metric units. Since the three control points delivered three scales (although a single one would already have been sufficient), there was a redundancy allowing us to estimate a scaling error. This scaling error is a relative error that affects distances measured from the DEM, and was estimated to be 0.1%.
The DEM (Fig. 2) was derived from a set of five oblique images plus the first frame of each video sequence, allowing us to determine the exterior camera orientation (i.e., camera position and orientation) during data acquisition. Combining exterior camera orientation and DEM, a so-called depth map was generated, which assigns the distance between camera and imaged object to each pixel. Together with the camera file (camera-specific information), the depth map serves as input for the scaling and georeferencing of deformation trajectories.

Data processing
Three main steps were followed to extract a scaled and georeferenced 2D displacement vector field from the previously described video sequences: preprocessing, processing, and postprocessing.

Preprocessing
The preprocessing workflow consists of four individual tasks: (I) Preparation of image data Frames must be extracted from the pre-explosion part of the video recording as a preparatory step to be processed by the Environmental Motion Tracking software package (EMT; Schwalbe and Maas 2017).

(II) Definition of fix points and object points
The determination of motion curves (trajectories) from image sequences is based on the assignment of identical points in successive frames. Since the imagery examined here shows a superposition of pre-explosion dome deformation and camera motion, the latter must be measured individually in order to be removed from the measured displacements. For this purpose, two sets of point clouds are defined, each within a specifically selected frame of the sequence. One point cloud must include stationary attributes (fix points) to determine camera motion (Fig. 3), while the other (object points) defines a point grid on the surface of the investigated object (Fig. 4a, b). The quantification of valid camera motion parameters requires a sufficient number of fix points framing the investigated object. For our image analysis, we assigned 3749 (eruption 1) and 4274 (eruption 3) fix points.
In contrast to fix points, object points are specified within the dynamic image section that is investigated regarding motion and deformation. The extent of this dynamic area is only accessible in the last frame of the pre-explosion image sequence, immediately prior to the explosion when deformation is strongest. Hence, the last frame is used for object point definition, allowing us to specify an object point grid entirely enclosing the fully developed lava dome (see Fig. 4a, b).

(III) Selection of relevant object points
The object point cloud is further classified into relevant and irrelevant points. Only relevant points (i.e., points inside the dome boundary; see Fig. 4c,d) are used in the subsequent processing steps for two reasons: (1) Reducing the number of points to process decreases computational effort.
(2) Rather than the point itself, it is the area around each point (i.e., a patch) that is tracked throughout the image sequence. Since the measured displacement is interpolated over this patch, notable deformation can be incorrectly suggested for points located outside the dome area. This interpolation artifact is avoided by exclusively evaluating relevant object points.

Fig. 3
Fix points (green) as specified within the first frame of the pre-explosion recording; a eruption 1; b eruption 3. Note that only the image area covered with points is shown. In both cases, the lava dome (enclosed by the red ellipse, see Fig. 4a, b for a more detailed view) is framed by fix points from three sides, as pre-explosion degassing prevents point definition above the crater area (IV) Object point transformation Unlike fix points, object points are specified within the last frame prior to the onset of explosive activity. However, tracking starts with the first image of the recording, meaning that object point coordinates must be transformed into the coordinate system of the first frame of the sequence. For this purpose, the camera motion between both frames is corrected using a transformation tool provided in EMT.

Processing
Deformation measurement is performed using EMT, a stand-alone photogrammetric evaluation software primarily designed to examine glacier flow. The program provides a workflow to analyze monoscopic time-lapse imagery, i.e., image sequences recorded from a single and stationary position in defined time intervals (Schwalbe and Maas 2017). The first step of the workflow comprises the quantification of camera motion by fix point tracking and the subsequent derivation of camera motion parameters. The resulting parameters are applied during object point measurement to remove the effect of camera motion from object point trajectories. Following their transformation from 2D image space into 3D object space, the scaled and georeferenced trajectories are exported for further processing (for a detailed description of the processing workflow, see Online Resource 3).

Postprocessing
Examination of the measured point trajectories in Fig. 5 reveals the requirement for two postprocessing steps. The first of both operations, however, only has to be applied to eruption 3.
(I) In case of eruption 3, a notable number of trajectories is considerably affected by pre-explosion degassing (especially in the NW dome area). Since these trajectories represent a superposition of deformation and gas emission, they must be excluded from the evaluation process. For this purpose, every distorted trajectory is identified and deleted manually from the data set. The result of this manual data classification is indicated in Fig. 5c. (II) Due to tracking inaccuracies, the trajectories for both events are characterized by short-period oscillations partly culminating in noise-like behavior. In order to smooth the data without altering their start points, a trailing moving average filter is employed. The filter window duration is defined to decrease near the endpoints of the input array, so that only existing elements contribute to the average calculation. To ensure comparability of results, all window lengths (durations) are initially determined for eruption 1 and subsequently adapted to eruption 3. This is done by tripling the number of frames from each window used for erup- Fig. 4 a, b Object points (red) as defined within the last frame prior to the onset of explosions; a eruption 1; b eruption 3. Note that the object point polygon covers domains beyond the actual area of interest. Following the strategy described in (III), these unnecessary object points are discarded. c, d Dome boundary (red) interpolated from manually picked coordinates along with selected (black) and discarded (gray) object points; c eruption 1; d eruption 3 tion 1 when analyzing eruption 3, because the analyzed frames of eruption 3 have a temporal resolution three times finer than those of eruption 1 (see section Object motion tracking in Online Resource 3).

Pre-explosion deformation
The object point trajectories and the manually picked dome boundary visualized in Fig. 5 suggest that both lava domes are comparable in size prior to the onset of explosive activity. Concerning eruption 1, the dome attains a horizontal extension of about 27.0 m, while it achieves roughly 6.7 m in vertical direction. Similarly, a horizontal dimension of more than 29.3 m is reached in case of eruption 3 with a maximum vertical extent of approximately 6.1 m. Considering the maximum object point displacements during both effusive stages, though, significant differences in deformation magnitude are exposed. In case of eruption 1, a maximum total displacement of 4.07 m is observed, while eruption 3 exhibits a notably higher maximum value of 7.20 m (displacement error is of order centimeters, see Online Resource 3). Similar distinctions in magnitude are found for both the horizontal and the vertical displacement component (see Table 1).
In order to investigate whether the dome margins behave differently to the central area, we divide each lava dome into three different parts with equal horizontal width for individual evaluation. We decided to group trajectories that way since we only measured the 2D displacement in the image plane, meaning that the horizontal displacement directed towards or away from the camera remains unknown. Hence, statements we can make about the spatial dome deformation are limited to the two dimensions of the image plane and only represent tendencies. Still, we can use these tendencies to estimate a simplified spatial deformation pattern.
Grouping of trajectories is carried out prior to deformation, i.e., using the trajectory start point coordinates (Fig. 6). In case of eruption 1, the dome is divided into a southwestern (SW), a central, and a northeastern (NE) area, each containing a similar number of trajectories (~ 200). However, for The black line marks the manually picked dome boundary, while the acronyms of the cardinal points specify the orientation of the image plane. Note the noisy oscillation of trajectories due to tracking inaccuracies. In case of eruption 3, strong degassingrelated distortion is observable especially within the NW dome area. The affected trajectories are deleted manually, resulting in the data set shown in c. Start points of discarded trajectories are marked by red points eruption 3, the trajectories are distributed more irregularly to the three distinct dome regions. This results from strong pre-explosion degassing mainly restricted to the marginal dome areas, especially the northwestern one. Most trajectories from these areas have been removed since they were significantly distorted (Fig. 5c). For both events, the average and maximum displacements are computed for each dome region.

Eruption 1
The dome's final horizontal dimension is approximately four times its height. Still, the observed deformation mostly consists of vertical displacement, while the horizontal component only contributes marginally (Fig. 7; Table 1). The maximum change in elevation measured for a single object point is 4.03 m, while maximum values of 1.31 m and 1.36 m are found for the eastern and western directions, respectively.
The temporal development of all deformation components (Fig. 7a-c) indicates an acceleration of dome deformation, which is confirmed by the increasing total expansion velocity shown in the inset of Fig. 7c. Moreover, a general spatial deformation pattern is discernible (Fig. 7d-f): Vertical deformation is exclusively positive (upward), with the maximum average uplift of 2.00 m occurring in the central dome region. Towards the dome margins, uplift decreases to mean values of 1.54 m (west) and 0.69 m (east), respectively (Fig. 7e). The horizontal displacement component, though, shows similar deformation magnitudes between 0.84 and 0.90 m for all dome areas, but in different directions. Both the eastern and central parts move to the east, whereas the western dome region moves in western direction. Overall, horizontal displacement is thus directed away from the dome center (Fig. 7d). Considering the amplitude ratio of vertical and horizontal deformation, we find that the horizontal proportion of total displacement is higher for the dome margins than for the central area. This must be the case since the horizontal motion of the dome center is restricted more than that of the dome margins. Hence, the dome center is dominated by uplift, and the outward-directed horizontal component increases towards the dome's marginal parts (Fig. 7f).

Eruption 3
The deformation measured during eruption 3 exhibits a notably higher magnitude than the deformation observed for eruption 1 (Fig. 8; Table 1). Nevertheless, the deformation process is similarly dominated by vertical displacement as reflected by a maximum object point translation of 6.20 m in vertical direction. In comparison, the maximum values obtained for the eastern and western components yield 3.55 m and 4.80 m, respectively.
In contrast to eruption 1, Fig. 8a-c suggest that the deformation process is not dominated by a continuous acceleration but shows two distinct phases. Velocities (see inset of Fig. 8c for total velocity) and accelerations (not shown here) show that the first ~ 6 s of the observation period are defined by accelerating deformation, while the remaining time span exhibits a constant deformation rate. The spatial deformation pattern, though, is similar to eruption 1 with the central uplift transforming into outwardmoving dome margins (Fig. 8d-f). However, the dome's central area shows a weaker and less definite trend to the east, along with a comparatively smaller uplift amplitude.

Explosive activity
Analyzing the explosion part of each video, we assess initial eruption velocities as well as mass eruption rates to identify a possible relation between pre-explosion deformation and explosive activity.

Initial eruption velocity
The estimation of initial eruption velocities is based on manual pulse evolution tracking. An automated approach is inapplicable mainly due to the insufficient temporal resolution of the video sequences. For the recording interval of 0.04 s between successive frames, most pulses are strongly de-correlated, meaning that the probability for incorrect point-matching is high. This error accumulates with each subsequent measurement and is further increased by ambiguous textures of the eruption cloud. Manual pulse evolution tracking is executed using EMT. Beginning with the frame where a pulse first appears, we trace the pulse tip throughout three consecutive images (Fig. 9). For this purpose, we place an object point on the pulse tip in each frame, giving two inter-frame velocity values v 1 and v 2 . While v 1 corresponds to the initial velocity of the pulse, v 2 is only determined to verify the reliability of v 1 . Since the pulse tips rise above the volcanic edifice, their 3D coordinates cannot be determined by ray intersection with the DEM surface. Thus, we measure the metric size of a pixel on the DEM and apply the result to scale the pulse tip distances between consecutive frames. Dividing these distances by the recording interval of 0.04 s, velocities are obtained. However, similar to deformation measurement, the horizontal velocity component perpendicular to the image plane remains unknown.
Analyzing nine pulses of eruption 1 as well as five pulses of eruption 3, we find that eruption 3 exhibits slightly higher velocities, with averages of 118 m/s (eruption 3) and 107 m/s (eruption 1; see Table 2). Moreover, as indicated by the larger maxima and lower minima, the velocities measured for eruption 3 are distributed over a wider range.

Mass eruption rate
The mass eruption rate (MER) is defined as the amount of volcanic material (i.e., tephra and gas) released into the atmosphere per unit time. It is considered to be one of the most important parameters for hazard assessment of explosive volcanic eruptions (Woods 1988;Glaze and Baloga 1996;Mastin 2014). In order to estimate the MER of individual eruption pulses from video sequences, Dürig et al. (2015b) developed a photogrammetry-based method termed the pulse velocity-derived model (PVDM).
Within the PVDM, the volumetric flow rate of a pulse is approximated as that of an ash package through a cylindrical cross section. In addition, minimum (Q tot_min ) and maximum (Q tot_max ) estimates of MER are obtained by assuming an exponential decay or a step function for the peak mass flux (Q(t I )), respectively (for more details, see Dürig et al. (2015b)). To convert the volumetric flow rate into the mass eruption rate, information about the pulse density is required. Hence, each pulse is analyzed at the so-called The manually picked dome boundary is shown by the black line, while the acronyms of the cardinal points specify the orientation of the image plane. In case of eruption 3, note the absence of trajectories in the northwestern dome region. These trajectories have been removed since they were strongly affected by degassing transitional level, where a negatively buoyant jet transforms into a convective buoyant plume. At this pulse-specific height, pulse density approximates the density of ambient air, allowing all parameters necessary for MER estimation to be derived.
The application of the pulse velocity-derived model to our eruption sequences requires the assessment of several parameters mostly based on distances picked in the images. Similar to velocity estimation, these picks cannot be georeferenced using the DEM. Therefore, we again use the metric size of a pixel to scale our picked distances.
Both analyzed eruptions exhibit a high pulsation frequency with most of the pulses being obscured by subsequent ones. As a result, only a few individual pulses allow the estimation of all required PVDM parameters. In case of eruption 1, only two pulses prove suitable for MER assessment, while the PVDM can be applied to four pulses of eruption 3. We note that the two pulses in eruption 1 differ significantly in strength, whereas for eruption 3, all distinguishable pulses are approximately comparable.
Deploying the PVDM to our video sequences, mass eruption rates of order 10 3 -10 4 kg/s are obtained (see Table 3). The minimum and maximum estimates (i.e., Q tot_min and Q tot_max ) for both eruptions give a range from 0.7 × 10 4 to 3.6 × 10 4 kg/s for the total mass flux. However, the individual values for Q tot_min and Q tot_max differ significantly between both events. For eruption 3, both Q tot_min and Q tot_max exceed the corresponding values for eruption 1 by more than a factor of three. Similarly, our analyses reveal stronger preexplosion deformation as well as higher initial eruption velocities for eruption 3. displacement; c total displacement. All measured trajectories are displayed, with their gray value scaling with the horizontal distance between the corresponding start point and the dome center. Shades close to black indicate trajectory start points near the dome's horizontal center point, whereas trajectories visualized in light shades originate close to its margins. Addi-tionally, mean values calculated for the three different dome areas and for the entire lava dome are shown. The inset of c shows the average total expansion velocity of the dome. Please note that the axes in this figure are scaled such that they can be easily compared to the results of eruption 3 shown in Fig. 8. d-f Average displacement for all dome parts and deformation components. Please note that amplitudes are only comparable within each deformation component

Pre-explosion deformation
Comparison to previous studies Watts et al. (2002) report the rapid extrusion of a pancake-shaped lobe at Soufrière Hills volcano, Montserrat, on December 25, 1996. The lobe was defined by a circular plan-form along with a slightly raised central summit. Both lava domes examined in this study exhibit a similar, although more irregular profile. We note that since they were only recorded from a single perspective, no statement can be made about the shape of both domes in respect of the horizontal plane. Still, the similarity identified in two dimensions can be investigated in further detail by comparing all three structures regarding their width-to-height-ratios (WH ratios). On December 28, 1996, the lobe at Soufrière Hills attained a horizontal extension of 155 to 169 m, while its height ranged from 21 to 33 m, giving an average WH ratio of 6.0. In case of both lava domes analyzed at Anak Krakatau, evaluation of the manually picked dome boundary yields smaller WH ratios of 4.0 (eruption 1) as well as 4.8 (eruption 3). Thus, the lateral flow of lava appears to have been lower at Anak Krakatau, suggesting that the magma viscosity was higher than the one at Montserrat during the observations described above. However, considering the different time scales of the two extrusion events, they may certainly also have been of similar magnitude.
At Santiaguito volcano (Guatemala), Johnson et al. (2008) measured dome uplift magnitudes between 20 and 50 cm occurring within a few seconds and coinciding with displacement; c total displacement. All measured trajectories are displayed, with their gray value scaling with the horizontal distance between the corresponding start point and the dome center. Shades close to black indicate trajectory start points near the dome's horizontal center point, whereas trajecto-ries visualized in light shades originate close to its margins. Additionally, mean values calculated for the three different dome areas and for the entire lava dome are shown. The inset of c shows the average total expansion velocity of the dome. d-f Average displacement for all dome parts and deformation components. Please note that amplitudes are only comparable within each deformation component the onset of lateral flow. The latter was directed from the dome center towards the crater periphery, which is similar to the motion pattern identified in this study. Since dome uplift was accompanied by frequent explosive eruptions emerging once or twice per hour, Johnson et al. (2008) suggest a close relation between surface deformation and explosive activity. In the current study, we find an equivalent connection for both Anak Krakatau events expressed by dome deformation prior to explosions.
In order to explain their observations, Johnson et al. (2008) assume that each upheaval event is preceded and consequently triggered by the steady accumulation of gas beneath the viscoelastic lava dome. The dome itself has become impermeable during the 20-40-min inter-eruption interval, allowing gas to accumulate in fissures and/or voids beneath the dome's edifice and in the magma conduit periphery. At critical pressurization, the dome starts to detach and accelerate upward. Within a few seconds, its 20-80-m-thick edifice inflates by several tens of centimeters, causing high strain rates and thus brittle failure of segments of the lava carapace. This leads to explosive gas emission. Furthermore, Johnson et al. (2008) suggest that the lateral dome flow concomitant with uplift is gravity-driven and facilitated by a temporary decoupling of the dome from the underlying Fig. 9 Manual pulse evolution tracking for an exemplary pulse of eruption 1 (pulse front indicated by solid line). The pulse fronts of other pulses for which the picking is not shown here are indicated by dashed lines. The position of the pulse tip is picked in three consecutive frames using object points. The first pick is set in the frame where the pulse first appears. From the resulting velocity values, v 1 represents the initial velocity of the pulse, which we need for our analysis. The second velocity value v 2 is only determined to verify the reliability of v 1 surface. Scharff et al. (2014) showed that the dome exhibits a damped oscillatory motion with opening and closing paths for gas to escape. A comparable but slightly different scenario is suggested by Maeno et al. (2013) to explain a Vulcanian explosion of Shinmoedake volcano (Japan) following a dome-building stage. Similar to Johnson et al. (2008), the dome is assumed to seal the underlying vent efficiently, preventing the emission of magmatic gas for the most part. As a consequence, gas accumulates inside the dome or near its surface. In this case, however, an isolated gas pocket or a magma-gas mixture is suggested instead of gas-filled fissures and/or voids. By the time the accumulated pressure exceeds the bulk tensile strength of the lava dome, an explosion occurs.
The imagery of Anak Krakatau examined here reveals a major role of gas-induced pressurization and is consistent with the model concepts provided by Johnson et al. (2008) and Maeno et al. (2013). In the next section, both scenarios are combined to evaluate our findings.

Dome deformation model
For each video recording available, we only have a short preexplosion sequence as well as the explosive event itself, with the footage ending at some time into the eruption. We cannot assess how much of the dome disintegrated throughout the explosive phase. Most probably, some of the material released during the eruption fell back into the vent area, but its thickness is most likely insignificant as the video sequences suggest trajectories that carried the ejecta away from the vent.
The deformation pattern identified, comprising a transition from centered upward to lateral outward displacement, suggests that the dome's core consists of hot magma, which is deformed in a ductile way. The dome largely prevents the emission of exsolved magmatic gas by acting as an efficient seal for the underlying conduit material. Consequently, significant pressure can build up underneath the dome's surface. The combined effects of pressure accumulation due to magmatic gas exsolution and magma inflow into the core of the dome result in the observed surface deformation. The point in time when the accumulated pressure exceeds the bulk tensile strength of the dome-forming material eventually marks the onset of explosive activity.
However, since we found considerable differences between the pre-explosion degassing and deformation characteristics of the two eruptions, we need to refine the above deformation model. The average eastern and western displacements for eruption 3 are about three times as large as the ones for eruption 1 (see Table 1). The vertical deformation measured in the central dome part, however, is only 1.5 times larger for eruption 3. Figures 7 and 8 further reveal that eruption 1 shows significant acceleration in all directions during the last 5 s prior to the onset of explosive activity, which is lacking in case of eruption 3 (compare insets of Figs. 7c and 8c). We also note that during eruption 3, notable gas emission was observed in the western part of the dome, whereas degassing from the dome was nearly absent during eruption 1. All these differences indicate that for the two eruptions, deformation was not dominated by the same mechanism.
In the following, it is assumed that the pressure needed to exceed the bulk tensile strength of the dome-forming material is similar for both events. The strong acceleration phase observed prior to the explosive onset of eruption 1 can be attributed to an effectively sealed dome. In contrast, the non-accelerated deformation during eruption 3 indicates that pressure did not increase as fast as during eruption 1 but over a longer time span. We therefore infer that for eruption 1, the deformation of the dome was mainly the result of pressurized gas, while for eruption 3, dome deformation was dominantly controlled by the lateral visco-elastic flow underneath. This interpretation is supported by some straightforward numerical simulations discussed below.

Numerical model verification
The numerical model simulates dome deformation due to an overpressure at some depth below the dome surface by solving for nonlinear elasticity in three dimensions. Using a deformable mesh, we approximate the plug inside the volcanic conduit by a truncated and upturned cone with a thickness of 14 m (static pressure at the bottom of this sheet ~ 0.25 MPa) as well as upper and lower radii of 28 and 24 m, respectively (Fig. 10). The plug material is defined Fig. 10 A truncated and upturned cone is simulated as a plug filling the conduit. Taking the final horizontal dimension of both lava domes into account, the cone's upper radius is estimated as 28 m, while its lower radius is set to 24 m. Based on assessments of the minimum erupted volume, a thickness of 14 m is assumed. The plug material exhibits a Young's modulus of E = 100 MPa (Chadwick et al. 1988), while the Poisson ratio is set to µ = 0.25 (Bonaccorso and Davis 1999) by a Young's modulus of 100 MPa (Chadwick et al. 1988) and a Poisson ratio of 0.25 (Bonaccorso and Davis 1999). For the boundary between the plug and the conduit wall, all displacement components are set to zero (no slip boundary condition). On the lower edge of the plug, a pressure is applied and altered until the EMT-derived maximum vertical surface deformation is obtained. Model calculations are performed with a FEM-based software called Elmer (CSC -IT Center for Science, https:// www. csc. fi/ web/ elmer).
The resulting pressures of 8.72 MPa (eruption 1) and 15.10 MPa (eruption 3) correspond to the lower range attributed to Vulcanian-type activity (Ishihara 1985;Robertson et al. 1998;Formenti et al. 2003), which seems plausible in view of the observed eruptive strength with MER ~ 10 4 kg/s. Contrary to our previous assumption, however, these pressure values differ significantly, with the value for eruption 3 being nearly twice the value for eruption 1. Together with the fact that the modeled horizontal displacements for eruption 3 are significantly smaller than those from the photogrammetric measurements (see Table 4), this strongly supports our inference that for eruption 3, lateral flow of dome material played a dominant role in pressurization, leading to explosive activity. As the numerical model disregards flowing material, pressure is overestimated for eruption 3, while horizontal deformation is underestimated. Thus, more complex modeling that includes visco-elastic flow is necessary to completely explain the measured displacement, but this is beyond the scope of this paper.

Initial eruption velocity
Since the numerical model described above indicates pressure values in the lower range of Vulcanian-type activity, it seems plausible to compare our velocity assessments to results from other studies on Vulcanian-type activity. By applying a ballistic trajectory model along with a scaling law for impact crater formation, Maeno et al. (2013) analyzed the size of impact craters created by ballistic ejecta from a Vulcanian explosion of Shinmoedake volcano (Japan). As a result, they obtained initial velocities from 240 to 290 m/s, while other studies on Vulcanian-type activity present values between 40 and 180 m/s (Ishihara 1985;Robertson et al. 1998;Formenti et al. 2003). Thus, the average initial velocities of 107 m/s (eruption 1) and 118 m/s (eruption 3) estimated for Anak Krakatau correspond to the lower range of measurements of Vulcanian activity at other volcanoes. Consequently, our velocity assessments are in good agreement with the numerically derived pressure values discussed above, which also correspond to the lower range attributed to Vulcanian activity.

Mass eruption rate
When implementing the PVDM, the main source of error is the manual parameter picking. Additionally, the limited number of evaluable pulses restrains the statements we can make about MER variability between individual pulses. It is therefore noteworthy that the results for all four mass fluxes of eruption 3 are consistent. In contrast, the two MER estimates for eruption 1 differ by an order of magnitude from each other. To further evaluate the reliability of the PVDMderived MER estimates, we compare them with findings of previous studies on the explosive activity of Anak Krakatau.
Applying remote sensing techniques, Gouhier and Paris (2019) investigated a phase of paroxysmal activity following the December 22, 2018, sector collapse of Anak Krakatau volcano. During the first 15 days after the collapse event, the authors identified three distinct main eruptive periods. While phase III was mostly defined by degassing, phases I and II formed a 7-day period of active tephra emission. For this period, the authors estimated an average mass flux of 1.2 × 10 5 kg/s. This value is considerably lower than the mass eruption rates they suggest for the first hours of the paroxysmal event, though; according to the authors, the eruption was initiated by the sector collapse and began with 40 min of Vulcanian activity at a mass flux of 9 × 10 5 kg/s. At least the subsequent 10 h showed less intense Surtseyan activity, which resulted in a reduction of mass eruption rate to 5 × 10 5 kg/s. Based on these findings, the authors conclude that MER peaked right after the collapse event and decreased drastically during the following days.
Consequently, the mass fluxes determined by Gouhier and Paris (2019) exceed those assessed in this study by more than an order of magnitude. However, the drastic decrease in MER suggests that the days following the collapse event were defined by mass eruption rates of the same order of magnitude as the PVDM results obtained here. However, since Gouhier and Paris (2019) cannot exclude phases of quiescence during the averaged 7-day period, it is difficult to draw a final conclusion. Moreover, Gouhier and Paris (2019) state that the paroxysmal nature of the post-collapse eruption contrasts with the moderate activity observed during the preceding six months. Thus, the October 2018 videos examined here show a lower degree of activity, for which mass fluxes of order 10 3 -10 4 kg/s seem plausible. Interestingly, similar mass eruption rates were found by Dürig et al. (2015b) regarding the 2010 Eyjafjallajökull eruption. By applying the PVDM to 30 pulses, Dürig et al. (2015b) obtained a mean mass flux of (2.9 ± 0.9) × 10 4 kg/s, which is only slightly higher than the average value for Anak Krakatau derived in this study (2.3 × 10 4 kg/s).

Conclusions and outlook
Photogrammetric analyses revealed that both examined lava domes attained a comparable final extension of almost 30 m in width and more than 6 m in height. In both cases, dome deformation over the ~ 30 s preceding explosive activity was dominated by uplift, which was strongest in the dome center and transformed into outward-moving dome margins. Using a simple numerical model, we confirmed that (I) the importance of lateral visco-elastic flow for dome deformation and pressurization differed significantly for both events over this time scale, and (II) the correlation between pre-explosion surface displacement and subsequent eruptive strength (expressed by initial eruption velocity and MER) can be explained by an overpressure below the sealing dome. For a dome thickness of 14 m, the model suggested pressure values between 8 and 16 MPa.
In conclusion, this study has shown that photogrammetry is suitable to quantitatively analyze both effusive and explosive volcanic activity. Future work could further optimize the acquisition of photogrammetric data to enhance accuracy and automation. For example, a horizontally aligned camera would increase the precision of the scaling and georeferencing process, while the application of at least two synchronized cameras would allow 3D trajectory measurements. In addition, by recording the imagery at a higher frame rate, pulse tips could be tracked automatically. This would not only improve the accuracy of the estimated initial eruption velocities but also further automate the workflow. By using a UAV with GNSS-driven autopilot, data acquisition itself could be automated as well. Our approaches can be used as a basis to develop automated monitoring strategies providing near real-time information on both effusive and explosive activity, which are crucial during a volcanic crisis.