Diego Coppola

Laiolo Marco

Francesco Massimetti

Sebastian Hainzl

Alina V Shevchenko

René Mania

Nikolai M Shapiro

Thomas R Walter

Rene Mania

Thermal remote sensing reveals communication between volcanoes of the Klyuchevskoy Volcanic Group

come   L'archive ouverte pluridisciplinaire

Introduction

Closely located or clustered volcanoes may become conjointly active and are hence considered especially hazardous, yet robust evidence for their connectivity remains sparse. Examples of such a synchronized volcanic activity are discussed for neighboring volcanoes in Iceland 1 , Alaska 2 , Kamchatka [START_REF] Eichelberger | Eruption of andesite triggered by dyke injection: contrasting cases at Karymsky Volcano, Kamchatka and Mt Katmai, Alaska[END_REF][START_REF] Gusev | Great explosive eruptions on Kamchatka during the last 10,000 years: Self-similar irregularity of the output of volcanic products[END_REF] , Italy [START_REF] Walter | Possible coupling of Campi Flegrei and Vesuvius as revealed by InSAR time series, correlation analysis and time dependent modeling[END_REF] , and elsewhere [START_REF] Miklius | Interaction between Kilauea and Mauna Loa[END_REF][START_REF] Biggs | The lateral extent of volcanic interactions during unrest and eruption[END_REF] , although larger time-scale synchronicity has been also reported for global volcanism [START_REF] Kutterolf | A detection of MILANKOVITCH frequencies in global volcanic activity[END_REF] . Reasons for the linked activity of adjacent volcanoes are only poorly understood and may be locally different, including triggering by large tectonic earthquakes and associated stress changes within the crust [START_REF] Nostro | Two-way coupling between Vesuvius eruptions and southern Apennine earthquakes, Italy, by elastic stress transfer[END_REF][START_REF] Hill | Earthquake-Volcano Interactions[END_REF][START_REF] Walter | The influence of great earthquakes on volcanic eruption rate along the Chilean subduction zone[END_REF][12][START_REF] Kennedy | What effects do earthquakes have on volcanoes?[END_REF][START_REF] Seropian | A review framework of how earthquakes trigger volcanic eruptions[END_REF] , and the competition of volcanoes for common reservoirs [START_REF] Eichelberger | Eruption of andesite triggered by dyke injection: contrasting cases at Karymsky Volcano, Kamchatka and Mt Katmai, Alaska[END_REF][START_REF] Dzurisin | Magma supply to Kīlauea Volcano, Hawai'i, from inception to now: Historical perspective, current state of knowledge, and future challenges[END_REF] . Conjoint unrest and deformation activity at clustered volcanoes occurs with temporal delays of days to months (or even more) and appears to be distancedependent [START_REF] Biggs | The lateral extent of volcanic interactions during unrest and eruption[END_REF] . Magmatic sources spaced less than about 10 km apart tend to interact, whereas those spaced over 25 km do not [START_REF] Biggs | The lateral extent of volcanic interactions during unrest and eruption[END_REF] . However, interactions over longer distances (>20 km) have been hypothesized for volcanoes that share a common deep source, or in response to large dike intrusions or subduction earthquakes [START_REF] Biggs | The lateral extent of volcanic interactions during unrest and eruption[END_REF] . While observations suggest positively correlated feedback (one volcano triggers the other one), only a few examples [START_REF] Walter | Possible coupling of Campi Flegrei and Vesuvius as revealed by InSAR time series, correlation analysis and time dependent modeling[END_REF][START_REF] Gonnermann | Coupling at Mauna Loa and Kīlauea by stress transfer in an asthenospheric melt layer[END_REF] underline the existence of anti-correlated activity (i.e., volcanic unrest may shut down activity in the neighborhood).

Previous studies of correlated and anti-correlated volcanism are mainly based on poor data and reduced or biased eyewitness accounts. This is because reports often describe a singular, possibly sporadic occurrence of conjoint activity change, where eruptions are more likely to be reported than periods of quiescence and decreasing activity. Robust and statistically significant testing of repeat observations was not yet achieved and may overcome some of the reporting limitations.

Here we employed a unique time-series of satellite thermal data (from 2000 to 2020) derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor [START_REF] Coppola | Thermal Remote Sensing for Global Volcano Monitoring: Experiences From the MIROVA System[END_REF] , to study potential volcano interactions within the Klyuchevskoy Volcano Group (KVG) in Kamchatka (Russia). The KVG hosts three adjacent (~10 to 30 km distant) active volcanoes (Klyuchevskoy, Bezymianny, and Tolbachik), with contrasting eruptive products that possibly indicate different magmatic sources and a compound multilevel feeding system [START_REF] Koulakov | Three different types of plumbing system beneath the neighbouring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka[END_REF][START_REF] Kayzar | Deciphering petrogenic processes using Pb isotope ratios from time-series samples at Bezymianny and Klyuchevskoy volcanoes, Central Kamchatka Depression[END_REF] .

Whether this composite feeding system forms a single interconnected trans-crustal magmatic system [START_REF] Cashman | Vertically extensive and unstable magmatic systems: A unified view of igneous processes[END_REF] , and in which degree the three volcanoes interact remains so far an open question [START_REF] Fedotov | The magmatic system of the Klyuchevskaya group of volcanoes inferred from data on its eruptions, earthquakes, deformation, and deep structure[END_REF][START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF] .

The Klyuchevskoy Volcanic Group (KVG) has been monitored for >50 years by the Kamchatkan Branch of the Geophysical Survey of the Russian Academy of Sciences (KBGS). Due to harsh climate conditions, the monitoring system mainly consists of seismic stations and webcams. Nowadays, 17 telemetric stations transmit the data to the central monitoring office, where it is being analyzed in real-time [START_REF] Chebrov | Development of digital seismic network in Kamchatka, in Geofizicheskii monitoring Kamchatki[END_REF] and allowed to investigate the deep structure of KVG [START_REF] Khubunaya | Magma chambers beneath the Klyuchevskoy volcanic group (Kamchatka)[END_REF] . Although localized and dense seismic networks were installed for short-term times [START_REF] Shapiro | Understanding Kamchatka's Extraordinary Volcano Cluster[END_REF][START_REF] Koulakov | Anatomy of the Bezymianny volcano merely before an explosive eruption on 20.12[END_REF] , these experimental networks do not permit decade-long analysis and robust statistics, so that details on eruption occurrence and eruption rates remained hidden. As this work shows, satellite thermal remote sensing reveals complex interactions under the KVG, whereby each volcano is able to influence its neighbors. Distinct eruptive activity appears sometimes in chorus and sometimes separate, comparable to an orchestra. [START_REF] Shevchenko | The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse[END_REF] , and (d) Tolbachik [START_REF] Belousov | Overview of the precursors and dynamics of the 2012-13 basaltic fissure eruption of Tolbachik Volcano, Kamchatka, Russia[END_REF] . All volumes are recalculated for DRE using: (i) a magma density of 2800 kg m -3 and an average density of lava flows of 2500 kg m -3 for Klyuchevskoy and Tolbachik 57 ; (ii) a magma density of 2500 kg m -3 and an average density of lava dome of 2000 kg m -3 for Bezymianny [START_REF] Zharinov | Assessing the Volumes of Material Discharged by Bezymyannyi Volcano during the 1955-2009 Period[END_REF] . A standard error of ±20% (colored fields) takes into account uncertainties in the estimates of bulk volumes and densities of erupted products. Red bold lines within insets, correspond to satellite-derived volumes (this work; see Methods). The grey bars outline the timing of the Bezymianny eruption in 1955-56 (BEZ55), the Great Tolbachik Fissure Eruption in 1975 (GTFE) and the Tolbachik Fissure Eruption in 2012 (TFE). Shaded relief map derived from ArcticDEM digital elevation model (https://www.pgc.umn.edu/data/arcticdem/) and elaborated using QGIS version 3.16.3 (http://qgis.osgeo.org); Time series generated using MATLAB software (www.mathworks.com).

The Klyuchevskoy Volcanic Group (KVG)

The KVG is a prominent volcanic massif located in the northern part of the Central Kamchatka Depression. Dozens of volcanic centers were built during the construction of the massif, which currently has three active volcanoes: Klyuchevskoy, Bezymianny, and Tolbachik (Fig. 1a). The KVG is a very active and relatively young volcanic massif mainly developed during the last 300-400 ka [START_REF] Melekestsev | Volcanism and relief formation[END_REF] . Volcanism is fed by sub-arc mantle, melted under an influx of melts and fluids from the subducting Pacific plate [START_REF] Churikova | Sources and Fluids in the Mantle Wedge below Kamchatka, Evidence from Across-arc Geochemical Variation[END_REF][START_REF] Münker | Behaviour of high field strength elements in subduction zones: constraints from Kamchatka-Aleutian arc lavas[END_REF][START_REF] Green | Magmatic and sedimentary structure beneath the Klyuchevskoy Volcanic Group, Kamchatka, from ambient noise tomography[END_REF][START_REF] Koulakov | Mantle and crustal sources of magmatic activity of Klyuchevskoy and surrounding volcanoes in Kamchatka inferred from earthquake tomography[END_REF][START_REF] Levin | Slab portal beneath the western Aleutians[END_REF] . Additional influx of hot mantle following recent slab detachment [START_REF] Levin | Seismic evidence for catastrophic slab loss beneath Kamchatka[END_REF] , and interaction with metasomatized mantle [START_REF] Dorendorf | Hydrated sub-arc mantle: a source for the Kluchevskoy volcano, Kamchatka/Russia[END_REF][START_REF] Portnyagin | Geochemistry of primitive lavas of the Central Kamchatka Depression: Magma generation at the edge of the Pacific Plate. Volcanism and Subduction[END_REF] contribute to the exceptional level of volcanic activity in the area and the very diverse volcanic manifestations and products.

Seismic activity of the KVG volcanoes is abundant and includes long periods of sustained tremors as well as numerous volcano-tectonic (VT) and long-period (LP) events. The latter mostly occur at two depth ranges: above 5 km and close to 30 km [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF] . Geophysical and petrographic data have been used to infer that all the KVG volcanoes are fed by a common parental magma [START_REF] Ozerov | Petrological-geochemical model for genetic relationships between basaltic and andesitic magmatism of Klyuchevskoi and Bezymyannyi volcanoes, Kamchatka[END_REF] . However, different isotope compositions of rocks from Klyuchevskoy and Bezymianny [START_REF] Kayzar | Deciphering petrogenic processes using Pb isotope ratios from time-series samples at Bezymianny and Klyuchevskoy volcanoes, Central Kamchatka Depression[END_REF][START_REF] Dorendorf | Hydrated sub-arc mantle: a source for the Kluchevskoy volcano, Kamchatka/Russia[END_REF] do not support such view, in favor of multiple magma sources, with only limited interaction. Klyuchevskoy volcano (4750 m a.s.l) is the highest in the group and one of the world's most active volcanoes.

Its recent activity is characterized by the effusion of voluminous basaltic andesite lava flows, often associated with moderate to violent explosive activity. Between 1930 and 2005, the volcano has erupted an estimated ~1.5 × 10 9 m 3 of lava [START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] (dense rock equivalent; DRE), with an average magma output rate (DRE) of ~0.67 m 3 s -1 . The output rate accelerated since 1978, associated with a change in the eruptive pattern that shifted from flank-to summit-dominated eruptions [START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] (Fig. 1b). After nine years of rest, in 2003, the volcano began a new activity phase characterized by nine summit eruptions until December 2019 [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] . Seismic, geodetic, and petrographic data [START_REF] Koulakov | Three different types of plumbing system beneath the neighbouring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka[END_REF][START_REF] Fedotov | The magmatic system of the Klyuchevskaya group of volcanoes inferred from data on its eruptions, earthquakes, deformation, and deep structure[END_REF][START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF][START_REF] Shapiro | Understanding Kamchatka's Extraordinary Volcano Cluster[END_REF][START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] suggest that Klyuchevskoy's eruptions are fed through a sub-vertical, pipe-like conduit extending to a depth of 30-50 km below the volcano, where the primary magma reservoir is located.

On its way to the surface, the magma is stored at a depth of 15-25 km and then transported further upwards to a shallow (3-5 km deep) peripheral reservoir [START_REF] Khubunaya | Magma chambers beneath the Klyuchevskoy volcanic group (Kamchatka)[END_REF] . During ascent, the magma evolves from high-Mg low-Al basalt to low-Mg high-Al basaltic andesite [START_REF] Bergal-Kuvikas | A petrological and geochemical study on time-series samples from Klyuchevskoy volcano, Kamchatka arc[END_REF] , making it different from the eruptive products of the other active KVG volcanoes [START_REF] Ozerov | Petrological-geochemical model for genetic relationships between basaltic and andesitic magmatism of Klyuchevskoi and Bezymyannyi volcanoes, Kamchatka[END_REF] . Eventually, before the conduit reaches the summit crater, numerous radial dikes depart and feed eruptions at the mid-and lower-volcano flanks [START_REF] Fedotov | The magmatic system of the Klyuchevskaya group of volcanoes inferred from data on its eruptions, earthquakes, deformation, and deep structure[END_REF][START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] .

Bezymianny is an andesitic volcano (2,882 m a.s.l.) that reawakened in 1955-56 with a paroxysmal eruption (VEI 5) that disrupted the old cone forming a large horseshoe-shaped crater [START_REF] Gorshkov | Gigantic eruption of the volcano Bezymianny[END_REF] . Since then, near-continuous lava dome growth was accompanied by mostly explosive activity [START_REF] Zharinov | Assessing the Volumes of Material Discharged by Bezymyannyi Volcano during the 1955-2009 Period[END_REF][START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] . The greatest rate of dome growth occurred during the first two decades until 1977 (Fig. 1c), when lava flows were observed for the first time marking a pivotal change in the volcano's dome growth mechanism [START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF][START_REF] Bogoyavlenskaya | Eruption of Bezymianny volcano in 1977[END_REF][START_REF] Shevchenko | The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse[END_REF] . Ever since then, Bezymianny's eruptions showed a recurrent cyclical behavior depicted by extrusive-explosive-effusive activity [START_REF] Zharinov | Assessing the Volumes of Material Discharged by Bezymyannyi Volcano during the 1955-2009 Period[END_REF][START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] . Previous works [START_REF] Carter | Detection of a new summit crater on Bezymianny Volcano lava dome: satellite and field-based thermal data[END_REF][START_REF] Manen | Satellite thermal observations of the Bezymianny lava dome 1993-2008: Precursory activity, large explosions, and dome growth[END_REF][START_REF] Manen | Forecasting large explosions at Bezymianny Volcano using thermal satellite data[END_REF][START_REF] Girina | On precursor of Kamchatkan volcanoes eruptions based on data from satellite monitoring[END_REF][START_REF] Reath | Predicting eruptions from precursory activity using remote sensing data hybridization[END_REF] outlined how this cyclic activity was accompanied by precursory thermal radiation preceding the explosive events by few weeks to days. Until 2017, more than 55 distinct episodes of dome-growth filled most of the 1956 collapse amphitheater [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF][START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] gradually developing a stratocone with an average growth rate [START_REF] Shevchenko | The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse[END_REF] of ~0.30 m 3 s -1 (blue line in Fig. 1c). Geophysical and petrological data suggest a multi-level magma plumbing system beneath Bezymianny volcano with at least three crustal reservoirs located between 10-18 km, 5-8 km, and < 2 km depth [START_REF] Koulakov | Three different types of plumbing system beneath the neighbouring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka[END_REF][START_REF] Koulakov | Anatomy of the Bezymianny volcano merely before an explosive eruption on 20.12[END_REF][START_REF] Ivanov | Magma source beneath the Bezymianny volcano and its interconnection with Klyuchevskoy inferred from local earthquake seismic tomography[END_REF][START_REF] Davydova | Petrology of mafic enclaves in the 2006-2012 eruptive products of Bezymianny Volcano, Kamchatka[END_REF][START_REF] Shcherbakov | Plagioclase zoning as an indicator of magma processes at Bezymianny Volcano, Kamchatka[END_REF][START_REF] Turner | The magma plumbing system of Bezymianny Volcano: Insights from a 54year time series of trace element whole-rock geochemistry and amphibole compositions[END_REF][START_REF] Thelen | Seismic characterization of the fall 2007 eruptive sequence at Bezymianny Volcano, Russia[END_REF] .

The Tolbachik massif comprises two large stratocones, Ostry ("Sharp") Tolbachik (3682 m a.s.l.) and Plosky ("Flat") Tolbachik (3085 m a.s.l.), in the southernmost part of the KVG [START_REF] Churikova | The Tolbachik volcanic massif: A review of the petrology, volcanology and eruption history prior to the 2012-2013 eruption[END_REF] at approximately 30 km distance to Kluchevskoy and ~20 km distance to Bezymianny (Fig. 1a). A 70 km long zone of monogenetic basaltic cones extends across the Plosky Tolbachik cone; whose southern branch was the place of the 1975-1976 Great Tolbachik Fissure Eruption (GTFE) [START_REF] Fedotov | The large fissure eruption in the region of Plosky Tolbachik volcano in Kamchatka, 1975-1976[END_REF] . This eruption produced extensive lava fields composed of highmagnesium and high-aluminum basalts, from northern and southern vents, respectively [START_REF] Churikova | The Tolbachik volcanic massif: A review of the petrology, volcanology and eruption history prior to the 2012-2013 eruption[END_REF] . With a total DRE volume of ~1.5 × 10 9 m 3 (Fig. 1d), it was one of the largest basaltic eruptions in Kamchatka during historical times [START_REF] Fedotov | The large fissure eruption in the region of Plosky Tolbachik volcano in Kamchatka, 1975-1976[END_REF] . After the GTFE, no signs of activity were recorded until November 2012, when increased seismic activity heralded the beginning of a new eruption [START_REF] Belousov | Overview of the precursors and dynamics of the 2012-13 basaltic fissure eruption of Tolbachik Volcano, Kamchatka, Russia[END_REF] . The 2012-2013 eruption took place at the south flank of the Plosky Tolbachik cone and was dominated by Hawaiian-style activity associated with an emplacement of a large lava field [START_REF] Dvigalo | The first quantitative estimates of parameters for the Tolbachik Fissure Eruption of 2012-2013 from aerophotogrammetric observations[END_REF] . During the 205 days of activity, a lava volume of ~0.5 × 10 9 m 3 was erupted, with a gently declining trend throughout the whole eruptive period [START_REF] Dvigalo | The first quantitative estimates of parameters for the Tolbachik Fissure Eruption of 2012-2013 from aerophotogrammetric observations[END_REF] . Satellite geodesy could reveal the intrusion of a 6.1 km long dike intrusion, opening up to 8 m, adding almost 10% to the total eruption volume [START_REF] Lundgren | Dike model for the 2012-2013 Tolbachik eruption constrained by satellite radar interferometry observations[END_REF] .The activity ceased entirely by the end of August 2013. According to Koulakov et al. [START_REF] Koulakov | Three different types of plumbing system beneath the neighbouring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka[END_REF] , one magmatic pathway of Tolbachik appears to be connected with the marginal part of the Klyuchevskoy deep reservoir, and another seems to originate from an independent mantle source located to the south of Tolbachik. s -1 used for automatized detection of eruptions (see Methods section). The start and end of each eruptive period (labelled on the top of each plot) are marked by the vertical dashed lines. The error for every single datapoint is ± 50%. The time series is generated using MATLAB software (www.mathworks.com).

Remote sensing of eruption effusion rates

We calculated the Time-Averaged lava Discharge Rate (TADR) and the erupted lava volumes at the three volcanoes by using MODIS infrared data acquired between March 2000 and December 2019 (bold lines in insets of Fig. 1b,c,d) determined with the MIROVA system [START_REF] Coppola | Thermal Remote Sensing for Global Volcano Monitoring: Experiences From the MIROVA System[END_REF][START_REF] Coppola | Enhanced volcanic hot-spot detection using MODIS IR data: results from the MIROVA system[END_REF] . Details of the methodology and associated limits are described in the Methods section accompanying this paper. A TADR threshold of 0.25 m 3 s -1 is used to automatically recognize the main eruptive periods at each volcano (Fig. 2), and to quantify the eruption parameters summarized in Table 1.

Klyuchevskoy

Nine eruptions occurred at Klyuchevskoy between 2003 and 2020 (Fig. 2a). Of these, eight were automatically recognized (see Methods), and one was manually selected, based on observations of Ozerov et al. [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] (eruption #6; Table 1). Most of the eruptions (#2,3,4,5,7,8,9; Table 1a) produced lava flows along the flanks of the volcano [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] and created lava volumes ranging from ~10 to 150 × 10 6 m 3 each, with a Mean Output Rate (MOR: total volume of eruption/duration) ranging between 2.5 and 10 m 3 s -1 and a maximum TADR often higher than 30 m 3 s -1 (Table 1a). Only two eruptions (#1,6; Table 1a) were limited to moderate explosive activity inside the summit crater [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] characterized by much lower volumetric output (< 10 Mm 3 ) and discharge rates (maximum TADR < 2.5 m 3 s -1 ; Table 1a). For some eruptions (#3,4,5; Fig. 2a), the onset of lava effusion was preceded by a precursory phase of several weeks, identified by increased fumarolic activity and degassing [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] .

In other cases, the beginning of the eruption was rather rapid, without any apparent thermal precursory phase (#2,7,8,9; Fig. 2a). The eruptive trends of Klyuchevskoy are often characterized by a TADR that increases with time to reach values of 10-100 m 3 s -1 immediately before the effusion suddenly ceases (Fig. 2a). The volumetric output of the 20 years (Fig. 2) defines the most recent period of intense activity of Klyuchevskoy characterized by a steady-state output rate (Qss) of 1.36 m 3 s -1 (1.21 m 3 s -1 DRE; Fig 3a), which is almost twice the average output since 1930 (Fig. 1b1). Notably, the cumulative curve in Fig. 3a shows a clear sawtooth pattern typical for steady-state volcanism [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] whereby each step is either produced by (i) an unbuffered arrival and eruption of discrete magma batches, or (ii) a partial or complete discharge of a shallow reservoir that is fed by a constant magma supply. A similar pattern could be also explained if the arrival of discrete magma batches is controlled by a steady-state destabilization of magma reservoirs, produced by passive degassing during quiescence, which can trigger magma ascent from depth [START_REF] Girona | On depressurization of volcanic magma reservoirs by passive degassing[END_REF][START_REF] Girona | Degassing during quiescence as a trigger of magma ascent and volcanic eruptions[END_REF][START_REF] Mittal | Volatile degassing from magma chambers as a control on volcanic eruptions[END_REF] . Whatever the model, the two lines, which envelope the sawtooth curve (parallel to the linear trend ± 2σ; Fig. 3a), define the maximum size (maximum eruptible lava volume) and maximum response time of the individual eruption, respectively [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] . For Klyuchevskoy, these values are approximately 143 × 10 6 m 3 and 1214 days. The analysis of inter-eruption time distribution (Methods) suggests a relatively strong periodicity (Fig. 4a). The degree of periodicity can be quantified by the coefficient of variation (CV), which is zero for perfect periodicity, one for randomness, and larger than one for clustering. In the case of Klyuchevskoy's eruptions, we found a CV equal to 0.38 and the inter-eruption time positively correlated to the size of the last event (correlation coefficient of 0.65), as expected for time-predictable systems (Fig. 4a). A load and discharge model is thus envisaged for Klyuchevskoy (Fig. 3), whereby an eruption starts when the upper, critical volume threshold is accumulated in the shallow reservoir [START_REF] Cruz-Reyna | Poisson-distributed patterns of explosive eruptive activity[END_REF] . The resumption of eruptive activity on November 2019 [START_REF]Global Volcanism Program. Report on Klyuchevskoy[END_REF] further supports a time-predictable behavior (Methods), which is in agreement with the achievement of a critical volume as shown in Fig. 3a.

Bezymianny

Thermal data acquired over Bezymianny (Fig. 2b) are indicative of an open-vent volcano, persistently emitting hot volcanic products. The retrieved long-term eruptive pattern can be subdivided into two distinct regimes:

(1) a continuous low-level regime, associated with passive degassing and possibly related to "endogenous growth," and ( 2) an intermittent high-level regime, associated with short-term (days to weeks) extrusiveexplosive-effusive cycles. A TADR threshold of 0.25 m 3 s -1 separates the two regimes and automatically recognizes 28 out of the 29 major eruptive cycles between 2000 and 2019 [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF][START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] (Table 1b). The only undetected event occurred on 5 November 2007, when strong cloud coverage over the volcano prevented the detection of this short-lived event (Table 1). Each eruption cycle is characterized by erupted volumes ranging from ~0.15 to ~15 × 10 6 m 3 and peak TADRs between 0.35 to 52 m 3 s -1 (Table 1b). The average duration of each eruptive cycle is 26.7 (± 20.7; 1σ) days, much shorter than the average inter-eruption time of 222.7 days. The cumulative volume curve of Bezymianny is essentially controlled by the sudden steps associated with the eruptive cycles detected by MODIS (Fig. 3c). Between 2000 and mid-2002, eruptions reached higher TADR peaks, causing a steeper cumulative volume curve than in the rest of the time series (Fig. 3c). Although this may reflect a higher magma output rate in this period, it is also possible that the dataset is biased by the fact that only one satellite was operating in that period (see Methods). Similar to Klyuchevskoy, the cumulative curve of the volumes erupted between May 2002 and September 2012 suggests steady-state volcanism for Bezymianny [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] , which is characterized by an average output rate of ~0.21 m 3 s -1 , a maximum eruptible volume of ~7.6 × 10 6 m 3 , and a maximum repose time of 410 days (Fig. 3c and3d). A notable lack of eruptive cycles occurred between September 2012 and December 2016 (Fig. 2b). This anomalously long rest period (low thermal regime) is also visible in Fig. 3b and3c, where the cumulative volume curve diverges horizontally from the steady-state model. According to Wadge [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] , this pattern occurs at steady-state volcanoes when magma is not being supplied into the shallow reservoir, here either because magma is not being generated or a neighboring volcano is capturing it. Bezymianny's activity resumed at the end of 2016 and continued intermittently with an output rate similar to the 2002-2012 period. The analysis of inter-eruption times (Methods) suggests a quasi-periodic behavior (CV=0.5) for Bezymianny's activity until 2012 (Fig. 4b), which is completely lost when including the whole dataset (Fig. 4c). No correlation is found between the intereruption times and volumes released during the last or the next eruption (Fig. 4). Although the lack of correlation can be due to the significant uncertainties affecting the Bezymianny time series (Methods), we may not exclude the role of a time-varying upper threshold (strength) of the shallow magmatic system [START_REF] Cruz-Reyna | Poisson-distributed patterns of explosive eruptive activity[END_REF] . Bezymianny. The variation from the steady-state is calculated as the residual between the linear growth model and the observed cumulative volume curves. A two-fold standard deviation (± 2σ) of residual is used to define the upper and lower volume limits for the load-discharge model. In this plot, a time-predictable system would have a near-constant upper threshold of the volume load when eruptions occur. A volumepredictable system would have a near-constant lower threshold. An entirely predictable system would have both. The gray bar indicates the timing of the 2012 Tolbachik Fissure Eruption (TFE). Time series generated using MATLAB software (www.mathworks.com).

Tolbachik

The eruption of Tolbachik (November 2012 -August 2013) commenced suddenly on 27 November 2012, producing an initial TADR peak of about 300 m 3 s -1 . This initial activity resulted in the emplacement of a lava flow that reached a length of about 15 km in a few days [START_REF] Ramsey | The influence of emissivity on the thermosrheological modeling of the channelized lava flows at Tolbachik volcano[END_REF] . Effusion rates decrease roughly exponentially during the following ten months of continuous effusive activity. The eruption stopped between 23 and 27 August 2013 when the TADR suddenly lowered from 7-9 to less than 0.25 m 3 s -1 . The eruption emplaced a volume of ~0.5 × 10 6 m 3 in 302 days [START_REF] Dvigalo | The first quantitative estimates of parameters for the Tolbachik Fissure Eruption of 2012-2013 from aerophotogrammetric observations[END_REF] , which is up to 4000 times the volume commonly erupted at Bezymianny volcano (>0.15 ×10 6 m 3 , see above). (2000-2012). The coefficient of variation (CV), a measure for the eruption's periodicity, is provided in the title of each plot. Symbols show the eruption volume (right axis) in all three panels as a function of the inter-eruption time, where points and crosses refer to the volume of the last and next event, respectively. The correlation coefficient (r) between volume and interevent-time with its p-value is provided in the legends.

Interactions between Klyuchevskoy, Bezymianny, and Tolbachik

We statistically explore if the three volcanoes interacted on more than one occasion and in different ways. Specifically, we found various degrees of interactions that become best observable by analyzing the data at the time scales from weeks to decades.

Conjoint activity and pattern's change before and after the 2012 Tolbachik eruption

A first indication of how volcanoes are interconnected with each other is revealed by the detailed analysis of the mutual activity of Klyuchevskoy and Bezymianny (and pattern's change) before and after the Tolbachik eruption (Fig. 5).

Before the latter (Fig. 5a1), we observe a simultaneous activation of Bezymianny and Klyuchevskoy several times (i.e., eruptions KLY#1,2,6) while no simultaneous activation is found afterward (Fig. 5a2). In particular, the onset of Klyuchevskoy's eruptions #1,2,6 coincided with the maximum activity of Bezymianny (BEZ#6,8,22) in the same period (a time window of ±10 days is considered to avoid the effect of clouds).

Although less significant, Bezymianny's activity also showed some synchronous activation (increase of TADR relative to the previous trend) with the onset of the other Klyuchevskoy eruptions (KLI#3,4) before the Tolbachick eruption.

More specifically, we found that in the days-to-weeks following each of Klyuchevskoy eruptions, the average TADR of Bezymianny increased, on average, by a factor of four (Methods). This increase suggests that before 2012, the eruptions of Klyuchevskoy were able to "galvanize" also the activity of Bezymianny. In contrast, after the eruption of Tolbachik, the two volcanoes no longer have erupted simultaneously.

Additionally, before the Tolbachik eruption, most of the eruptions of Klyuchevskoy were characterized by a precursory phase marked by a gradual increase in thermal activity and estimated TADR (Fig. 5b1). This preeruptive pattern is typical of open-vent volcanoes, in which the rise of the magma column causes the appearance and growth of fumaroles or weak explosive activity [START_REF] Laiolo | Space-and Ground-Based Geophysical Data Tracking of Magma Migration in Shallow Feeding System of Mount Etna Volcano[END_REF] . However, the precursory pattern disappeared after the eruption of Tolbachik (Fig. 5b2), and all the three subsequent eruptions of Klyuchevskoy showed a sudden beginning of activity more typical of closed-vent systems [START_REF] Laiolo | Space-and Ground-Based Geophysical Data Tracking of Magma Migration in Shallow Feeding System of Mount Etna Volcano[END_REF] . Similarly, at Klyuchevskoy, the precursory activity is evident for the eruptions that occurred before the Tolbachik eruptions (b1), but not afterward (b2). Time series generated using MATLAB software (www.mathworks.com).

Reactivation of Tolbachik in 2011-2012

The comparison between the surface activity, retrieved from satellite, and the long-period (LP) earthquakes occurred within the KVG during the reactivation of Tolbachik [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF][START_REF] Frank | Progressive reactivation of the volcanic plumbing system beneath Tolbachik volcano (Kamchatka, Russia) revealed by long-period seismicity[END_REF] , provides further indications of mutual communication between these volcanoes. This last was preceded in 2011 by an increase of the Deep Long Period (DLP) seismic activity reaching its maximum level in May 2012 (Fig. 6c). It reflected the gradual pressurization of the whole KVG plumbing system [START_REF] Frank | Progressive reactivation of the volcanic plumbing system beneath Tolbachik volcano (Kamchatka, Russia) revealed by long-period seismicity[END_REF] , possibly in response to a pulse of volatile-rich basaltic magmas rising from the mantle [START_REF] Melnik | Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas[END_REF] . At Bezymianny, this gradual pressurization may have triggered three consecutive shallow LP swarms, each preceding an eruption, the last one being in September 2012 (Fig. 6a,b).

Similarly, LP seismicity also migrated shallow below the Klyuchevskoy volcano in September 2012 (promptly triggering the onset of eruption #6), and later, on October-November 2012, LPs occurred below Tolbachik, just before the onset of its voluminous flank eruption (Fig. 6a,b). Interestingly, during eruption #6 of Klyuchevskoy, the TADR trend almost mirrors the shallow LP seismicity, both reaching a maximum in November 2012 and then declining in correspondence with an acceleration of the Tolbachik seismic swarm (Fig. 6a,b). Eruption #6 was somehow unusual for Klyuchevskoy since it produced only weak Strombolian activity [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] , with a TADR always below 0.25 m 3 s -1 and a volume of less than 1 × 10 6 m 3 (Table 1). Moreover, unlike the other Klyuchevskoy eruptions (cf. Fig. 2a), it never culminated in effusive activity, which is atypical for this volcano. Together with a waning trend of surface and seismic activity since mid-November 2012, these peculiar features suggest a sort of partial depletion of the shallow magma supply of Klyuchevskoy, precisely in correspondence with the acceleration of seismic swarms below Tolbachik (Fig. 6a,b). It is worth noting that the eruptions of Klyuchevskoy stopped in 1975-1976 during the GFTE and were renewed in 1977-1978 after the GFTE [START_REF] Fedotov | Seismic Processes and Migration of Magma during the Great Tolbachik Fissure Eruption of 1975-1976 and Tolbachik Fissure Eruption of 2012-2013, Kamchatka Peninsula[END_REF] .

Reactivation of Klyuchevskoy and cessation of Tolbachik eruption in August 2013

Even more intriguing is the resumption of the activity at Klyuchevskoy (eruption #7) and the almost concurrent cessation of activity at Tolbachik on 22 August 2013 (Fig. 6a). The beginning of eruption #7 occurred suddenly on 14 August 2013, with the onset of Strombolian explosions, which evolved in few days into summit effusive activity [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF] fed with a TADR of about 10 m 3 s -1 (Fig. 6a). Lava discharge rates increased rapidly in the following months to reach a maximum value of ~100 m 3 s -1 on 18 October 2013, just before the abrupt cessation of surface activity on 25 October 2013.

The onset of eruption #7, which also occurred abruptly on 14 August 2013, preceded the end of the Tolbachik eruption by eight days (Fig. 6a). Our data suggest that the Tolbachik eruption ended when the TADR-values were still moderately high (7-9 m 3 s -1 ), shutting down the monthly-long, (almost) exponential decay.

Reactivation of Bezymianny in 2016

Strong evidence for volcano-volcano interactions is the lack of the typical extrusive-explosive-effusive cycles of Bezymianny for four years after the eruption of Tolbachik [START_REF] Girina | Satellite high-resolution data used to clarify the position of fault zones within the Klyuchevskaya volcanic group of Kamchatka. Sovremennye problemy distantsionnogo zondirovaniya[END_REF] (Fig. 7a). This rest period was unusually long for Bezymianny (1550 days) and started already on 11 September 2012 (~3 months before Tolbachik). As discussed above, the September 2012 eruption of Bezymianny represents the superficial response of its plumbing system to the main deep magma pulse revealed by DVLP, which heralded a few months later the eruption of Tolbachik. The following lack of activity at Bezymianny persisted for four years during which continuous thermal anomalies were likely related to passive degassing (Fig. 7). In early 2016, a viscous, crystallized, cold plug started to be extruded from the summit crater [START_REF] Mania | Deformations and Morphology Changes Associated with the 2016-2017 Eruption Sequence at Bezymianny Volcano[END_REF] . This slow, cold extrusion was undetected by MODIS but, according to Mania et al. [START_REF] Mania | Deformations and Morphology Changes Associated with the 2016-2017 Eruption Sequence at Bezymianny Volcano[END_REF] , accelerated in September-November 2016 (right at the end of Klyuchevskoy eruption #9) until the effusion of a viscous lava flow on 9 December 2016 (eruption #23; Fig. 7). The extrusion of solid plugs at the onset of eruptive cycles is a typical feature of Bezymianny [START_REF] Ozerov | Eruptions in the Northern Group of Volcanoes, in Kamchatka, during the Early 21st Century[END_REF][START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] .

However, that of 2016 represented an abnormally long precursory phase for this volcano. It was followed by eruptions #24 and #25, both characterized by a gentle effusion of two lava flows with increasingly stronger explosivity [START_REF] Mania | Deformations and Morphology Changes Associated with the 2016-2017 Eruption Sequence at Bezymianny Volcano[END_REF] . This peculiar dynamic after four years of rest seems to be consistent with an interruption (or decrease) of the magma supply after the TFE that favored the formation of a cold crystallized plug in the shallow conduit of Bezymianny.

Influence of the 1955-56 eruption of Bezymianny and the 1975 Great Tolbachik Fissure Eruption in the long term eruptive pattern of Klyuchevskoy and Bezymianny

Bezymianny eruption in 1955-56 (BEZ55) was the largest in the recent history on this volcano [START_REF] Gorshkov | Gigantic eruption of the volcano Bezymianny[END_REF] and its occurrence may have perturbed the whole KVG in some way. Some evidence appears from the analysis of the long term volumetric output of Klyuchevskoy (Fig. 1b) which shows an evident reduction in the eruption' frequency after 1956, passing from 0.28 events/year, between 1930 and 1953, to 0.17 events/year between 1956 and 1973 (Fig. 1b). Even more indicative is the fact that this decrease was also associated with an evident geochemical change in the products erupted by Klyuchevskoy after 1960 [START_REF] Bergal-Kuvikas | A petrological and geochemical study on time-series samples from Klyuchevskoy volcano, Kamchatka arc[END_REF] (a few years after the unrest of Bezymianny), which has been ascribed to the injection of new type of primary magma that was not produced beneath the volcano previously.

In 1977, immediately after the GTFE, the eruptive regime of Bezymianny changed considerably to give rise, for the first time, to the effusion of lava flows and the establishment of extrusive-explosive-effusive cycles [START_REF] Bogoyavlenskaya | Eruption of Bezymianny volcano in 1977[END_REF][START_REF] Shevchenko | The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse[END_REF] . Simultaneously, the volumetric output rate decelerated (Fig. 1c). The erupted magma became more and more primitive, indicating the arrival of deeper mafic magma components at the surface [START_REF] Turner | The magma plumbing system of Bezymianny Volcano: Insights from a 54year time series of trace element whole-rock geochemistry and amphibole compositions[END_REF] . For Klyuchevskoy, the available data and observations point to a significant change in its eruptive regime starting 1-3 years after the GTFE, when summit eruptions began to dominate over flank eruptions [START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] (Fig. 1b). In contrast to Bezymianny, a significant increase of the output rate accompanied this change likely associated with an increased magma supply at shallower levels. Geodetic measurements 37 carried out between 1979 and 2005 also suggest that the magma feeding system of Klyuchevskoy was accumulating considerable amounts of magma before the eruptions of this period, which is coherent with a gradual rise of the effective pressure source, from the probable region of deep magma storage (25 km) to shallow levels (5 km). Also the seismicity of the entire KVG showed a dramatic change after the GTFE eruption [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF][START_REF] Fedotov | Seismic Processes and Migration of Magma during the Great Tolbachik Fissure Eruption of 1975-1976 and Tolbachik Fissure Eruption of 2012-2013, Kamchatka Peninsula[END_REF] with several remarkable earthquake swarms that occurred during 1977-1978. All these pieces of evidence make it plausible that both the BEZ25 and GTFE have perturbed the activity of the neighboring volcanoes, compatibly with a process of general rejuvenation of the whole KVG magma system. Whether a similar rejuvenation process occurred during the 2012 Tolbachik eruption is still unclear. However, the occurrence of the major swarms of deep very-long period events (DVLP) in 2011 and 2012 (Fig. 6) point toward a process of pre-eruptive reload of the shallow magmatic reservoirs from depth [START_REF] Chebrov | Development of digital seismic network in Kamchatka, in Geofizicheskii monitoring Kamchatki[END_REF] .

Discussion

Our new satellite data suggest that the three volcanoes of the KVG are related to each other on various timescales. The mode and directivity of the relation vary, showing correlated and anti-correlated activity changes. This observation probably reflects a complex response to changes occurring in a seismically inferred common magmatic source and/or at the associated hydrothermal system.

Conjecturing the presence of crustal magma chambers at the volcano systems, we may develop a simple conceptual model to explain some of the modulations and concurrent activity changes observed in our data. At Klyuchevskoy, the magma supply within the crustal plumbing system follows a general steadystate load and discharge model. The frequent but intermittent arrival of magma batches is buffered by the elastic deformation of the subvolcanic reservoir [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] . Eruptions occur when the stored amount of magma exceeds a specific threshold (time-predictable behaviour [START_REF] Cruz-Reyna | Poisson-distributed patterns of explosive eruptive activity[END_REF] ) with the maximum eruptible volume (~150 ×10 6 m 3 for Klyuchevskoy; Fig. 3a2) strictly connected to the capacity of the reservoir to buffer the arrival of magma [START_REF] Turner | The magma plumbing system of Bezymianny Volcano: Insights from a 54year time series of trace element whole-rock geochemistry and amphibole compositions[END_REF] . It is interesting to note that during the steady-state regime, the magma ascent feeding the activity at these volcanoes could be driven by processes occurring at depth 61 (down-top mechanism), but also by the passive degassing during quiescence [START_REF] Girona | On depressurization of volcanic magma reservoirs by passive degassing[END_REF][START_REF] Girona | Degassing during quiescence as a trigger of magma ascent and volcanic eruptions[END_REF][START_REF] Mittal | Volatile degassing from magma chambers as a control on volcanic eruptions[END_REF] , which induces the opening of pathways connecting deep and shallow magma reservoir (top-down-mechanism).

The eruptive behavior of Bezymianny is also compatible with a steady-state magma supply. However, in this case, the smaller capacity of the reservoir(s) and the lower magma supply rate (compared to Klyuchevskoy) give rise to much more frequent but less voluminous eruptions (maximum eruptible volume ~7 ×10 6 m 3 ; Fig. 3b2). In this steady-state framework, the volcanoes' conjoint activation indicates that both systems responded to a common perturbation, possibly sourced at lower crustal levels.

On the other hand, significant large swarms of DVLP (Fig. 6), ascribed to deep magma pulses, can reactivate the Tolbachik magmatic path [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF][START_REF] Frank | Progressive reactivation of the volcanic plumbing system beneath Tolbachik volcano (Kamchatka, Russia) revealed by long-period seismicity[END_REF][START_REF] Melnik | Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas[END_REF][START_REF] Fedotov | Seismic Processes and Migration of Magma during the Great Tolbachik Fissure Eruption of 1975-1976 and Tolbachik Fissure Eruption of 2012-2013, Kamchatka Peninsula[END_REF] , which in turn modify the properties of the of nearby magmatic systems and perturb their steady-state regime.

Deviation from the steady-state cumulative volume curve indicates a change in the magma supply rate [START_REF] Wadge | Steady state volcanism: Evidence from eruption histories of polygenetic volcanoes[END_REF] , as occurred after the BEZ55 and the GTFE at both Klyuchevskoy and Bezymianny.

The GTFE eruption directly affected Bezymianny's activity, causing a reduction of the magma output rate since 1977 (Fig. 1b2) and producing a radical change in Bezymianny's eruptive regime and a rejuvenation of its eruptive products [START_REF] Turner | The magma plumbing system of Bezymianny Volcano: Insights from a 54year time series of trace element whole-rock geochemistry and amphibole compositions[END_REF] . Similarly, but in the opposite direction, the GTFE led to an increase of the long-term magma output rate of Klyuchevskoy and promoted a change in its eruptive pattern, switching from lateral to summit eruptions (Fig. 1b1).

To a lesser extent, the reactivation of the Tolbachik in 2012 inhibited the steady-state magma supply of Bezymianny for several years. It caused the interruption of its surface activity until the extrusion of a crystallized plug in 2016 (Fig. 7). During this period, multiple interactions between Tolbachik and Klyuchevskoy were also observed, supporting the existence of a very efficient hydraulic connection between the plumbing systems of the three volcanoes. We note that the details on the presence of a common primary magma feeding all volcanoes in KVG as well as the location and geometry of crustal magma chambers are still debated [START_REF] Kayzar | Deciphering petrogenic processes using Pb isotope ratios from time-series samples at Bezymianny and Klyuchevskoy volcanoes, Central Kamchatka Depression[END_REF][START_REF] Ozerov | Petrological-geochemical model for genetic relationships between basaltic and andesitic magmatism of Klyuchevskoi and Bezymyannyi volcanoes, Kamchatka[END_REF] , which is why our conceptual model remains speculative. Shapiro et al. [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer[END_REF] proposed a model based on fluid-pressure propagation through porous rocks to explain the migration of LP events and infer the existence of such hydraulic connections below the KVG volcanoes. Our data supports and reinforces this hypothesis, although other processes can also influence the eruption cycles and dynamics [START_REF] Sulpizio | Influence of stress field changes on eruption initiation and dynamics: a review[END_REF] . Independent of any model assumption, our data show that the magmatic systems below the KVG are interconnected, and eruptions of individual volcanoes can be the direct consequence of their neighbors' activity.

To what extend magmatic systems are connected and if one eruption can trigger another volcano are essential questions for assessing volcanic hazard. In the case of interacting volcanoes, such as in the case of KVG, a volcano's behavior can be the direct consequence of its neighbor's activity. In these cases, traditional hazard assessments of isolated volcanoes have to be replaced by a comprehensive assessment involving the whole volcanic group. In addition to its eruptive history, the volcano's hazard assessment has to account for its neighboring volcanoes' eruptive history, which may influence its current state.

Methods

Satellite Thermal data

Satellite thermal data were processed using the MIROVA system 60 (www.mirovaweb.it), which is based on the analysis of the images acquired by MODIS. The two MODIS sensors, launched in March 2000 and May 2002, provide approximately six infrared images per day over Kamchatka (three night-times and three daytimes) with a nominal ground resolution of 1 km. MODIS images are processed at each volcano to quantify the Volcanic Radiative Power (VRP in Watts), a combined measurement of the area and integrated temperature of the hot (>200°C) volcanic features with a standard error of ± 30% over every measurement [START_REF] Coppola | Enhanced volcanic hot-spot detection using MODIS IR data: results from the MIROVA system[END_REF] .

We used only the night-time MODIS dataset, consisting of approximately 19,500 images acquired over the Klyuchevskoy Volcanic Group (KVG). Thermal anomalies detected by MIROVA were geolocated (errors in geolocation are less than 0.5 km for nadir acquisition [START_REF] Coppola | Enhanced volcanic hot-spot detection using MODIS IR data: results from the MIROVA system[END_REF] ) to discriminate the hotspots sourced by the three distinct volcanoes. All the images were visually analyzed to discard the data contaminated by clouds, ash plumes, or poor viewing conditions (i.e., high satellite zenith), which preclude a correct estimation of VRP [START_REF] Dzurisin | Magma supply to Kīlauea Volcano, Hawai'i, from inception to now: Historical perspective, current state of knowledge, and future challenges[END_REF][START_REF] Shcherbakov | Plagioclase zoning as an indicator of magma processes at Bezymianny Volcano, Kamchatka[END_REF] .

Finally, the supervised dataset consists of 2,139 images for Klyuchevskoy, 2013 images for Bezymianny, and 219 images for Tolbachik, which have been used to reconstruct the time-series of VRP (Fig. S1 -Supplementary Material). For each volcano, the cumulative Volcanic Radiative Energy (VRE) in Joules is calculated as the trapezoidal integration of the supervised VRP time series (Fig. S1 -Supplementary Material).

Erupted volume and Time-Averaged lava Discharge Rate

We used a simplified approach, which has been expressly developed to derive Time Averaged lava discharge Rate (TADR) directly from MODIS-derived VRP [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF] . This approach assumes that during an eruption, the energy radiated by a lava body (i.e., VRE) is linearly correlated to the bulk erupted volume (Vol), 1), where c rad (in J m -3 ) is the best-fit coefficient that describes the ability to radiate thermal energy by unit volume of the observed lava body. Thus the c rad value can be determined retrospectively by measuring the energy radiated during an eruption (or during an eruptive period) and the bulk volume of the lava flow(s) or domes emplaced during the same time interval (measured independently).

𝑐 "#$ = &'( &)* (eq.
Once calibrated, the c rad coefficient is used to retrieve the TADR for each single VRP measurements according to 𝑇𝐴𝐷𝑅 = &'/ 0 123 (eq. 2).

Note that this approach does not take into account the volume of magma erupted explosively (i.e., ash plumes, pyroclastic density currents). It accounts only for magma erupted during effusive/extrusive periods, that is, when sufficient thermal radiation is detectable from the satellite.

To estimate the c rad -value of Klyuchevskoy, we considered the period between 2002 and 2009, during which about 231 x 10 6 m 3 of lava erupted [START_REF] Fedotov | The magmatic system of the Klyuchevskaya group of volcanoes inferred from data on its eruptions, earthquakes, deformation, and deep structure[END_REF] . Assuming an average fraction of tephra equal to 15% in volume [START_REF] Fedotov | On the eruptions, deformation, and seismicity of Klyuchevskoy Volcano, Kamchatka in 1986-2005 and the mechanisms of its activity[END_REF] , the cumulative volume of lava flows erupted between 2002 and 2009 become ~196 x 10 6 m 3 . This activity produced a VRE of 1.6 x 10 16 J (Fig. S1c1 -Supplementary Material), which results into an average c rad -value of 8.16 × 10 7 J m -3 .

For Bezymianny volcano, we calibrated the c rad , by considering the dome volume's growth between 31 July 2006 and 9 September 2017 [START_REF] Shevchenko | The rebirth and evolution of Bezymianny volcano, Kamchatka after the 1956 sector collapse[END_REF] . Given a total volume of ~ 69 x 10 6 m 3 and a VRE of 1.17 x 10 15 J (Fig. S1c2 -Supplementary Material), we estimated c rad = 1.88 x 10 7 J m -3 . Note that the TADR and inferred volumes do not include the contribution of the explosive activity, which in the case of Bezymianny may be relevant.

According to Girina et al. [START_REF] Girina | Chronology of Bezymianny Volcano activity, 1956-2010[END_REF] , each extrusive-explosive-effusive cycle produces volumes up to ~10 7 m 3 , in the form of pyroclastic flows. Although the amount of juvenile material inside these deposits is unknown, the large amount of material erupted explosively, together with the short duration of each cycle, introduces a significant level of noise into our time series and an uncertainty possibly higher than 100% in the volumes reported in Table 1b.

Equations 1 and 2 have been successfully applied to estimate the TADRs of the 2012-2013 Tolbachik eruption [START_REF] Ramsey | The influence of emissivity on the thermosrheological modeling of the channelized lava flows at Tolbachik volcano[END_REF] , where a c rad equal to 1.08 × 10 8 J m -3 has been calculated based on a final lava flow volume [START_REF] Dai | Measuring lava flows with ArticDEM: Application to the 2012-2013 eruption of Tolbachik, Kamchatka[END_REF] of 573 × 10 6 m 3 and a corresponding VRE equal to 6.07 × 10 16 (Fig. S1c2 -Supplementary Material).

As described by Coppola et al. [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF] , this approach provides single TADR measurements with an associated error of ±50%. Errorbars are not shown for graphical convenience.

Statistical testing of correlated activity

The frequency plot of inter-eruption time (dt_es in Tab 1) for Bezymianny and Klyuchevskoy is shown in the left axis of Fig. 4. The peaked distribution for Bezymianny data (gray bars in Fig. 4a,b) can be reasonably fitted by a Brownian-Passage Time (BPT) distribution (blue line). This models assumes a fixed eruption threshold and volume release, plus a constant loading rate with noise. The coefficient of variation (CV), also called the aperiodicity parameter, measures a signal's periodicity, where CV=0 refers to a perfect periodicity, CV=1 to a random Poisson occurrence, and CV>1 to clustering. When considering the whole dataset of Bezymianny (Fig. 4a), the CV value is 1.20, indicating a random occurrence of eruptions. However, the CV value decreases to 0.50 when post-Tolbachik eruption data are excluded (Fig. 4b), thus indicating a quasi-periodic behavior until the Tolbachik eruption. For Klyuchevskoy (Fig. 4c), the CV value is even lower (CV=0.38), indicating a rather clock-wise recurrence of eruptions.

The relation between inter-eruption time and volume release of the last or next event is shown on the right axis of Fig. 4. The data shows no correlations for Bezymianny (p-values > 0.25 in Fig. 4a,b), while p<0.05 would indicate a statistically significant correlation. In contrast, Klyuchevskoy (Fig. 4c) shows a positive correlation between the inter-eruption time and the volume of the last eruption (r=0.65, p=0.081), which become statistically significant (r=0.73; p=0.025) when including the timing of the last eruption (not included in our study) started in November 2019 [START_REF]Global Volcanism Program. Report on Klyuchevskoy[END_REF] . The weaker correlation with the next events' volume found for the Kluchevskoi volcano (Fig. 4c) suggests that its eruption periodicity is consistent with a time-predictable rather than a volume-predictable model.

Conjoint activity of Klyuchevskoy and Bezymianny volcanoes

We tested the hypothesis that before the Tolbachik eruption, the activation of Klyuchevskoy (eruptions #1 to 6) affected the Bezymianny activity (Fig. S4). To perform this test, we first calculated the average TADRvalue (of Bezymianny) within T days (from 10 to 30 days) relative to each Klyuchevskoy eruption and averaged those six values. Then we calculated the ratio between the averaged TADR-value in the T days after the eruption and the corresponding value in the T days before the eruption to measure the average activation (blue points in figure 8). Finally, we compared the observed ratio (as a function of T) with the corresponding result obtained after randomizing the activation times of the six Klyuchevskoy eruptions within the period between 2002 and the Tolbachik eruption. The fraction of randomized data with a ratio similar or larger than the observed one (green line in Fig. 8) shows that the observed activation value can only be reached in less than 5-10% of the randomized data. Although the results are close to the significance threshold, these data suggest that the result is significant for the shortest time intervals (i.e., T=10 days; p<0.05) with a 4-5 times increase of the averaged TADR of Bezymianny after the onset of a Klyuchevskoy eruption.

Note that we have not analyzed T<10 days because of missing Bezymianny measurements during short periods before/after the Klyuchevskoy eruptions.

Tables

Figure 1 .

 1 Figure 1. (a) Distribution of volcanic centers within the Klyuchevskoy Volcanic Group (KVG). Right panels show the Cumulative Dense Rock Equivalent (DRE) volumes of lavas erupted between 1930 and 2020 for (b) Klyuchevskoy 21 , (c) Bezymianny 44 , and (d) Tolbachik 57 . All volumes are recalculated for DRE using: (i) a magma density of 2800 kg m -3 and an average density of lava flows of 2500 kg m -3 for Klyuchevskoy and Tolbachik 57 ; (ii) a magma density of 2500 kg m -3 and an average density of lava dome of 2000 kg m -3 for Bezymianny 41 . A standard error of ±20% (colored fields) takes into account uncertainties in the estimates of bulk volumes and densities of erupted products. Red bold lines within insets, correspond to satellite-derived volumes (this work; see Methods). The grey bars outline the timing of the Bezymianny eruption in 1955-56 (BEZ55), the Great Tolbachik Fissure Eruption in 1975 (GTFE) and the Tolbachik Fissure Eruption in 2012 (TFE). Shaded relief map derived from ArcticDEM digital elevation model (https://www.pgc.umn.edu/data/arcticdem/) and elaborated using QGIS version 3.16.3 (http://qgis.osgeo.org); Time series generated using MATLAB software (www.mathworks.com).

Figure 2 -

 2 Figure 2 -Time-series of the Time-Averaged lava Discharge Rate (TADR; logarithmic scale) for (a) Klyuchevskoy, (b) Bezymianny, and (c) Tolbachik. Single satellite measurements (dots) are interpolated and smoothed (line) to provide continuous data. The horizontal dashed lines represent the threshold of 0.25 m 3 s -1 used for automatized detection of eruptions (see Methods section). The start and end of each eruptive period (labelled on the top of each plot) are marked by the vertical dashed lines. The error for every single datapoint is ± 50%. The time series is generated using MATLAB software (www.mathworks.com).

Figure 3 -

 3 Figure 3 -Upper panels. Cumulative volume curves for (a1) Klyuchevskoy and (b1) Bezymianny, derived from satellite measurements (see Methods). The linear fits of each cumulative curve represent the linear growth model (blue thick dashed lines) and provide the steady-state output rate (Qss). The two parallel lines (linear model ± 2σ) define the maximum eruptible volume (Vol max ) and the maximum response time (t max ) for the inferred steady-state condition. Lower panels. Load and discharge model for (a2) Klyuchevskoy and (b2)Bezymianny. The variation from the steady-state is calculated as the residual between the linear growth model and the observed cumulative volume curves. A two-fold standard deviation (± 2σ) of residual is used to define the upper and lower volume limits for the load-discharge model. In this plot, a time-predictable system would have a near-constant upper threshold of the volume load when eruptions occur. A volumepredictable system would have a near-constant lower threshold. An entirely predictable system would have both. The gray bar indicates the timing of the 2012 Tolbachik Fissure Eruption (TFE). Time series generated using MATLAB software (www.mathworks.com).

Figure 4 -

 4 Figure 4 -Inter-eruption times distribution (bars) for (a) Klyuchevskoy, (b) Bezymianny -whole period (2000-2020), (c) Bezymianny -pre-Tolbachik eruption(2000-2012). The coefficient of variation (CV), a measure for the eruption's periodicity, is provided in the title of each plot. Symbols show the eruption volume (right axis) in all three panels as a function of the inter-eruption time, where points and crosses refer to the volume of the last and next event, respectively. The correlation coefficient (r) between volume and interevent-time with its p-value is provided in the legends.

  Below we provide evidence of interactions related to (i) conjoint activity of Klyuchevskoy and Bezymianny throughout 2003-2020, (ii) the reactivation of Tolbachik in 2011-2012, (iii) the reactivation of Klyuchevskoy and the cessation of Tolbachik in August 2013, (iv) the reactivation of Bezymianny in 2016-2017, and (v) changes in the long-term magma output rate after the Bezymianny eruption, in 1955-56, and after the Great Tolbachik Fissure Eruption (GFTE), in 1977.

Figure 5 -

 5 Figure 5 -Stacked TADR time-series relative to the onset of the Klyuchevskoy eruptions, before (left column) and after (right column) the Tolbachik eruption. Upper panels: TADRs of the Bezymianny volcano. Lower panels: TADRs of the Klyuchevskoy volcano. The simultaneous activation of Bezymianny at the time of the eruption onset of Klyuchevskoy is obvious before the Tolbachik eruption (a1) but absent afterwards (a2).Similarly, at Klyuchevskoy, the precursory activity is evident for the eruptions that occurred before the Tolbachik eruptions (b1), but not afterward (b2). Time series generated using MATLAB software (www.mathworks.com).

Figure 6

 6 Figure 6 Reactivation of Tolbachik and responses at other volcanoes. -(a) Stacked TADR time-series for Klyuchevskoy, Bezymianny, and Tolbachik between 2011 and 2014. (b) Shallow LP seismicity (normalized number of events per day) below the Klyuchevskoy, Bezymianny, and Tolbachik during the same period (modified from Shapiro et al., 2017); (c) Deep LP seismicity below the KVG (seismic data modified from Shapiro et al., 2017); (refer to the color version of the figure). Time series generated using MATLAB software (www.mathworks.com).

Figure 7 -

 7 Figure 7 -(a) Stacked TADR time-series for Klyuchevskoy, Bezymianny, and Tolbachik volcanoes between 2011 and 2018. Eruptive periods for each volcano are labeled according to Table 1. Gray fields outline periods of Bezymianny activity characterized by frequent explosive eruptions (stars); (b) The reactivation of Bezymianny started in early 2016, through the extrusion of "cold" crystallized plug, undetected by MODIS. This anomalous precursory phase culminated in December 2016 (#BEZ23) with the effusion of a new lava flow that marks the resumption of activity of Bezymianny after the Tolbachik eruption. Time series generated using MATLAB software (www.mathworks.com).

Figure 8 -

 8 Figure 8 -(Left Axis) Ratio between the averaged Bezymianny TADR values within T days after and before the Klyuchevskoy eruptions (#1 to 6). (Right axis). Mean fraction in randomized data that exceed the observed ratio (see the text for more explanations).

Table 1a -

 1a Eruptions and parameters of Klyuchevskoy volcano retrieved from MODIS data. Volumes, TADR and MOR are calculated as bulk values (see Methods) Eruption #6 was detected by MODIS but with TADR < 0.25 m 3 s -1 . Eruption start and end, are selected manually based on first and last thermal anomalies. ** Start of Eruption#10 from 36

	4	15/10/2008	30/11/2008	07/01/2009	86.20	485.3	45.2	38.2	6.1
		16:20	10:50	15:55
	5	20/09/2009	28/08/2010	28/10/2010	441.60 255.8	151.5	54.1	4.0
		10:10	11:10	16:25
	6*	09/09/2012	04/11/2012	07/01/2013	119.80 681.97 0.15	0.1	0.01
		15:35	10:40	10:40
	7	13/08/2013	16/10/2013	25/10/2013	75.10	218.2	94.6	93.9	14.6
		15:10	11:05	11:00
	8	01/01/2015	14/01/2015	23/02/2015	55.80	433.0	12.3	10.2	2.6
		10:05	15:15	10:25
	9	06/04/2016	05/09/2016	01/11/2016	211.80 408.2	137.8	24.8	7.5
		15:20	14:25	16:00
	10**	11/11/2019	n.d.	n.d.	n.d.	1104.3	n.d.	n.d.	n.d.
						33
	erupti	start	peak	end	durati	dt_es	Vol	TADRma	MOR
	on				on			x
	#	dd/mm/yyyy	dd/mm/yyyy	dd/mm/yyyy	days	days	x 10 6 m 3 m 3 s -1	m 3 s -1
		hh:mm	hh:mm	hh:mm
	1	23/07/2003	07/11/2003	24/01/2004	186.90 N.d.	7.3	2.7	0.5
		16:10	10:10	16:05
	2	14/01/2005	07/02/2005	28/03/2005	73.60	356.0	52.3	40.7	8.2
		15:40	16:25	11:30
	3	18/03/2007	22/05/2007	18/06/2007	103.50 720.2	85.4	58.2	9.5
		15:35	10:35	10:15

*

Table 1b -

 1b Continue. Eruptions and parameters of Bezymianny retrieved from MODIS data. Volumes, TADR and MOR are calculated as bulk values (see Methods) Date of explosion from Ozerov et al. 38 **Event #15 not detected by MODIS but reported in Ozerov et al. 38

	14	14/10/	16/09/200	15/10/200	05/11/200	50.6	155.	101.	2.3	1.7	0.5
		2007	7 08:57	7 14:40	7 22:31		9	0
	15*	05/11/	N.d.	05/11/200	N.d.	N.d.	N.d.	N.d.	N.d.	N.d.	N.d.
		2007		7
	16	19/08/	02/08/200	11/08/200	06/09/200	35.4	300.	270.	1.8	3.7	0.6
		2008	8 08:08	8 11:45	8 17:04		9	4
	17	16/12/	03/12/200	17/12/200	07/01/201	35.7	493.	452.	5.8	8.2	1.9
		2009	9 00:28	9 15:20	0 16:41		1	3
	18	16/02/	26/01/201	09/02/201	27/02/201	31.3	54.0 19.1 2.3	2.7	0.9
		2010	0 18:46	0 14:40	0 01:27
	19	31/05/	10/05/201	29/05/201	15/06/201	36.2	109.	72.8 2.6	2.5	0.8
		2010	0 19:28	0 15:50	0 23:40		0
	20	13/04/	31/03/201	15/04/201	02/05/201	32.4	320.	288.	2.3	4.1	0.8
		2011	1 14:07	1 10:50	1 23:11		8	6
	21	08/03/	14/02/201	08/03/201	21/03/201	36.3	328.	287.	7.3	34.2	2.3
		2012	2 06:58	2 15:40	2 13:33		2	3
	22	01/09/	24/08/201	04/09/201	15/09/201	21.7	180.	156.
		2012	2 16:26	2 15:15	2 10:18		0	1
	even	Date	start	peak	end	durat	dt_p	dt_e	Vol	TADR	MOR
	t	Exp.*				ion	p	s		max
	#					days	days days x 10 6	m 3 s -1 m 3 s -1
									m 3
	1	13/03/	04/03/200	18/03/200	24/04/200	51.1	nan	nan	3.5	4.2	0.8
		2000	0 05:30	0 11:35	0 08:01
	2	01/11/	27/09/200	01/11/200	16/11/200	50.1	228.	156.	14.2 51.9	3.3
		2000	0 15:42	0 11:05	0 18:10		0	3
	3	06/08/	22/07/200	06/08/200	23/08/200	32.5	278.	247.	10.9 30.9	3.9
		2001	1 09:52	1 10:25	1 21:20	18.2	0 272.	7 253.	0.6	0.4	0.4
	4	16/12/	28/11/200	13/12/200	01/01/200	34.3	129. 0	96.8 8.1 1	20.7	2.7
		2001	1 17:20	1 11:05	2 23:58		0
	5	25/12/	12/12/200	25/12/200	07/01/200	26.0	377.	344.	4.8	50.9	2.1
		2002	2 20:43	2 14:55	3 20:17		2	9
	6	26/07/	13/07/200	26/07/200	10/08/200	28.7	213.	186.	5.3	16.1	2.1
		2003	3 00:25	3 15:15	3 16:45		0	2
	7	13/01/	06/01/200	14/01/200	23/01/200	17.2	171.	149.	0.6	3.2	0.4
		2004	4 17:04	4 10:00	4 22:43		8	0
	8	18/06/	03/06/200	19/06/200	30/06/200	26.7	157.	132.	1.6	2.0	0.7
		2004	4 22:01	4 15:10	4 13:43		2	0
	9	11/01/	10/01/200	16/01/200	20/01/200	10.4	211.	194.	0.3	0.4	0.3
		2005	5 13:02	5 15:40	5 22:28		0	0
	10	30/11/	16/11/200	30/11/200	16/12/200	29.4	317.	299.	6.4	12.8	2.5
		2005	5 15:08	5 11:50	5 00:15		8	7
	11	09/05/	29/04/200	08/05/200	27/05/200	28.7	159.	134.	1.1	0.9	0.5
		2006	6 07:39	6 11:05	6 23:59		0	3
	12	24/12/	15/12/200	23/12/200	03/01/200	18.7	229.	201.	0.8	1.5	0.5
		2006	6 18:15	6 11:25	7 10:45		0	8
	13	11/05/	29/04/200	12/05/200	07/06/200	38.8	140.	116.	6.1	18.5	1.8
		2007	7 13:39	7 15:55	7 08:13		2	1

*

Acknowledgments

This research was funded by the Italian Ministry for Universities and Research (MIUR). Part of this study was supported by the European Research Council under the European Union Horizon 2020 research and innovation program (Grant Agreement 787399-SEISMAZE) and by the Russian Ministry of Education and Science (Grant 14.W03.31.0033). We acknowledge the LANCE-MODIS system (http://lancemodis.eosdis.nasa.gov/) for providing Level 1B MODIS data.

Data Availability

The satellite datasets are available as Supplementary Material at:

Author Contributions

D.C. designed the study and conducted the analysis. M.L. and F.M. contributed to the analysis of MODIS data. N.S. provided the seismic data. S.H. conducted the statistical analysis. A.S., R.M., and T. W. reconstructed the eruptive history of volcanoes and interpreted satellite data in the light of independent observations. The ideas in the manuscript were developed through group discussions and written up principally by D.C. with contributions from all authors.

Additional Information

Supplementary information accompanies this paper:

(1) Supplementary_Material_KVG_Coppola.doc;

(2) Supplementary_TableS1.xls

Competing Interests

The authors declare no competing interests.