The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse

Continued post-collapse volcanic activity can cause the rise of a new edifice. However, details of such edifice rebirth have not been documented yet. Here, we present 7-decade-long photogrammetric data for Bezymianny volcano, Kamchatka, showing its evolution after the 1956 sector collapse. Edifice rebirth started with two lava domes originating at distinct vents ~400 m apart. After 2 decades, activity became more effusive with vents migrating within ~200 m distance. After 5 decades, the activity focused on a single vent to develop a stratocone with a summit crater. We determine a long-term average growth rate of 26,400 m3/day, allowing us to estimate the regain of the pre-collapse size within the next 15 years. Numerical modeling explains the gradual vents focusing to be associated with loading changes, affecting magma pathways at depth. This work thus sheds light on the complex regrowth process following a sector collapse, with implications for regrowing volcanoes elsewhere. Regrowth of the Bezymianny volcano, Kamchatka developed from dispersed domes to a focussed stratocone due to loading changes beneath the edifice over the course of seven decades, according to photogrammetric observations throughout the period.

A ctive volcanoes are subject to flank instability that may cause large-scale sector collapses 1 . The collapses may occur slowly and creep-like 2 or fast and catastrophic 3 , some with precursory eruptive activity 4 and others without any apparent change in behavior 5 . The volumes mobilized by sector collapses vary considerably, from small scale 6 to large edifice collapses involving 10 8 −10 11 m 3 7 , with the larger collapses being rarer 8 . The consequences of sector collapses are often dramatic and associated with catastrophic directed blasts and pyroclastic flows, devastation over distances as far as tens of kilometers, and potential for triggering far-reaching secondary hazards such as lahars and tsunamis, with an estimated 20,000 death toll caused by historic volcano flank collapses and associated hazards 9 . The mobilization of large rock masses at the surface is felt beneath the volcano: at depth, unloading forces cause changes in magma reservoir pressure 10 , its location and geometry 11 , and are associated with compositional changes of the erupted products 12 . Simulations of magma paths following a sector collapse suggest that they may deflect 13 and bifurcate from their pre-collapse course 14 , leading to a relocation of the eruption site(s). At the surface, following a collapse the regrowth of a new volcanic cone was observed at Mount St. Helens (USA) 15 , Soufriere Hills (Montserrat) 16 , Santa Maria (Guatemala) 17 , Ritter Island (Papua New Guinea) 18 , Shiveluch (Russia) 19 , and Bezymianny (Russia) 20 , and was geologically and experimentally reconstructed 21 for prehistoric collapse events such as evidenced at Socompa (Argentina and Chile) 22 , Parinacota (Chile and Bolivia) 23 , Vesuvius (Italy) 24 , and Tungurahua (Ecuador) 25 . To our knowledge, the Bezymianny case presented here, however, is providing the first detailed chronological evolution of a rebuilding cone morphology. The rebuilding of a new cone may continue until another sector collapse occurs, with similar direction, dimension, and effects 21,26 .
Volcanic sector collapse and regrowth are recurrent and geologically inferred at volcanoes worldwide, e.g., by structural, geologic, and stratigraphic analysis 27 . Changes in the direction of feeding magma paths may favor repeated collapses, and a possible feedback mechanism between magmatic activity and topographic changes may be given 13 . The pathways of magma may adjust as a result of the stress field redistribution owing to sudden topographic changes 13,14 , causing a shift of the eruption locations and their pattern 28 . However, detailed observations of shifting and then gradually focusing vents at active volcanoes have not been documented, so far. Worldwide there are only a few cases where cone regrowth after sector collapse had been continuously monitored over decades so that direct observations of geomorphologic and structural changes are limited.
After a sector collapse, the regrowth of the volcanic edifice may start immediately or with a delay and, with few exceptions 29 , commonly nucleates in the center of the newly exposed amphitheater. Regrowth may last days to years, as particularly observed at Mount St. Helens in the 1980s 30 and 2000s 31 . Despite the episodic growth nature, long-term regrowth is assumed to be rather regular, forming characteristic dome height-to-width ratios of 0.2-0.3 30 . To study the topographic evolution of a central cone regrowing after a sector collapse, continuous long-term observations are therefore required.
Bezymianny volcano (see Fig. 1a, b), which is one of the most active andesitic volcanoes in the world 33 , experienced a sector collapse in 1956 32 . Before the collapse, Bezymianny was a conically shaped stratovolcano that reached an elevation of 3113 m similar to that of today (Fig. 1c, f). The collapse removed over 0.7 km 3 of material from the former edifice 34 and led to a catastrophic eastward directed blast eruption 32 (Fig. 1d). Immediately after this climactic episode, construction processes commenced inside the collapse amphitheater 20,35 and continued until today (Fig. 1e, f). The dome growth style was characterized as extrusive-explosive until 1977, when it changed to extrusive-explosive-effusive 36 . This transition was accompanied by a gradual change in rock chemistry. In fact, comprehensive petrographic analysis 37 showed that Bezymianny's eruptive products became continuously more mafic since 1956.
Here, we employ a unique photogrammetric data set 38 that allows an unprecedented study of the surface changes associated with a seven-decade-long period of the post-collapse edifice regrowth at Bezymianny. With high-resolution digital elevation models (DEMs) we illuminate the morphological and structural evolution of the edifice, and that volcanic activity changed progressively from destructive to constructive processes, which eventually led to the transition from a lava dome into a stratocone with a centralized vent. Ultimately, we make use of a numerical model to show how the sector collapse caused a shift in eruption center location and the edifice regrowth promoted consecutive vents centralization.

Results
Evolution after sector collapse. The first aerial survey (Supplementary Fig. 1) from 1949 depicts the pre-collapse topography of Bezymianny volcano and reveals two collapse scars that are related to Holocene activity 34,39 (Fig. 1c, Supplementary Fig. 2).
The processed data of the first post-collapse aerial survey from 1967 shows the presence of a steep headwall amphitheater (1300 m in width, 2700 m in axial length, and up to 400 m in visible depth) facing ESE. Within the amphitheater, scattered volcanic activity (see section "Vents migration and focusing") formed a broad NW-SE elongated edifice (azimuth 120°) that was formed between 1956 and 1967 35,[40][41][42] . This new edifice is separated by a NE-SW linear depression (azimuth 50°), 200 m wide and 110 m deep, that delineates two consecutive endogenous 43 domes located in the southeast (1st) and northwest (2nd) of the amphitheater, and that coalesced through time 35,41 (Fig. 2a,  Fig. 3a, b) (see section "Morphological mapping" in "Methods"). The base of this dome complex covered an area of 1975 km 2 , the dome reached an elevation of 2801 m (~600 m in relative height) and occupied a volume of 0.239 km 334 . The resulting average growth rate between 1956 and 1967 was 56,600 m 3 /day. Moreover, the 1967 survey revealed a shear lobe (i.e., a portion of highly viscous lava 43 ) that extruded from an open crater at the summit of the north-western dome that indicates a transition to exogenous growth 43 at Bezymianny. The comparison of the 1967 and 1968 aerial images shows that the lobe increased from 220 m to 380 m in diameter and from 70 m to 110 m in thickness (Fig. 2b, Fig. 3c, d). This lobe was described by eyewitnesses and named "Nautilus" 41 . Overall, the lobe might resemble the "whaleback" spine that extruded 2.5 decades after the 1980 sector collapse at Mount St. Helens 31 . The calculated differential volume between 1967 and 1968 is 0.014 km 3 , resulting in a derived average extrusion rate of 43,600 m 3 /day (see further morphometric characteristics in Table 1, Supplementary Fig. 3; see errors in Supplementary Table 1).
The subsequent aerial data from the 1976 survey revealed that new shear lobes piled up at the south-eastern sector of the dome complex (Fig. 2c, Fig. 3e, f) suggesting that the growth dynamics has not changed significantly since 1967. Two open craters at the top of these lobes indicate that explosive activity increased. The oldest of the lobes ("Oktyabr") formed between 1969 and 1973 41 , had a thickness of 60 m, and an average diameter of~450 m. The latest of them, a 70-m-thick lobe, emerged from the youngest crater, covered an area of 150 × 300 m, and was bisected by a 190m-long crease structure (i.e., negative "crack-like" topographic feature of silicic domes 44 ).
In 1977, a new, up to 180 m wide and 300 m long, portion of lava poured out onto the north-eastern flank of the dome complex ( Fig. 2d and Fig. 4a, b). In contrast to previously observed shear lobes, it was much thinner (15 m), and flowed over a comparatively longer distance toward the foot of the dome. Therefore, we can consider it as a first post-collapse lava flow at Bezymianny, which is also in agreement with eyewitness accounts 36 . We link the different characteristics of shear lobes and lava flows with two distinct exogenous growth modes: extrusion of higher viscous shear lobes and effusion of lower viscous lava flows. Hence, dome development between 1956 and 1976 was dominated by first endogenous and then extrusive growth, whereas from 1977 onwards the dome growth was mainly marked by extrusive and effusive behavior. This change in growth style will eventually lead to a significantly different shape of Bezymianny's new edifice.
Images from the 1982 overflight unveiled a semi-circular remnant of a 60-m-thick shear lobe in the central part of the dome ( Fig. 2e; Fig. 4c, d) that was destroyed during the 11.02.1979 explosive eruption 34 . In the excavated crater, we identified remnant of a lava flow and 50-m-high remnant of a cylindrical lava plug (i.e., lava extruded from a vent after being solidified 43 ), which was also partially destroyed by additional explosive eruptions. Originating from this new crater, another 700-m-long lava flow was emplaced on the eastern flank of the dome. Aerial data from 1994 signify a new semi-circular shaped, 80-m-thick remnant of a shear lobe that was extruded around the top-center of the dome. Subsequently, this lobe was pierced by multiple lava flows (maximum length 1500 m) that emanated from an excavated crater ( Fig. 2f; Fig. 4e, f). Thus, the type of eruptions remained the same since the 1980 s.
In the 2006 images, only the north-western and southern sectors of the previous dome complex are visible: most parts-including the base of the amphitheater-were buried by numerous lava flows and pyroclastic deposits ( Fig. 2g; Fig. 5a, b). The deposits reached a maximum thickness of 140 m in the western sector of the amphitheater (relatively to the 1967 base surface). In addition, the summit of the edifice is characterized by a system of coalesced craters: in its north-eastern sector, they were identified by fragments of two concentrically located craters of 300 m and 270 m in diameter, whereas the central part was occupied by one main crater, 200 × 280 m, that was formed during the 09.05.2006 eruption. The position of this crater, most strikingly, remained relatively stable from this point onwards ( Supplementary Fig. 4).
The 2013 overflight images show that all flanks of the former dome complex were covered by lava flows and pyroclastic deposits ( Fig. 2h; Fig. 5c, d). In particular, the lava flows from the 2009 to 2012 eruptions covered the southern and south-eastern slopes. Pyroclastic deposits add another 30-55 m in thickness in the western sector of the amphitheater. As the edifice acquired a rather symmetric and conical shape owing to repeated coverage by lava and pyroclastics, we can conclude that by 2013, Bezymianny's former dome complex eventually had converted into a typical stratocone 45 , which hosts a summit crater and has flanks made up of alternating effusive and explosive deposits.
At last, the data set from 2017 revealed that new lava flows filled the summit crater, and poured onto the western and southwestern slopes during the 2016-2017 eruptions 46   new stratocone started to merge with the flanks of the precollapse edifice. Moreover, the height of the cone reached 3020 m (~820 m in relative height), which is~90 m lower than the precollapse height in 1949.
Over the entire period of Bezymianny's regrowth, the total volume of the rebuilt central cone amounts to 0.591 km 3 ( Supplementary Fig. 5), which suggests that the regrowth of Bezymianny between 1956 and 2017 occurred at an average rate of 26,400 m 3 /day. The missing volume between the current and the pre-collapse edifice is~0.147 km 3 . Thus, assuming that activity at Bezymianny will continue with the detected average rate, the complete refilling of the amphitheater up to the precollapse height and volume may be reached between 2030 and 2035.
Vent migration and focusing. After the 1956 collapse, the first endogenous dome (Fig. 2a) emerged from a vent with a location shifted~500 m to the ESE of the pre-collapse summit ( Fig. 6; Supplementary Fig. 4). The inferred vent of the second dome emerged 400 m north-westward from the vent of the first dome, and~230 m north-eastward from the location of the pre-collapse vent. In the period 1967-1976, the transition to exogenous dome growth was associated with vents concentrating, now having a maximum distance (D max ) between two furthest vents of~210 m. During 1977-2006, which was dominated by lava flow activity, the vents focused further, now with a D max of~150 m. Since 2006, when the lava dome started to convert into a stratocone, D max has dropped further to~80 m. The current eruption center is located 270 m ESE of the pre-collapse vent location. Therefore, our observations provide evidence of the gradual vent centralization at the regrowing volcanic edifice. Although a few hundred meters distances between the vent locations along the northwest-southeast axis of the amphitheater are observed (Fig. 6c), lesser variation is visible perpendicular to the sector collapse direction (Fig. 6b).
Based on studies analyzing the effects of loading and unloading at volcanoes and their magma pathways 14,47,48 , we designed numerical models to investigate (i) how the 1956 sector collapse may explain the observed shift in vent location, and (ii) how regrowth of the central edifice relates to the focusing of the vents. We use a 2D boundary element model to simulate magma propagation paths. Our model refers to the cross-section of a mixed-mode magma-filled crack in plain-strain approximation. The most favorable direction for the magma pathway is computed iteratively, following the criterion of maximum energy release 14,47 . This approach is based on fracture mechanics principles, however, it neglects the dynamics of magma within the fracture so that it cannot account for magma viscosity 49 . Further detail on the modeling technique and assumptions are provided in the Method section.
Our model domain is a vertical~WNW-ESE cross-section of Bezymianny volcano (Fig. 6c) and extends from a depth (z) of −0.1 km up to 1.9 km above sea level, which is the approximate altitude of the post-collapse amphitheater's floor (dashed profile in Fig. 6c). We simulate the stress field owing to topographic changes using loading and unloading force distributions applied at the top of the model 14,47,48 (z = 1.9 km). These force distributions have simplified shapes based on the profiles in Fig. 6c (triangles and trapezoids in Fig. 7). Magma paths start vertically, from the base of the model domain (z = −0.1 km), within a 3 km wide region centered below the summit of the precollapse edifice (−1.5 < x < 1.5 km). We choose the depth and extension of the region from where the intrusions start in order to cover the whole area affected by the stress changes associated with the morphological changes we study. This area should not be considered as the location of a magmatic reservoir (which should be deeper at Bezymianny, c.f. Discussion section). Rather, our assumption is that an intrusion would rise from deeper depths and enter vertically our model domain. These models, although being highly simplified in terms of topographic geometry and magma dynamics, are constrained by the Bezymianny morphologies as summarized in Fig. 6c, and the resulting magma paths can be compared with observations of post-collapse vent-shift and the final vent focusing. In fact, magma pathways are influenced by the local stress, whereby changes in the direction of principal stress affect the favorite direction of propagation (see Fig. 7; Supplementary Fig. 6). We find that the unloading owing to the 1956 sector collapse favors magma rising towards and beneath the collapse amphitheater~300-800 m east of the precollapse vent location (0.28 < x < 0.83 km, Fig. 7a), which is in agreement with our photogrammetric observations at Bezymianny and with previous results for the shift of vent location as a consequence of a flank collapse 14 ( Supplementary Fig. 7). Our simulations for the post-collapse stress scenario results in magma paths that are rather scattered beneath the collapse embayment (highlighted in red in Fig. 7a). In contrast, the paths beneath the rim of the collapse scarp are more focused-pointing at the location of the rim-and get arrested beneath the rim, because of higher compressive horizontal stress due to the rim topography and the horizontalization of some paths, which reduce the effective buoyancy of the intrusion (black paths in Fig. 7a). Also, our photogrammetric observations at Bezymianny suggest that initial dome growth was spread out, causing the formation of two main coalescing domes elongated along the collapse amphitheater. When considering a stress scenario that accounts for the regrowth of the central cone (Fig. 7b), the number of magma pathways rising toward the collapse amphitheater increases, and paths focus, pointing at the location of the center of the new load (Fig. 7b). The distribution of magma pathways beneath the regrowing edifice resembles the distribution of paths below the scarp rim, however, the compressive stress beneath the growing cone is still lower (as in this simulation the new dome topography did not reach the height of the scarp, yet). In addition, a zone of lower compressive stress persists between the two lobes of maximum horizontal compression (in 0 < x < 0.3 km). Such a "corridor" of lower compressive stress may further help magma rising within the new cone (dark red path in Fig. 7b). This is also   The colored contours represent the horizontal stress owing to the unloading and reloading force distributions that simulate the effect of edifice collapse and regrowth, respectively (gray triangles and trapezoid shades, respectively). Gray dashes indicate the direction of maximum compression. The trajectories and the shape of magmatic intrusions are plotted in red when they are deflected toward the collapse embayment. Magma pathways that converge toward the rim of the collapse scarp are marked by solid black lines. a Post-collapse scenario: the pre-collapse vent location (white reverse triangle), appears to be inhibited (paths diverge to both sides away from it). The most favorable position for post-collapse vents is within the collapse embayment (red trajectories), because of collapse unloading forces. New vents may be shifted up to several hundreds of meters towards the collapse direction. b Regrowth scenario: more paths are deflected toward the collapse embayment, and they tend to converge toward the center of the loading force distribution associated with the new dome (gray trapezoid). Even though the compressive stress below the collapse embayment is increased by the new dome, it is still lower than the horizontal compression underneath the rim. A region of relatively lower compressive stress (paths highlighted in dark red), may represent a more favorable location for magma paths to intrude within the new dome. in agreement with our observations, with vents progressively focusing towards the center of the new cone.

Discussion
Repeated decadal acquisition of photogrammetric data allows us to distinguish stages of the Bezymianny regrowth ( Fig. 8; Supplementary Table 2). The stratocone edifice was destroyed by the 1956 sector collapse, decapitating the central conduit, and triggering a catastrophic lateral explosion. These events left a pronounced amphitheater morphology. New magmatic activity resumed inside this amphitheater further eastward of the former conduit zone. The post-collapse activity was initially in the form of endogenous growth from two main eruptive centers but changed to primarily exogenous growth with shear lobe extrusions from migrating vents. Then, effusive activity manifested with lava flows emplacement and has continued to the present. As the morphology of the dome built up, it gradually developed into a symmetric stratocone with alternating lava and pyroclastics deposition (see Supplementary Movie 1).
The results from numerical modeling show that the stress change associated with edifice destruction and rebuilding causes the shift of new eruptive vents towards the collapse embayment, (as previously shown for other volcanoes 14 ), and magma pathways pointing toward the center of the area where new topography has formed. These results are consistent with the observation of scattered endogenous dome growth after the collapse, and the successive dome-to-cone transition. However, our numerical models necessarily simplify a much more complex system: several aspects may affect the propagation paths of magma, and their likelihood to feed eruptions. Among others, the interaction with rock heterogeneities, as well as changes in magma density and viscosity (owing to degassing, for instance), were not considered in our simulations. Although it is known that changes in the volcanic edifice loading may affect magma composition and eruptive style 10 , we still know very little about the effect of magma viscosity on the propagation paths of magmatic intrusions 50 . The effect of these processes should be investigated with dedicated modeling approaches in order to achieve better insights on their relevance for the observations considered here. Nevertheless, a shift of post-collapse vent locations has been observed at several volcanoes, and Bezymianny fits well within the trend previously shown in ref. 14 (Supplementary Fig. 7), suggesting a common mechanism-related to the amount of mass redistribution-which may be primarily due to the local stress state induced by large topographic changes.
However, mechanisms of regrowth after a collapse likely depend also on factors other than stress changes, such as eruption rate variations (Table 1, Fig. 9, Supplementary Fig. 3). During Bezymianny's edifice reconstruction, the growth rate decreased noticeably from an initial~56,000 m 3 /day (1956-1967) to 33,000 m 3 /day (1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977) and then to 15,500-17,000 m 3 /day . The constant growth rates coincided with the occurrence of lava flow emplacements that successively covered the older dome complex. We conjecture that the lava flows may have acted as an armor preventing partial flank collapses and stabilizing the new edifice, which eventually facilitated the formation of a stable centralized vent through gradual load increase. Thick and unevenly distributed shear lobes, in turn, may contribute to slope instability leading to partial dome collapses, as was observed at Bezymianny in 1985 51 and frequently at Shiveluch 52 .
Similar to Bezymianny, the post-collapse dome growth at Shiveluch during the 1980-1981 extrusive period was characterized by an initially higher rate of 186,000 m 3 /day that later decreased to~60,000 m 3 /day 19 . Over the 2001-2012 period of Shiveluch regrowth, the initial rate of resumed extrusive activity was 700,000 m 3 /day that gradually decreased to 150,000 m 3 / day 52 . After the collapse of Santa Maria volcano, regrowth of Santiaguito dome complex was also characterized by multiple episodes 53 , the highest rate 178,000 m 3 /day was observed during the first post-collapse episode in 1922-1925. Further the growth of each new dome began with a 3-5-year period of higher rate (50,000-130,000 m 3 /day), followed by a decrease to Helens), the prevalence of destructive over constructive processes (Shiveluch), and/or because these volcanoes have not produced any lower viscosity lava flows that covered and thus stabilized the dome. For Bezymianny, we argue that the transition from the extrusion of shear lobes to the effusion of lava flows in 1977-2006 may mirror a contemporaneous decrease of the viscosity of its eruptive products. In fact, petrologic studies 37,56 showed that the silica content decreased constantly between 1956 and 2012 from 60.4% to 56.8% but without strong kinks (Fig. 9b). This suggests that the rather rapid change in the growth style and extrusion rates was accompanied by a smooth decrease in the silica content.
Moreover, in ref. 37 , it was pointed out that the gradual compositional changes may be related to a change of the source depth. The magma feeding system at Bezymianny includes at least two reservoirs 37,57 . The shallower reservoir is located at~7 km depth, and the deeper reservoir at~18 km; they provide possible endmember mixing possibilities explaining gradual compositional and morphologic transitions. In the concept 37 , the beginning of the dome growth in 1956 was fed by the shallow magma reservoir, which was successively replenished and fed by the deeper reservoir. Since the 1977 eruption, when Bezymianny produced the first lava flow, the deeper magma source became dominant.
The hypothesis of the activation of two different feeding sources, cannot be supported-neither is in contrast-with our magma path simulations: in fact, our model cannot account for viscosity changes. Also, accounting for the mechanical coupling with a magmatic reservoir (not considered here) may not add further constraints, as the unloading and loading force distributions affect magma pathways only at shallower depths (because of the relatively small collapse volume with respect to other cases previously studied with this approach 14 (Supplementary Fig. 7).
Nevertheless, an effect of unloading and loading stress on the activation (or inhibition) of reservoirs at depth, may be possible 10 . At the moment, it is not clear for how long the effect of topographic loading may last, but given the regrowing load of the new cone is approaching pre-collapse conditions, we conjecture that also the stress relationship at depth may become close to the condition as it was before the 1956 sector collapse.
In the past decades, near-continuous activity and recurrent emplacement of lava flows have almost fully restored Bezymianny's conical pre-collapse shape. This example shows that the process of post-collapse regrowth can be relatively fast, possibly displaying an evolution of volcanoes in general. The regrowth at Bezymianny indeed compares with other similar type andesitic volcanoes, although time scales may significantly vary.
The Santiaguito dome complex has been formed since 1922, thus started two decades after the Santa Maria volcano collapse in 1902. During the early stages of volcanism, endogenous and exogenous lava domes formed at Santiaguito, gradually changing to lava flow dominated regime today 53 . Santiaguito also experienced the distribution of eruptive centers with a D max of~1300 m and subsequent centralization at the (currently most active) Caliente vent 17 . On Montserrat Island, Soufriere Hills volcano regained its conical shape in two decades after the 1997 collapse 58 , but the scale of the collapse and regrowth is much smaller than at Bezymianny. On Anak Krakatau (Indonesia), a stratocone was rebuilding after the 1883 catastrophic eruption, and had repeated collapses to the SW, most recently in December 2018 4 . For the majority of other basaltic-andesite and andesitic volcanoes, the post-collapse regrowth takes thousands of years. For example, cone regrowth after sector collapse at Avachinsky volcano (Russia) took 4000-5000 years 59 , at Colima volcano (Mexico) 4000 years 60 , at Merapi volcano (Indonesia)~1900 years 61 , and at Parinacota volcano~2000 years 18 .
As Bezymianny almost gained the pre-1956 collapse height, a new catastrophic collapse might be prone 21,26 . During the Holocene, Bezymianny repeatedly produced large-volume pyroclastic flows up to 0.3-0.4 km 3 39 that may indicate the occurrence of large destructive events. However, each of those powerful eruptions were followed by a long period of effusive and moderate explosive activity that lasted from hundreds to thousands of years 39 . The last destructive episode before 1956, was associated with partial flank collapses that left behind two noticeable scars on the western and eastern flanks 34 , yet their extent is not even close to that of the 1956 (Fig. 1c, e, Supplementary Fig. 2). Moreover, as the 1956 collapse was preceded by a long period of quiescence 39 and was triggered by intrusion of a cryptodome (i.e., shallow magma body destabilizing the strength of the edifice) 32 , we cannot say whether the next catastrophic edifice collapse will occur shortly after reaching the pre-collapse size. What we can say, is that if the current activity will go on (the latest two eruptions occurred in 2019), causing further steepening of the flanks, the volcano's upper slopes may eventually partially collapse, even in the near future. This clearly poses a hazard to the visitors of Bezymianny. Furthermore, the area of the potential hazard has increased since 2017, when the northern and southern sectors of the amphitheater's rim were covered with pyroclastics, which led to an unhindered inflow of eruptive material toward the northern and southern feet of the volcano.  Table 1

Methods
Aerial photogrammetry. Aerial photography of the volcano was first performed in 1949 (before the 1956 collapse). The post-collapse topography was estimated in previous research 34 from the 1956 ground-based images, which were taken when the first dome occupied a small area of the amphitheater's floor (Fig. 1e). The base of the amphitheater was estimated to be~2200 m above sea level, and its vertical sections across the collapse direction appeared to be close to parabolic. The amphitheater inner walls in the 1967 elevation model were extrapolated according to the estimated parameters 34 . Information about the morphology after the 1956 sector collapse was also derived from ground-based images 20,41 , and further by aerial surveys, which were regularly realized since 1967. We used the 1949 archive images to extract the pre-collapse DEM (Supplementary Fig. 2). Then we focused this study on the volcano regrowth by processing eight sets of near nadir aerial photographs acquired in 1967, 1968, 1976, 1977, 1982, 1994, 2006, and 2013. The photographs were taken with dedicated topographic analog cameras AFA 41-10 (1967, 1968, 1976, and 1977; focal length = 99.086 mm), TAFA 10 (1982 and1994; focal length = 99.120 mm), and AFA TE-140 (2006 and 2013; focal length = 139.536 mm). The cameras have an 18 × 18 cm frame size. The acquisition flight altitude above the mean surface of Bezymianny varied from 1500 to 2500 m in order to capture the entire 1956 collapse amphitheater and surrounding slopes during one pass. For photogrammetric processing, we used 3-4 consecutive shots that provided a 60-70% forward overlap.
The analog photo negatives were digitized by scanning with a resolution of 2400 pixels/inch (approx. pixel (px) size = 0.01 mm). The mean scale within a single photograph depends on the distance to the surface and corresponds on average to 1:10,000-1:20,000. Thus, each px in the scanned image represents~10-20 cm resolution on the ground.
The software packages Erdas Imagine 2015 v15.1 62 and Photomod 5 63 were used for processing, allowing to perform interior and relative orientation of the aerial photographs 64 , exterior orientation of the stereo models, triangulation, and DEM extraction and correction.
For the interior orientation, the analog cameras' focal length, frame size, lens distortion, and the position of the main point and fiducial marks were included. The relative orientation of adjacent photographs was performed automatically based on 25 tie points (root mean square errors (RMSE) = 0.1 pixels).
The coordinates of 12 ground control points (GCPs) were derived from a Theo 010B theodolite data set collected at geodetic benchmarks during a 1977 fieldwork. These benchmarks were established on the slopes of Bezymianny before the 1977 aerial survey and then captured with the AFA 41-10 aerial camera. RMSE of the GCPs did not exceed 0.06 m. From this data, we created the 1977 stereo model referenced to the USSR State Geodetic Network 65 , which itself served as a reference for the following acquisitions similar to how it was performed for Mt. St. Helens 31 .
Those benchmarks covered by volcanic deposits made it necessary to extract new coordinates similar to ref. 31 of six clearly distinguishable and stable (not affected by eruptions) topographic prominences in the vicinity of Bezymianny (peaks or large rocks). These locations were identified in the georeferenced 1977 stereo model and defined as GCPs used for the exterior orientation of the previous (1949, 1967, 1968, and 1976) and succeeding (1982,1994,2006, and 2013) stereo models. Thus, all other stereo models were oriented according to the 1977 model that allowed to perform triangulation for each of them with RMSEs varying from 0.4 to 1.9 m (Supplementary Table 1) depending on the age of images and snow coverage.
The oriented stereo models were used to automatically extract points in the Erdas Enhanced Automatic Terrain Extraction (eATE) module, which applies a normalized cross-correlation algorithm with a window size set to 11 × 11 px, and with a correlation range of 0.2-0.7 for the highest pyramid level and the last pyramid, respectively. The obtained point clouds in LAS format were filtered in CloudCompare v2.9.1 (https://www.danielgm.net/cc/), with the noise filter tool (spherical radius = 0.75 neighbors), which resulted in an average number of points per point cloud of~300,000 for 5 km 2 (0.06 points/m 2 ).
As intensive fumarolic activity prohibited automatic detection of the ground surface, some areas were not processed properly producing gaps in the point clouds. To solve this problem, we imported each point cloud in the corresponding stereo models using the Photomod DTM module and then performed further manual extraction of the missing points by placing a floating mark on the surface in anaglyph stereo mode and storing XYZ coordinates. This manual approach allowed us to visually identify the surface through the light steam that is not possible for the automatic algorithms, similar to performing in ref. 31 . This allows also a validation of the automatically extracted points and to process some of the dome elements with better resolution. The resolution of final point clouds varies from 2 m to 30 m depending on the complexity of the surface.
Photographs from 2002,2005,2009, and 2010 aerial surveys were considered for visual interpretation and identification of recent lava flows.
Satellite and unmanned aerial vehicle photogrammetry. To complement the temporal coverage of our aerial time-series, we additionally tasked tri-stereo highresolution Pléiades satellite imagery acquired on 09.09.2017. We used Erdas Imagine to process the three overlapping monochromatic images (spatial resolution = 1 m) in a manner similar to the processing chain described in ref. 66 . We identified 45 automatically and manually tracked tie points and used them for the relative orientation of the images (RMSE = 0.1 pixels). The exterior orientation was calculated automatically by employing the provided Rational Polynomial Coefficient data. Eventually, we extracted DEM with the Erdas eATE module and subsequently filtered it with the CloudCompare noise filter. The final LAS point cloud has~10 million points for 60 km 2 (0.16 points/m 2 ).
Yet, the Pléiades images captured strong degassing and shadows, which caused a gap in the point cloud along the western dome flank. But we were fortunate to gather optical images of this area with an Unmanned Aerial Vehicle (UAV) (DJI Mavic Pro) in July 2017. We processed the UAV data set with Agisoft Metashape v1.6.4 (https://www.agisoft.com/downloads/installer), which is carried out by using the default settings for photo alignment (high quality with a 40,000 key-point and 1000 tie-point limit), dense cloud generation (high quality, aggressive depth filtering), and DEM building (WGS84) based on the dense cloud. Eventually, we merged both the Pléiades and the UAV derived point clouds with CloudCompare (RMSE = 0.8 points) and subsampled it to a spatial resolution of 2 m.
Data alignment and referencing. The investigated aerial, satellite and drone data have the same spatial scale but shifted in geo-position, depending on the coordinate system (USSR State Geodetic Network for the aerial data and WGS84 for the merged Pleiades-UAV data) and precision. In order to allow quantitative comparison and analysis of the differently located data, we aligned the aerial point clouds relatively to the merged Pléiades-UAV point cloud with eight distinct points on the 1956 amphitheater's rim using the CloudCompare software. To avoid alignment errors, we first aligned the 2013 point cloud to the 2017 point cloud (RMSE = 1.3), so that the effect of degradation of the amphitheater rim is minimal. The remaining aerial point clouds were then successively referenced to the next younger point cloud (RMSE = 1.3-1.7). Finally, we obtained nine LAS point clouds referenced to the WGS84 (UTM zone 57 N).
Volumes and growth rates estimation. To quantify changes throughout the evolution of Bezymianny's new edifice, we calculated volumetric differences between consecutive DEMs with the 2.5D volume calculation tool in CloudCompare, which allowed us to automatically estimate the added and removed volumes by comparing two surfaces. As material may be added (e.g., intrusions, lava, pyroclastics, and tephra deposition) and/or removed (e.g., erosion, explosive excavation, small flank collapses deposited outside the area of the dome) during growth, volume differences may be represented by positive and negative changes that yield in a net volume change. The positive changes were then used to determine growth rates by dividing the volume by the time interval between the dates of the survey.
The volume of the dome over the 1956-1967 period was calculated in previous research 34 . The 1967-2017 volumes were estimated only within the areas affected by the eruption (dome areas) ( Table 1). CloudCompare provides a possibility to preselect an area for volumetric comparison. Each time comparing consecutive DEMs, we outlined a clearly visible border between the dome's feet and the amphitheater's inner walls on both DEMs relatively to a later one.
To avoid errors in volume calculation caused by the point clouds alignment, volumes for the aerial DEMs were derived from initial point clouds (unaligned to the 2017 Pleiades-UAV point cloud). Only for the 2017 volume estimation, we used the base surface from 2013 aligned to the 2017 point cloud.
As the exact acquisition dates for the 1967 and 1968 aerial data are unknown, we estimated the acquisition dates based on the angle of incidence of the sun. For each stereo model, we identified the respective pairs of points (peaks on the amphitheater's rim and shadows from these peaks on the slopes of the dome). The angle data sets from each stereo model were averaged; the azimuths were 209.1°a nd 239.8°, and the angles of the sun's fall relative to zenith were 71.3°and 65.3°for the 1967 and 1968, respectively. Using NREL's Solar Position Algorithm 67 , we made an iterative search for the date options correspond to these values. For 1967, the appropriate options were 19 Feb. and 24 Oct.; for 1968, they were 3 Apr. and 9 Sept. The estimation accuracy was ±1 day. As both 1967 and 1968 data show only partial snow coverage, the dates of 19 Feb. and 3 Apr. were rejected as too early. The region has snowy winters with a continuous snow cover until late May. For the dates of 24 Oct. and 9 Sept., the snow situation is quite possible, as the snow starts to fall here in the mid-autumn.
Error estimation. The volume estimation errors (Table 1, Supplementary Table 1) in this study mainly depend on the triangulation errors (TRMSE), which are calculated automatically in Erdas Imagine and presented in the triangulation reports. Possible sources of the TRMSE can be camera lens optical distortion, deformation of old analog films, scanner distortion, and manual GCPs assignment. Furthermore, the resolution of the extracted point clouds can also contribute to the volume error. Finally, the RMSE of the point clouds alignment (ARMSE) affects the volume estimation accuracy.
As Bezymianny's dome is relatively circular symmetric, a slight XY shift of a DEM does not contribute to the volume error as much as Z shift, since a volume adding at one side of the dome is compensated by the same amount of volume removed from the other side. Thus, we consider only Z TRMSE. To determine the contribution of the TRMSRE in volume estimation, we distributed each TRMSE Z value over the respective area affected by eruption (dome area) ( Table 1) similar to ref. 31 . Depending on the TRMSEs, volume uncertainty varies from 0.5 million m 3 (in 2017) to 3.2 million m 3 (in 1967) that is 0.1-1.4% of the dome volume.
To estimate volume uncertainty caused by point clouds resolution we used equation (14) from ref. 68 . The base parameters for the calculation are the mean distance between points in each cloud, which varies from 2 m (in 2017) to 7 m (in 1977 and 2006), and the standard deviation of the error in the spatially distributed points, which varies from 0.2 m (in 2017) to 0.5 m (in 1967 and 1968). The resolution uncertainties in volume increase between two consecutive point clouds vary from 3900 m 3 (in 2017) to 7900 m 3 (in 1982) that is >0.001% from the dome volume and might not be taken into account.
As the 2017 DEM was combined from the two point clouds (Pleiades and UAV) with the 0.8 m RMSE in alignment, we distributed this error over the area of the aligned UAV point cloud (430,000 m 2 ). The contribution of this ARMSE to the 2017 volume estimation is 344,000 m 3 (0.06%). The 2013-2017 point clouds alignment RMSE (1.3 m) distributed over the 2017 cone area contributes 2,900,000 m 3 (0.5%) to the volume error. Thus, the total ARMSE contribution in the 2017 cone volume is 0.6%, and the total 2017 volume error is 0.7% (Supplementary Table 1).
The errors of the growth rates were calculated by dividing the root sum squared volume errors for neighboring dates by the time between them. This model of inaccuracy 69 was used because errors at different dates have an equal effect and are independent of each other. The smallest errors of the rates (2.1-4.3%) turned out to be for the periods between 1968-1976, 1982-1994, 1994-2006 because of their longer length and, thus, a high averaging of the values. The largest relative rate error (31.2%) occurs for the interval 1967-1968 owing to the small increase in the volume, which is almost equal to the absolute values of the volume estimation error.
Morphological mapping. To map the features of the Bezymianny new edifice we performed visual interpretation of the high-resolution stereoscopic aerial images in anaglyph mode using the Stereo Analyst module of Erdas Imagine. This method allows an operator, wearing anaglyph glasses, to see the detailed morphology of the studied object in 3D and perform all measurements of the identified features. Each stereoscopic image was carefully compared to the next one to trace the same features and to analyze their development or decay. The results of the interpretation of the dome's features are based on commonly accepted classification of volcanic landforms and their elements 43,70 , which reveals the nature and sequence of their formation. The exogenous growth was distinguished from the former endogenous growth by a presence of separate shear lobes, which extruded from the vents or open craters at the different sectors of the dome. The lava flows were recognized by their low thicknesses (from 15 to 30 m), and length to width ratios (from 1.8 to 2.1) in comparison with the thick (from 60 to 110 m) and almost circle-symmetrical shear lobes. Morphologically lava flows differ from shear lobes having a drop shape; their near-vent parts are narrower and thinner than fronts, whereas shear lobes have the highest thickness at their top and near-vent parts as they become thinner at their margins owing to the intensive crumbling of solidified material. There are also differences in the directions of deposition: lava flows spread toward the foot of the dome, and shear lobes are accumulated at the place of extrusion. The lava plugs were identified as cylindrical bodies extruded at the vent areas of the previously emplaced lava flows.
To estimate the migration of the vents in the early decades of edifice regrowth (1956-1967) we consider a vent as a center of an endogenous lava dome, which corresponds to ref. 70  The changing location and migration of the vents are visualized in profiles (Fig. 6) and along sections created in directions NE-SW and NW-SE across and along the collapse amphitheater. We project the location of the pre-collapse central vent, and of the post-1956 collapse vents from each period of regrowth onto those sections. This allows us to indicate the location of the two most expressed distal vents for each period of regrowth in corresponding vertical lines with remaining maximum distances (D max ) from each other according to the corresponding profile (Fig. 6d, e).
Magma pathway simulations. Numerical simulations were performed following a previously published approach 14,47,48 . A short description of the main characteristics and assumptions of this approach are given here.
Magma pathways are modeled as 2D boundary element mixed-mode cracks in plane strain approximation. The modeling plane is perpendicular to the crack plane, therefore the model refers to a cross-section of the intrusion. The crack opening depends on assigned normal and shear stress boundary conditions that are the magma overpressure and the shear component of the topographic stress, respectively. The overpressure along the crack is given by the difference between the magma pressure and the confining stress (superposition of the lithostatic pressure and the normal component of the topographic stress). Topographic stresses are computed using analytical formulas for loading and unloading forces at the surface of an elastic half-space 14,47,48 . Loading and unloading force distributions are used to model the loading of the pre-collapse volcanic edifice and the stress change induced by the flank collapse 14 . The model for the magmatic intrusion accounts for the magma buoyancy and compressibility, whereas it neglects magma viscosity. As a consequence, the magma overpressure profile within the crack is linear along the depth and proportional to the magma-rock density difference (hydrostatic overpressure profile). The trajectories within the crust are determined by testing the incremental elongation of the magma-filled crack in different directions. Our algorithm chooses the direction where the sum of elastic and gravitational energy release is maximal 47,48 .
Here, we applied this model to three stress scenarios, pre-collapse, post-collapse, and cone regrowth.
Pre-collapse scenario: we consider the loading due to the pre-collapse volcanic edifice based on the pre-collapse profile (black solid line) in Fig. 6c. The effect of topography on the crust beneath is modeled by vertical forces applied on a reference horizontal surface at 1900 m on the sea level (which is the approximate altitude of the post-collapse embayment, dashed profile in Fig. 6c). This reference surface represents the upper boundary for the stress calculation and for the magma paths in all simulations. The force distribution is triangular, centered in x = 0 (the center of the precollapse edifice), with height H prc = 1.2 km (so that the altitude of the pre-collapse volcanic edifice is 3100 m), and with base W prc = 4.4 km. The effect of the triangular loading force distribution is computed in the rocks beneath 14 . The magnitude of the loading force is proportional to shallow rock density (r1 = 2000 kg/m 3 , times the gravity acceleration) so that the loading forces represent the mass of the rocks standing over the reference surface. Free surface effects are neglected. We assume that during the emplacement of the pre-collapse edifice, non-elastic effects may have partially released elastic stress, this is accounted by testing different "Effective Volcanic Loading" (EVL) 14 . Conversely, the flank collapse and the dome growth are very recent events, and therefore we consider them as "purely elastic" 14 . Different values of EVL were tested (0.5-0.6-0.7), meaning that 50%, 40%, and 30% of the elastic loading had been released at the time of the collapse 14 . We used the triangular difference approach to simulate the stress relaxation: "Trapezoidal Loading" in ref. 14 . Finally, we obtain three scenarios for the "Pre-collapse", one for each of the EVL values ( Supplementary Fig. 6).
Post-collapse scenario: we introduce the unloading force distribution simulating the flank collapse (superposed on the stress owing to the pre-collapse edifice). Here we used a simple triangular profile to approximate the effect of the collapse. The maximum height of the unloading is H poc = 0.8 km (in x = 0), the base of the unloading is W poc = 2.8 km, from x = −0.6 km to x = 2.2 km. These geometrical constraints are based on the difference between the pre-collapse and post-collapse profiles in Fig. 6c (black solid and dashed lines, respectively). Similar to the precollapse scenario, the magnitude of the unloading force is proportional to the density of shallow rocks, but with opposite sign. We do not consider any stress relaxation effect for the unloading forces 48 , as the collapse event is sudden and recent.
Regrowth scenario: the reloading due to dome growth is introduced. We used a trapezoidal shape for the dome. Based on the 1956-1967 profile of Fig. 6c (first phase of dome growth), the trapezoid has a height H dg = 0.5 km, with a lower-base W dg = 1.1 km (from x = −0.25 km to x = 0.85 km) and an upper-base W dg = 0.37 km. Again, forces are proportional to the shallow rock density, and they are superposed to the post-collapse scenario.
Magma paths in interaction with the background stress are computed for each scenario. Magma paths start from 2 km below the upper boundary for stress calculation (z = 1900 m). For each scenario, we computed 19 independent magma path trajectories with paths starting vertically oriented from a depth z = −0.1 km, on an area as wide as 3 km, centered in x = 0 (center of the pre-collapse edifice).
Magma and rock parameters are as follows:

Code availability
The Fortran90 code used for the magma paths simulations and the instructions about how to compile the code and run the simulations are available via the Zenodo repository at https://doi.org/10.5281/zenodo.3957577.