# Correlated petrology and seismicity indicate rapid magma accumulation prior to eruption of Kizimen volcano, Kamchatka

Lea Ostorero, Hélène Balcone-Boissard, Georges Boudon, Nikolai Shapiro, Alexander Belousov, Marina Belousova, Andreas Auer, Sergey Senyukov, Svetlana Ya. Droznina

► To cite this version:

Lea Ostorero, Hélène Balcone-Boissard, Georges Boudon, Nikolai Shapiro, Alexander Belousov, et al. Correlated petrology and seismicity indicate rapid magma accumulation prior to eruption of Kizimen volcano, Kamchatka. Communications Earth & Environment, 2022, 3 (1), pp.290. 10.1038/s43247-022-00622-3. hal-03862228

HAL Id: hal-03862228

https://hal.science/hal-03862228

Submitted on 23 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Correlated petrology and seismicity indicate rapid magma accumulation prior to eruption of Kizimen volcano, Kamchatka

Lea Ostorero<sup>D 1✉</sup>, Hélène Balcone-Boissard<sup>2</sup>, Georges Boudon<sup>1</sup>, Nikolai M. Shapiro<sup>D 3,4</sup>, Alexander Belousov<sup>D 5</sup>, Marina Belousova<sup>5</sup>, Andreas Auer<sup>D 6</sup>, Sergey L. Senyukov<sup>D 7</sup> & Svetlana Ya. Droznina<sup>7</sup>

The increase in number and intensity of earthquakes during a pre-eruptive crisis is the main basis of seismic volcano monitoring. However, a strong understanding of how these seismic signals relate to magmatic processes in the magma plumbing systems prior to volcanic eruptions is crucial for these efforts. Here we compare the characteristics of a seismo-volcanic crisis prior to the 2010-2013 explosive-extrusive eruption of Kizimen volcano, Kamchatka with the timescales of processes in the magma plumbing system. These time-scales are inferred from the numerical modelling of iron-magnesium intracrystalline inter-diffusion in 88 zoned orthopyroxene crystals from dacites and silica-rich andesites collected after the eruption. We find that the eruptible magmas were assembled rapidly during a magma mixing process beginning around 1.5 years before the eruption, which is well correlated with the onset of the seismic crisis. We conclude that the observed seismic re-activation marked the onset of magma mixing and led to destabilization of the reservoir, followed by the eruption.

<footnote n="1"><sup>1</sup>Université Paris Cité, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France. <sup>2</sup>Institut des Sciences de la Terre de Paris (ISTeP), UMR 7193, CNRS-Sorbonne Université, Paris, France. <sup>3</sup>Institut des Sciences de la Terre (ISTerre), Université Grenoble Alpes, CNRS, 38400 Grenoble, France. <sup>4</sup>Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, 10 B. Gruzinskaya St., Moscow 123995, Russia. <sup>5</sup>Institute of Volcanology and Seismology, 9 Piip Boulevard, Petropavlovsk-Kamchatsky 683006, Russia. <sup>6</sup>Department of Geoscience, Shimane University, 1060 Nishikawatsu, Matsue 690-8504, Japan.</footnote>

<footnote n="7"><sup>7</sup>Kamchatka Branch of the Geophysical Survey, Russian Academy of Sciences, 9 Piip Boulevard, Petropavlovsk-Kamchatsky 683006, Russia.</footnote>

<footnote>✉email: ostorero@ipgp.fr</footnote>
During volcanic unrest, the information about the magmatic processes eventually leading to an eruption can be obtained from geochemical, geodetic, and seismic observations¹⁻⁴. The latter are particularly valuable because of their depth and time resolution. An increasing number and intensity of volcanic earthquakes is considered the hallmark of an impending eruption⁵. While the correlation between the preparation of eruptions and the increased seismic activity is well-established empirically, relating the latter to specific magmatic processes remains problematic. At the same time, the temporal evolution of magmatic processes in the plumbing system can be deciphered based on modern petrological and geochemical methods⁶⁻⁹. Therefore, a more direct link between seismicity and magmatic process can be obtained by comparing petrological to seismological records.

Within a reservoir, crystals grow from melt incorporating major and trace elements following thermodynamic and kinetic laws¹⁰. In silica-rich melt, orthopyroxenes (Mg, Fe)₂Si₂O₆ are relevant candidates as their composition evolve from Mg-rich to Fe-rich composition, depending especially on melt temperature and composition, creating zoned crystals⁶. Diffusion in the zoned crystals is a time-dependent process (Fick's law) acting as a clock for measuring timescales of a variety of magmatic processes⁷, 8, 11, 12. Zoned orthopyroxenes in evolved mushy zones below volcanoes are thus time capsules allowing us to constrain the pre-eruptive timescales of magmatic processes⁷, 8, 11, 12.

Timescales modelled using (Fe–Mg) interdiffusion chronometry applied to well-monitored eruptions, may be correlated to monitoring data (seismicity, changes in the gas chemistry and/or flux, deformation) and then have implications for volcanic eruption forecasting models⁷, 13–18 and a high potential for petrological monitoring of volcanoes¹⁹. Up to now only syn-eruptive investigations of magmatic processes in correlation with seismic activity have been performed for differentiated systems¹⁷, 20, 21 compared to pre-eruptive correlations. To establish a temporal correlation between the magmatic processes and the seismic activity during a pre-eruptive re-activation of an arc plumbing system, we selected products from the first historical magmatic eruption of Kizimen in 2010–2013, in Kamchatka, involving differentiated magmas. This eruption is of particular interest because it occurred after several decades of quiescence and was well-monitored²²–²⁵. Seismic network operated by the Kamchatka Branch of the Geophysical Survey of the Russian Academy of Sciences resulted in a strong record of the pre-eruptive seismic re-activation. Furthermore, the petrological record in the orthopyroxene crystals that is very sensitive to temperature fluctuations²⁰, 26–29 has not been perturbed by previous volcanic activity.

Therefore, this 2010–2013 Kizimen eruption provides us with a unique opportunity to study magmatic-seismic correlation and better understand the processes leading to intermediate to silicic eruptions. We report a good correlation between the petrological timescales and the pre-eruptive seismic monitoring 18 months prior the eruption. The petrological timescales are thus a valuable warning-clock for volcano monitoring, particularly in subduction settings.

**Geological context and previous results about the 2010 eruption of Kizimen.** The Kamchatka Volcanic Arc, in Eastern Russia, results from the subduction of the northeast corner of the Pacific plate at a rate of 8–9 cm/yr beneath the Okhotsk plate (which is connected to the North American plate)³⁰, 31 (Fig. 1a). This arc is the most productive arc on Earth (~310 millions of tons of magma emitted per year from the late Pleistocene to the Holocene)³², with more than 25 active volcanic centers. One of

the world's largest clusters of subduction volcanoes, the Klyuchevskoy volcanic group, is located in the northern part of the arc, in the Central Kamchatka Depression³³, 34 where the largest magmatic production is recorded²², 31, 35 and where the volcanoes exhibit very sustained seismic activity (Fig. 1a, b). The particularly strong volcanic activity in this region is likely caused by the hot mantle material brought by the around-slab-edge asthenospheric flow³⁶ and related crustal extension³⁷. Nearly continuous occurrence of deep long-period earthquakes near the crust-mantle boundary³⁵, 38 is likely related to the underplating of mafic magmas and widespread seismic tremors mark an active transcrustal magmatic system connecting different volcanoes of this group³⁹. At the same time, well-developed magma chambers beneath some of most active volcanoes such as Bezymianny are well documented⁴⁰.

Kizimen is located on the eastern margin of the Central Kamchatka Depression and its magmas might have similar origin to those from the Klyuchevskoy volcanic group. Because of the less dense coverage by seismic stations, the deep structure of the crust beneath Kizimen is not well known and the reported magmatic and seismic activity is linked to the upper crust. At the same time, the existence of crustal scale magmatic system in analogy with the Klyuchevskoy volcanic group cannot be excluded. Kizimen is a Late Pleistocene - Holocene volcano, with four cycles of eruptive activity identified (between 12,000–11,000 to 8400 years ago, 8400–6400, 6400–3000, and 3000 years ago to the present³⁴).

A possible phreatic eruption occurred in 1927–1928³⁴, 41 but the first historical magmatic eruption of Kizimen occurred in 2010–2013, when magma volumes of 0.4 km³ Dense Rock Equivalent (DRE) were emitted⁴¹ (Fig. 1b; Supplementary Fig. 1). On 11 November 2010, explosive episodes generated pyroclastic density currents²², 41, 42, that were then followed by the extrusion of a lava dome which gradually formed a ~200 m thick lava flow⁴¹ (Fig. 1b; Supplementary Fig. 1). Numerous block and ash pyroclastic density currents were produced as a result of the gravity collapse of the lava dome and flow front. Dacites, silica-rich andesites and banded lavas, formed by an incomplete magma mixing event between two magmas were emitted⁴¹ (Fig. 1b–f; Supplementary Fig. 1 and Supplementary Data 1). For this eruption, a model of the pre-eruptive dynamics proposed that the emitted dacite came from the differentiation of a basaltic magma⁴¹. This basaltic magma was not erupted at the surface⁴¹. Several models have assessed the presence of a stable silicic reservoir under Kizimen, periodically refilled by magmatic injections⁴¹–⁴³. This silicic magma chamber is located between 5 and 11 km depth (based on earthquakes depths and amphibole and orthopyroxene compositions⁴¹, which is in agreement with a study on older eruptions from Kizimen with storage pressures of magmas between 1 to 3 kbars (up to 10 km depth)⁴², 43). Another seismic zone from the ground surface to ~5 km has been recognized during the 2010–2013 eruption and probably shows the conduit in which magma ascended to the surface⁴¹. An intensive seismic swarm²², 41 was recorded in 1963, with earthquakes depths reaching down to more than 35 km but without leading to an eruption⁴¹. This seismic swarm is interpreted as a mafic injection from depth, as the samples from the 2010 eruption that followed showed resorbed olivines and as it occurs in a similar geologic and tectonic situation at the volcano Tolbachik in the Klyuchevskoy group, where intrusion of dikes of mafic magmas generate multiple earthquakes at depths of up to 30 km⁴¹. An episodic increase of the fumarolic activity also followed the occurrence of large shallow earthquakes of this seismic swarm (magnitudes of 6.2 and 5.8)³⁴. However, due to the high viscosity of magma in the reservoir, a long time period can exist between the input of basic magma and the following
<figure><img src="image_3.png" /><figcaption>**Fig. 1** From macroscale to microscale: Kizimen location in Kamchatka, 2010–2013 eruption samples and orthopyroxenes zonations. **a** Kamchatka map with the three volcanic zones: Sredinny Range (SR), Central Kamchatka Depression (CKD) and Eastern Volcanic Front (EVF)³³,³⁴ and the boundary of the active subduction zone of the Kurile-Kamchatka Trench (subduction rate of 8 cm/yr³⁰,¹²²,¹²³). Image from the ArticDEM: Hillshade Gray from the ArcGIS Living Atlas of the World (by Esri Imagery, Digital Elevation Model created by the Polar Geospatial Center from Maxar Imagery; https://livingatlas2.arcgis.com/arcticdemexplorer/)¹²⁴ and modified according to other studies¹²⁵,¹²⁶. **b** Sampling sites on Kizimen (Supplementary Table 1) (Image credit: U.S. Geological Survey (USGS)/photo by Robert Simmon/Nasa Earth Observatory image from 17 June 2014, using Landsat 8 data from the USGS Earth Explorer). **c** Dacite (64 wt% SiO₂; 5.6 wt% Na₂O + K₂O). **d** Andesite (58–60 wt% SiO₂; 4.8–5 wt% Na₂O + K₂O). **e, f** Banded samples (alternation between dacitic and andesitic bands⁴¹). On the field, estimates of the percentages of these different samples were ~70% of banded lavas, ~10% of dacites and ~20% of andesites (even if it is unclear how the initial products looked like or what fraction remained in the chamber). **g, h** Backscattered electron images (BSE) of major zonations identified in the dacite and andesite: normal and reverse zonations, respectively. **i** Red-Green-Blue (RGB) chemical map of a normal-zoned orthopyroxene, with a Mg-rich core and Fe-rich rims ("Outer" on (i, j)). Mg is in blue and Fe in red. **j** Reverse-zoned orthopyroxene, with a Fe-rich core and Mg-rich rims, with the same colors as in (i).</figcaption></figure>

eruption⁴¹. This mafic injection led to magma hybridization of the lower part of the dacitic reservoir, forming an andesitic magma at the base of the magmatic reservoir⁴¹. A large portion of the dacitic magma was not affected by this mafic recharge, leading to the formation of a zoned reservoir.

A seismic crisis started ~1.5 years before the eruption with earthquakes located between 11 km and the surface²²,²³,⁴¹,⁴⁴. In addition, deformation of the volcanic edifice was detected by Interferometric synthetic aperture radar (InSAR)²⁴. The seismi data and deformation by InSAR were interpreted as the injection of a dike of 15 km long and about 20 cm thick in the magmatic reservoir²⁴. As no seismicity below 12 km has been recorded, this

interpretation is thought not to fit the existing geological and seismological data⁴¹. Prior to the eruption of 2010, another mixing event occurred between the dacitic and andesitic part of the reservoir due to some degassing of the dacitic magma leading to magma convection, with incomplete mixing between the two magmas⁴¹. Kizimen is crossed by an active tectonic fault, through which volatiles are escaping to the surface. This leakage is visible as a persistent fumarole on the NW flank of the volcano⁴¹, this long-term passive degassing could be responsible for the depressurization of the magma chamber and increase the viscosity of the magma⁴⁵,⁴⁶ (Fig. 1b). The temperatures of the fumaroles at the summit of Kizimen have been measured from 15 to 31 August 2009
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Enstatite (En) content histograms of orthopyroxenes of the 2010–2013 eruption. **a** Unzoned and **b** single-zoned orthopyroxenes (unzoned orthopyroxenes compositions shown in the dashed frame). The arrows underline the major path from the core to the rim recorded by the single-zoned orthopyroxenes. The number of orthopyroxenes analyzed is specified (n). The cores of the zoned orthopyroxenes are in the same range as the unzoned orthopyroxenes. Orthopyroxene is noted as opx.</figcaption></figure>

and a temperature of 340 °C was reached close to the central mouth of the outermost fumaroles on Kizimen which is higher than measurements prior to the onset of the unrest (maximum of ~270 °C)⁴⁷. However, this temperature increase of the fumaroles is not thought to be a consequence of seismic activation⁴⁷. One month before the eruption, visual signs of reactivation were observed with thermal areas and two new large fumaroles that appeared (16 October 2010)²²,⁴¹.

## Results

**Magma compositions and orthopyroxene crystals.** Magmas are differentiated, with whole rock compositions in the andesitic and dacitic field (Fig. 1c, d; Supplementary Fig. 1; Supplementary Note 1 and Supplementary Data 1). Orthopyroxenes are present in the dacite and andesite. 70% of orthopyroxenes are unzoned (Supplementary Figs. 2, 3). In the four fractions studied (see “Methods”), the proportions of unzoned orthopyroxenes are approximately the same (Supplementary Fig. 3). Zoned crystals were classified according to the type of zonation observed: either a single one (normal or reverse zonation) or a multiple one (Fig. 1g–j; Supplementary Figs. 2, 3; Supplementary Note 2). In the dacite, the majority of the zoned orthopyroxenes are normal-zoned while, in the andesite, reverse-zoned orthopyroxenes are more abundant (Supplementary Figs. 2, 3). Few orthopyroxenes show reactional borders or remelting textures (Supplementary Figs. 2, 3). Out of the 1841 orthopyroxene crystals mounted, 112 zoned orthopyroxenes compositions and 108 unzoned crystals compositions were investigated in the dacite and andesite (Fig. 2; Supplementary Data 2; Supplementary Figs. 4, 5).

Orthopyroxenes compositions are spread over an En$_{44–73}$ compositional range (Enstatite content (En = $\frac{Mg}{Mg+Fe}$)). Three compositional groups are identified: the first one corresponds to all of the cores of the orthopyroxenes (both unzoned and zoned orthopyroxenes of the dacite and andesite; En$_{62–64}$), then a group for the normal-zoned rims of the orthopyroxenes of the dacite (En$_{44–54}$) and finally a group for the reverse-zoned rims of the orthopyroxenes in the andesite (En$_{64–73}$) (Fig. 2; Supplementary Figs. 4, 5 and Supplementary Data 2; Supplementary Note 3). Harker diagrams also highlight that the unzoned orthopyroxenes and cores of the zoned orthopyroxenes have the same compositions (Supplementary Fig. 5).

**Timescales.** Out of 112 profiles on zoned crystals, 26 were suitable for diffusion modelling for the dacite, as well as 62 for the andesite (Supplementary Figs. 6–8; Supplementary Note 4 and Supplementary Data 2). The temperatures used to model the timescales are calculated on magnetite/ilmenite couples in the dacite and andesite, this equilibrium between these minerals has been widely used to estimate temperatures in igneous rocks⁴⁸. Mean temperatures of 832 ± 37 and 850 ± 57 °C are found for the andesite and dacite respectively⁴⁸ (Supplementary Note 5, Supplementary Fig. 6, and Supplementary Data 3, 4). These profiles have sigmoidal shapes due to diffusion effects, without bumps or peaks that could indicate a growth component²⁶ (Supplementary Figs. 7, 8 and Supplementary Data 2). Al$_{2}$O$_{3}$ and CaO gradients are also compared to the En content to verify that the gradients in En are due to diffusion⁴⁹ (Supplementary Fig. 7; Supplementary Data 2). For both the dacite and andesite, most of the timescales
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Linking monitoring data with timescales of magmatic processes. **a** Time vs depth plot of seismicity at Kizimen from 2003 to the onset of the eruption (11/11/2010)²²,⁴¹. The onset of the surface deformation detected by InSAR is specified²⁴,⁴¹. **b** Cumulative number of the timescales estimated on the major compositional changes evidenced for the dacite and andesite (normal zoning for the dacite (En₆₂⁻₆₄ to En₄₄⁻₅₄) and reverse zoning for the andesite (En₆₂⁻₆₄ to En₆₄⁻₇₃)); Supplementary Data 4) and cumulative seismic moment. **c, d** Focus on the 1.5 years pre-eruption. More crystal clocks are started closer in time to the eruption (steeper slope one year before the eruption). **e** Uncertainties associated to the individual timescales modelled shown on (c) (Supplementary Fig. 6 and Supplementary Data 4). These uncertainties are due to the diffusivity calculations (D₀, the interdiffusion coefficient, the activation energy Eₐ), temperature and resolution of grayscale values of the BSE images⁶,¹¹⁹. They are calculated using a Monte Carlo simulation. These uncertainties are asymmetric, with a larger positive error bar and a smaller negative error bar, due to the logarithmic effect of the uncertainties on mainly temperature and diffusivity calculations²⁶. Orthopyroxene is noted as opx.</figcaption></figure>

are young (<2 years) for single-zoned and outer rims of multiple-zoned orthopyroxenes (Supplementary Fig. 6, Supplementary Note 5, and Supplementary Data 4). Longer timescales associated to multiple zonations can be related to previous events in the magmatic reservoir, recorded by the orthopyroxenes. On the major En changes recorded by the zoned orthopyroxenes (En₆₂⁻₆₄ to En₆₄⁻₇₃ in the reverse zoned orthopyroxenes of the andesite and from En₆₂⁻₆₄ to En₄₄⁻₅₄ in the normal-zoned crystals of the dacite), 80% of the timescales modelled are in the 1.5-years period prior to the eruption whatever the type of zoning (normal or reverse) (Fig. 3b, d; Supplementary Figs. 6, 8 and Supplementary Data 4).

## Discussion

Magmatic processes modify the storage conditions of magmas before volcanic eruptions (changes in temperature, pressure, oxygen fugacity (fO₂) or volatile content)⁷,¹³,¹⁵,²⁶,⁵⁰⁻⁵² and Fe-Mg interdiffusion timescales modelling of the zonations in the orthopyroxenes help us to retrieve the timescales of these

processes by studying the start of the “clocks” of these crystals (Fig. 3b, d; Supplementary Fig. 6 and Supplementary Data 4).

70% of the orthopyroxenes are unzoned (Supplementary Figs. 2, 3). Different hypotheses can be evoked to explain this high proportion of unzoned orthopyroxenes: (i) the magmatic changes did not affect the entire magmatic reservoir; (ii) the unzoned crystals were older and their zoning patterns were homogenized by diffusion with time; (iii) they were formed after the last magmatic perturbation discussed below. The size of the unzoned crystals and their proportions compared to zoned orthopyroxenes are comparable in the different fractions studied (from 125 to 710 µm; Supplementary Figs. 2, 3). Their compositions are also in the same range as the cores of the zoned orthopyroxenes, in En₆₂⁻₆₄ and in the same domains in the CaO, FeO, MgO and MnO versus SiO₂ diagrams, compared to the rims compositions of the zoned orthopyroxenes (Fig. 2; Supplementary Figs. 4, 5). It was then more likely that they formed before the last magmatic event prior to the eruption and that the magmatic changes did not affect the entire magmatic reservoir (Fig. 3a).
<figure><img src="image_6.png" /></figure>

**Fig. 4** Pre-eruptive dynamics of the magmas in the timeframe given by the timescales modelled in zoned orthopyroxenes: crystal clocks mainly starting 1.5 years before the eruption. a Magma plumbing system below Kizimen before the 2010 eruption. The architecture of the plumbing system of a previous study was used⁴¹: a dacitic reservoir which has been intruded by a basaltic injection, possibly in 1963²²,²⁵,⁴¹, with the lower part of the dacitic reservoir becoming andesitic by hybridization with the mafic magma. During this hybridization, crystals of the dacite get mixed with the andesite (cores of all orthopyroxenes are the same as the dacite) (Fig. 2; Supplementary Figs. 4, 5). Then, the most likely hypothesis to explain orthopyroxenes zonations in the dacite and andesite is presented (see “Discussion”): b Magma mixing occurred 1.5 years before the eruption, at the contact between the dacite and the andesite, exchanging orthopyroxenes from the dacite with the andesite (explaining the formations of reverse-zoned orthopyroxenes in the andesite by increase of the temperature (T)). Orthopyroxenes crystals can also return in the dacitic magma, forming normal zonations, and get reincorporated in the andesite before the eruption (forming multiple-zoned crystals). As a result of mixing of magma with colder areas in the dacite, the orthopyroxenes in the dacite mainly record normal zonations. Not all the reservoir was affected by the mixing process as there is a majority of unzoned orthopyroxenes. Seismicity is recorded during this period. c Magma mingling shortly prior eruption. The rate of crystal diffusion starting rates (progressive increase of timescales from April 2009 to the eruption) and progressive increase of seismic moment and number of earthquakes are shown with the shading in the arrows. The deeper part of the transcrustal system may be composed of mafic lenses. The dacitic and andesitic reservoir have the same mineral assemblage: plagioclase, orthopyroxenes, amphiboles, Fe-Ti oxides and some resorbed olivines and quartz (T = 832 ± 37 °C and 850 ± 57 °C, respectively).

As all the unzoned and zoned orthopyroxenes have the same core compositions (Supplementary Figs. 4, 5 and Supplementary Data 2), the orthopyroxenes from the dacite and the andesite are likely to be originated from the dacitic magma present in the reservoir before the 1963 basaltic injection. The andesitic part of the reservoir was then probably formed after the injection of an aphyric basaltic magma in 1963 and mixed with the dacite in the lower part of the reservoir, as suggested in a previous study of the eruption⁴¹ and incorporated some of the crystals of the dacite (including orthopyroxenes) (Fig. 4a). This is also confirmed by a previous study on plagioclase crystals in banded rocks of the same eruption⁴¹.

The majority of the zoned orthopyroxenes studied indicate a short timescale (<1.5 years) and no zoned orthopyroxene reflects a longer timescale in the order of 50 years (that could correspond to a possible mafic injection in 1963) (Supplementary Fig. 6).

The uncertainties associated to the timescales modelled are mainly due to the temperature uncertainty and uncertainties due to diffusion coefficient determination⁶,²⁶,⁵³ (see “Methods”; Fig. 2; Supplementary Fig. 6; Supplementary Data 2). For most of the short timescales (<2 years), the positive error associated to the timescale is in the order of 1 year and the negative error is around 100 days (Supplementary Fig. 6; Supplementary Data 4). The young timescales have lower uncertainties than the timescales superior to 1 year before the eruption (Fig. 3e).

The temperatures used to model the timescales were determined using the compositions of magnetite/ilmenite pairs which can show the variability of the temperatures in the reservoir.

Varying the temperature and using the one determined in orthopyroxenes in the andesitic part of mingled samples by another study⁴¹ implies shorter timescales for the andesite (divided by 2, with T = 875 °C), showing the importance of the choice of the temperature in the timescale modelling. As no temperature was determined in orthopyroxenes of the dacite in the previous study⁴¹ and no melt inclusion with a sufficient size was found in our samples, the temperatures determined on the magnetite/ilmenite pairs were used.

On Fig. 3, considering only the major En changes of the zoned orthopyroxenes (normal zoning for the dacite (from En₆₂⁻₆₄ to En₄₄⁻₅₄) and reverse zoning for the andesite (from En₆₂⁻₆₄ to En₆₄⁻₇₃); Supplementary Data 4), all the zoned orthopyroxenes give timescales equal to or below 1.5 years. The cumulative number of timescales recorded mostly increase 1 year before the eruption (Fig. 3d). So, 1.5 years prior to the eruption of 2010, after ~50 years of dormancy, the orthopyroxenes in the dacitic and andesitic magmas begin to record a progressive magmatic disturbance accompanied by changes in the storage conditions, with different processes as their most external zonations are not the same (Fig. 3b, d; Supplementary Figs. 2–4). One year before the eruption, these magmatic processes become more dynamic (as the curve of the cumulative timescales is steeper than before; Fig. 3d). In the andesite, reverse zonations (~80%) in the zoned orthopyroxenes reflect heating processes²⁶,²⁷ or degassing⁵⁴ and normal zonations (~20%) reflect cooling processes¹⁵,¹⁷,²⁶,²⁷ or a high water content in the melt⁵⁴ (Fig. 2; Supplementary Figs. 2–4). In the dacite, dominant normal zonations (93%)
principally register cooling processes¹⁵,¹⁷,²⁶ or a high water content in the melt⁵⁴ (Fig. 2; Supplementary Figs. 2–4).

To explain the reversed-zoned rims of the andesite, four hypotheses can be evoked:

i. the injection of a mafic magma at the base of the reservoir, 1.5 years prior to the eruption. This mafic magma would have been aphyric as no orthopyroxenes have a different En core composition. But as no deep earthquakes (>11 km) have preceded and accompanied the 2010 eruption, this event is not likely to have been triggered by input of basic magma from greater depth⁴¹, unless magma ascent occurred through previous pathways, explaining the lack of deep seismicity. But if the last mafic injection occurred in 1963, it is likely that fracturing of the rocks needs to happen for the ascending magma, after ~50 years, generating seismicity. No enclave of basaltic or basaltic andesite compositions were found, unlike many of the earlier eruption products where mafic recharge before the previous eruptions have been inferred⁴²,⁴³. Furthermore, the melt inclusions data are in equilibrium with their hosting orthopyroxenes (Supplementary Fig. 9) and they are in the same domain as the residual glasses (Supplementary Fig. 10; Supplementary Note 6). The melt inclusions do not show different CaO or Al₂O₃ values compared to the En content compositions (Supplementary Fig. 10). So, either the mafic magma only heated the andesite without mixing with it (which is unlikely, with the chaotic nature of mixing processes⁵⁵⁻⁵⁷), generating 30% of zoned orthopyroxenes or it mixed in a part of the reservoir without reaching the surface during the eruption.

ii. A convective event in the andesitic part of the reservoir, zoned in temperature, between the walls and cores of the reservoir, generating the reverse-zoned rims in the andesite. However, this convective magma mixing process is unlikely to occur without affecting the dacite, as only one reservoir, or a reservoir with linked magma pockets, is more likely to be present at Kizimen (depths of 5–11 km⁴¹; locations of earthquakes and amphibole and orthopyroxenes compositions), even if magma pockets can remain isolated.

iii. A partial mixing in the reservoir at the contact between the andesitic and dacitic magmas, with crystal exchanges. Some orthopyroxenes from the dacite could have migrated in the hotter andesite, generating reverse zonations in the orthopyroxenes (rims of En₆₂⁻₆₄) (Fig. 4b). In parallel, some orthopyroxenes from the andesite could have been in contact with the dacitic magma generating normal zonations (rims of En₅₅⁻₆₁) (Fig. 4b). This mixing process only affected part of the two magmas in contact, explaining that 70% of the orthopyroxenes in the andesite and dacite are unzoned.

iv. Degassing processes in the andesite, with a decrease of the water content in the melt, which favors Mg-rich orthopyroxenes (producing reverse zoning)⁵⁴. Indeed, series of H₂O fluid-saturated phase-equilibrium experiments showed that the Fe²⁺-Mg exchange coefficient K_D, between a rhyolitic melt and minerals, such as orthopyroxenes, increases as a function of dissolved water in the melt⁵⁴. An increase of the Mg# of the orthopyroxenes was observed when the H₂O content in the melt decreases, which means that Mg-rich orthopyroxenes (reverse zoning) crystallizes at low water contents whereas Fe-rich orthopyroxenes (normal zoning) are favored at high water contents¹⁸,⁵⁴. Normal-zoned plagioclase (sodic rims) generally accompany reversely zoned orthopyroxenes in the case of loss of H₂O during degassing⁵⁴. However, plagioclase rims in the andesitic bands of banded samples of the same eruption are of

different types but calcium-rich rims are mainly found (reverse-zoned)⁴¹. Plagioclase in the andesites then need to be studied, to see if normal-zoned plagioclase crystals are dominant, alongside the reversely-zoned orthopyroxenes. The mean variation of Fe²⁺-MgK_D from the cores to the reverse-zoned rims of the orthopyroxenes, ΔFe²⁺-MgK_D, of the andesite was estimated at ~0.06 and was used to estimate the variation in H₂O content in the melt, from the core of the orthopyroxene to the rim. Using an experimental relationship¹⁸,⁵⁴, this ΔFe²⁺-MgK_D corresponds to a ΔH₂O₉₄⁻₉₅ of ~0.7 wt% H₂O, showing a decrease in water content in the melt when the orthopyroxenes rims formed (ΔH₂O₉₄⁻₉₅ = wt%H₂O₉₄⁻ - wt%H₂O₉₅) (see Supplementary Note 7). However, the initial H₂O content of the melt was not estimated, as the few melt inclusions studied here were too small to be analysed with the Secondary-ion mass spectrometry, so the initial dissolved water composition is not known and cannot be compared with the experiments, that studied the Fe²⁺-MgK_D variations as a function of dissolved water in the melt (from 2.7 to 5.6 wt %)⁵⁴. Still, we calculated an estimate of the H₂O content (2.4 wt% H₂O), using only one melt inclusion in a plagioclase from the dacite (plagioclase/liquid thermometer⁵⁸; and hygrometer⁵⁹ (Supplementary Note 6, 7 and Supplementary Data 1). Provided that all these estimates are correct and based on only one melt inclusion, the H₂O content decreased from 2.4 to 1.7 wt%, which is outside of the water contents investigated in the experimental study⁵⁴. We are then not able to conclude on this hypothesis, due to the scarcity of our melt inclusions data.

The hypothesis of magma mixing between the dacitic and andesitic magmas is then more likely, given the points discussed above. But if it explains the normal and reverse zonations in the orthopyroxenes of the andesite, it does not explain the normal zonations in the crystals of the dacite. The high viscosity of the dacitic reservoir probably prevented the magma mixing of all the dacite with the andesitic part of the reservoir⁴¹,⁶⁰⁻⁶⁶.

The normal zonations in the orthopyroxenes of the dacite can be explained by cooling¹⁵,¹⁷,²⁶,²⁷ or a high water content in the melt⁵⁴. As only 30% of the orthopyroxenes in the dacite are zoned and that normal zoning is dominant (93% of the zoned orthopyroxenes; Supplementary Fig. 2), an increase in the H₂O content in the melt by differentiation could have led to the development of normal zoning in the orthopyroxenes⁵⁴. Some movements of some orthopyroxenes of the dacite in colder areas in the dacitic part (for example with the walls of the dacitic reservoir, which are colder than the core), could also have led to the formation of the normal zonations (Fig. 4b). ΔFe²⁺-MgK_D between the cores and normal-zoned rims of the orthopyroxenes in the dacite is equal to ~-0.05, corresponding to a mean water content change of ~-0.5 wt% H₂O (see Supplementary Note 7). This means that there was a water content increase in the melt from the core to the rims of the orthopyroxenes (from 2.4 to ~3 wt% H₂O), which is in agreement with the experimental study, where normal zoning is favored at higher water content in the melt⁵⁴. However, in the light bands of the banded samples, normal-zoned plagioclase are dominant⁴¹, which is not coherent with the presence of normal-zoned orthopyroxenes in the dacite, according to the experimental study⁵⁴. Further investigation of the plagioclase in the dacite is then also needed. The hypothesis of some partial internal mixing within the dacitic pocket with colder areas is then preferred to explain the normal zonations in the dacite (Fig. 4b).

A period of unrest then follows, while the magmas reach an eruptible state⁴¹, with magma mixing processes still occurring, up to the time of the eruption.
Finally, magma mingling occurred between the andesitic and the lower dacitic part, during ascent shortly prior to the eruption, forming the banded andesites and dacites (Fig. 4c). Magma mingling events between different magma types are also recognized as very short events prior to an eruption, as shown combining magma mixing experiments, numerical simulations and statistical analysis before eruptions in the Campi Flegrei volcanic system in Italy, where mixing processes can be of the order of a few days to a few minutes<sup>57,67–69</sup>.

The eruption trigger could then correspond to the magma mixing process beginning 1.5 years before the eruption and increasing progressively until it reaches an eruptible state, with the last magma mingling phase and then the eruption. Another possibility for the eruption trigger could be a decompression of the dacite, generating the last magma mingling phase leading to the eruption. However, as no explosive products from the first phases of the eruption were found, it is difficult to conclude on a scenario for the eruption trigger.

CO₂ flushing could also participate in triggering the eruption, i.e. interaction between carbonic fluids and H₂O-rich magmas stored in the upper crust<sup>70</sup>. Due to fluid saturation and magma ascent, release of carbonic fluids can occur, increasing the magma volumes and generating overpressure that can destabilize the magmatic system, leading to the eruption<sup>70</sup>. A way to test this hypothesis is to study melt inclusions H₂O and CO₂ contents, however, as said before, no H₂O and CO₂ content measures were carried out in the samples, as the melt inclusions were too small (only one melt inclusion in a plagioclase gave a water content of 2.4 wt% H₂O; plagioclase/liquid thermometer<sup>58</sup>; and hygrometer<sup>59</sup>) (Supplementary Note 6, 7 and Supplementary Data 1). Data from the literature of the same eruption in the banded lavas give H₂O content between 5 and 5.7 wt% H₂O in the dacite and 4.9 and 5.6 wt% H₂O in the andesite (amphibole thermometers<sup>41,71,72</sup>) but no CO₂ content have been measured. From an older eruption, in 7550 <sup>14</sup>C BP, H₂O content are of 0.5–3.5 wt% in melt inclusions in plagioclase and magnetite, 0.6–4.5 wt% H₂O in amphiboles and 1.6–5 wt% in orthopyroxenes and CO₂ of 62–9 260 ppm (with 2 over 46 above 4000 ppm)<sup>73</sup>. As data from another eruption is not relevant to study the hypothesis of CO₂ flushing, we are not able to conclude on this hypothesis.

The distribution of timescales obtained on the orthopyroxenes crystals can be correlated with the seismic activity (Fig. 3). For this goal, we used the catalog of volcanic earthquakes produced by the Kamchatkan Branch of Russian Geophysical Service. The completeness magnitude for this catalog is difficult to estimate and can be variable in time depending on the number of stations and meteorological conditions. Therefore, instead of reporting the number of recorded earthquakes (that is strongly sensitive to this completeness magnitude), we use the cumulative seismic moment that is dominated by strongest earthquakes (Fig. 3a–d). Also, this quantity is closely related to the overall energy seismic release and therefore has a physical relationship with ongoing processes within a volcano. Because Kizimen remained quiescent for decades and is located far from other volcanoes, this catalog is exclusively related to the 2010 eruption (and to the magma recharge potentially in 1963 that did not lead to the eruption<sup>41,74</sup>). Comparison of the timescales modelled with the cumulative number of earthquakes is also shown (Supplementary Fig. 11), with an increase of the number of timescales that record magma mixing processes closer to the time of the eruption, even if the number of earthquakes is increasing faster than the rate of crystal clocks starting. This could be due to the increase of seismic unrest near the time of the eruption, with a relatively lower rate of magma mixing. An inflection of the cumulative seismic moment and number of timescales can be evidenced, from ~January

2010–November 2010 (approximately one year before the eruption; Fig. 3d; Supplementary Fig. 11). This inflection can correspond to an intensification of the magma mixing processes (between the dacite and the andesite). As this process becomes more dynamic, the orthopyroxenes rims are formed more and more rapidly.

80% of the modelled timescales cover the period between April 2009 and the onset of the eruption on 11 November 2010 and 60% from January 2010 to November 2010 (Fig. 3b, d). This increase in number of petrological timescales coincides with the onset of the seismic crisis recorded ~1.5 years before the eruption<sup>22,41</sup>, with earthquakes between ~3 and 13 km from April 2009 to January 2010 and then becoming more numerous and shallower from January 2010 to the eruption (Fig. 3a, c and Supplementary Fig. 8). This increase of the timescales is linked to magmatic reactivation, as discussed before, with magma mixing at the contact between the two magmas. Magma mixing probably generated an overpressure in the reservoir at the origin of the seismicity recorded during the same time period (Fig. 3a–d), April 2009–November 2010. As the timescales are also more numerous 1 year before the eruption, coinciding with the increase of the seismic moment around that same time (January 2010; Fig. 3c, d), increase in the magma mixing dynamics between the andesite and dacite could induce bubble nucleation (and overpressure), leading to an increase in the seismicity<sup>75</sup>. The cumulative number of timescales and cumulative seismic moment are perfectly overlaid around April 2010 (8 months before the eruption; Fig. 3c, d), probably because of the increase in the magmatic processes near the onset of the eruption. Finally, the very short seismic and timescales signals close to the eruption (in the order of hours to a few days; steeper curve in Fig. 3c, d) can correspond to the magma mingling event, prior to the eruption on 11 November 2010, as evidenced before several explosive eruptions<sup>76</sup>.

We note that the apparent migration of the seismicity around January 2010 becoming shallower (Fig. 3a, c), between –5 km to +3 km is likely due to the modification of the software used for the analysis of the earthquakes<sup>22</sup>. While this modification might affect the number of detected small earthquakes its influence on the cumulative seismic moment is negligible (Fig. 3b, d). Progressive ground deformation was also detected before the eruption, from September 2008 to September 2010 (InSAR data)<sup>24</sup>. These types of signals could be partly explained by the increase of fluid pore pressure along the plane of the tectonic fault<sup>41</sup> and be related to magma mixing and movements<sup>13,14</sup>. The increase of temperatures of the fumaroles measured in August 2009<sup>47</sup>, are compatible with the magmatic disruptive processes in the reservoir (Fig. 4).

The orthopyroxene chronometer is thus successfully used here to determine the links between pre-eruptive magma dynamics in the plumbing system and pre-eruptive monitoring data at the surface for an eruption involving differentiated magmas. This chronometer shows the reactivation of the magmatic feeding system of Kizimen 1.5 years before the eruption, well correlated with the monitoring signals recorded on the surface from April 2009 to the time of eruption. Most of the studies investigating petrology and monitoring signals are carried out in mafic systems, using olivines, such as Mt Etna (Italy)<sup>13,15</sup>, Kilauea (Hawai'i)<sup>77</sup>, Eyjafjallajökull (Iceland)<sup>14</sup>, Piton de la Fournaise (La Réunion)<sup>78,79</sup>, El Hierro (Canary Islands)<sup>2,80</sup>, and arc volcanoes such as Shishaldin volcano (Alaska)<sup>16</sup>, Merapi (Indonesia)<sup>81</sup> (using clinopyroxenes), Irazú (Costa Rica)<sup>82</sup> or Llaima (Chile)<sup>83</sup>. These studies obtained mixing-to-eruption or recharge events of the order of a few days to months (3–11 months<sup>2,14–16,78,80,83,84</sup>) to years before the eruption (1–2 yrs)<sup>13,77,81,82</sup>. Compared to monitoring signals, it depends on the eruption, but overall, the petrological timescales match well with seismic precursory
activity months to days before the eruption<sup>2,16,78,80,83</sup> or with the GPS displacements<sup>14</sup>. The timescales obtained on arc volcanoes such as Merapi or Irazú agree with the ones found on Kizimen (2 and 1.1 yrs in average, respectively), which are also consistent to monitoring signals (seismicity, deformation)<sup>81,82</sup>.

To the best of our knowledge, only three studies have used the orthopyroxene chronometer, but to investigate syn-eruptive monitoring<sup>17,20,21</sup>. On Mount St-Helens correlations were made with some eruptive phases during the 1980 eruption but no correlation was made with seismicity before the beginning of the eruption<sup>17</sup>. According to the diffusion timescales modelled, timescales of the order of 1 year before the eruption in 1980 were found<sup>17</sup>, without being correlated to the seismicity that began 2-3 months before the eruption. On Ruapehu (New Zealand)<sup>20</sup> as on Bezymianny (Kamchatka)<sup>21</sup>, the quasi-permanent high seismic activity of these volcanoes makes it difficult to clearly correlate timescales and seismicity, but magma perturbations occurred months to years before eruption. One recent study has used orthopyroxenes to decipher the processes and pre-eruptive conditions that led to the 1990 Kelud basaltic andesite eruption (Indonesia) and showed a relationship between the timescales (pre-eruptive fluid accumulation a few weeks to months prior to the eruption) and the change of hydroacoustics and water lake temperature but the relationship with seismic swarms was more diffuse<sup>18</sup>.

A comparison between the diffusion timescales found by the previous studies coupling diffusion modelling with monitoring signals with Kizimen ones modelled here is shown in Fig. 5. Compilation of all the data on these studies correlating these two types of information was done to reflect on other studies linking petrological data and monitoring signals (Supplementary Data 5). Kizimen timescales are in good agreement with other pre-eruptive timescales modelled in basaltic andesitic or dacitic volcanoes, such as Ruapehu (1969 eruption), Bezymianny (2012 eruption), which is also in Kamchatka, Llaima (2008 eruption) and the 2014 eruption of Piton de la Fournaise (basaltic) (Fig. 5; Supplementary Data 5). The number of timescales modelled here is high, compared to the other studies and the relationship with the monitoring signals is well correlated (Fig. 5; Supplementary Data 5). These eruptions are thought to have been preceded by magma mixing or magmatic injections before the eruption, which are the scenarios investigated here at Kizimen. These correlations could have implications in terms of monitoring, with unrest periods of more than a year that could lead to eruptions. This long period of unrest has been observed for high silica systems that have experienced long quiescent time where an eruption can occur after a long alert period (from days to years)<sup>85</sup>.

More broadly, a contrast between pre-eruptive timescales obtained on basaltic and silicic volcanoes can then be highlighted, with short timescales obtained in mafic systems compared to more evolved systems<sup>7</sup>. Indeed, magmatic perturbations prior to eruptions emitting basalt or basanite are generally occurring from days to months before the eruption (olivine studies)<sup>13–16,50,51,77–80,86–92</sup>. For more evolved magmas (basaltic andesite/andesitic/dacitic magma; on olivines, orthopyroxenes, clinopyroxenes, plagioclase or apatites)<sup>18,20,21,28,93</sup>, timescales are in the order of a few years<sup>21,81,82</sup> to decades or more<sup>26,94–98</sup>, even if some eruptions have shorter pre-eruptive timescales (~1 year). For dacitic to rhyolitic magmas, using the same crystals as for the previous systems but also quartz, timescales are between one year to longer timescales, such as centuries<sup>17,26,29,49,97,99–106</sup> and also correspond to magma recharge or heating experienced by the magmas. The timescales found for Kizimen 2010 eruption are then in a rather young range compared to these systems.

To conclude, acquiring this type of petrological data on recent and past eruptions on a given volcanic system is of great

importance. Petrological disturbances in a storage system are synchronous with geophysical or geochemical signals followed by monitoring institutions. The strong temporal correspondence between timescales and monitoring data (seismicity, deformation...) that we found at Kizimen (magma mixing that occurred in the reservoir 1.5 years before the eruption, correlated to the onset of a seismic crisis), shows the possible application of this approach on volcanoes whose past eruptions have not been monitored but on which it is possible to estimate timescales by petrology. It is often difficult, particularly on intermediate to silicic volcanic systems, to determine when the first warning signs are recorded and how long it will take before a possible eruption. Moreover, if one can establish a systematics in the timescales for several past eruptions, as it has been established at Montagne Pelée, Martinique<sup>107</sup>, it can be applied more easily for a future eruption. These results indicate that this type of data can be used to estimate the timescales between the first precursory signals and the beginning of an eruption in case of reactivation of a volcano.

## Materials and methods

**Sampling on the field: andesites, dacites, and banded samples.** Samples were collected in the pyroclastic deposits of the 2010–2013 eruption (Fig. 1a) during a three-week field mission which took place in August 2019. Samples coordinates can be found in Supplementary Table 1.

**Sample preparation.** The unaltered dacitic and andesitic samples were crushed (up to 3 mm). One half was then crushed again into fine powder for whole rock analysis while the other half was sieved into different size fractions (710–1000, 500–710, 355–500, 250–355, and 125–250 µm fractions). These fractions were washed in ultrasonic bath and dried at 80 °C for 48 h. They were then observed under the binocular microscope to select the fractions in which the crystals were the more abundant and automorphic. Orthopyroxenes, plagioclase, magnetites, and amphiboles were then handpicked under the binocular microscope in the different size fractions. They were mounted in epoxy resin and polished up to 0.3 µm to the middle part of the crystals. Orthopyroxenes were oriented with c-axis in a north-south direction as they have been used for intracrystalline diffusion modelling<sup>20,26,28,97,108</sup>. Before scanning electron microscopy (SEM) or electron microprobe micro-analysis (EPMA) investigations, selected mounts were all carbon-coated<sup>97</sup>.

**Textural observations: scanning electron microscope.** Orthopyroxenes were observed under scanning electron microscopes (SEM): Zeiss Supra 55VP (Sor-bonne Université, Institut des Sciences de la Terre de Paris (ISTEP), Paris) and the Carl Zeiss EVO MA10 SEM at the PARI platform at the Institut de physique du globe de Paris (Université Paris Cité) using an acceleration voltage of 20 kV and a beam current of 8 nA. An identification of possible chemical zonations in the crystals (grey levels of different intensities) was possible in order to position the rim-core profiles to be carried out with the electron microprobe for orthopyroxenes. Proportions of zoned/unzoned orthopyroxenes crystals were determined on all the crystals mounted by looking closely at the SEM images.

High resolution backscattered electron (BSE) images (mag ×755) were taken for diffusion modelling and intercalibration with chemical profiles of orthopyroxenes with the same parameters mentioned above, a high integration line (N = 7, 7 integrations of each image line to reduce the impact of noise) and a dwell time per pixel of 48 µs.

Melt inclusions in crystals and residual glass areas were also identified using SEM images. Melt inclusions that were of sufficient size for electron microprobe analysis (>6 µm), round in shape, suggesting the absence of crystallization after trapping in the crystals and not containing bubbles (which could indicate outgassing) were selected. The selected glasses areas had a size greater than 10 µm and without microcrystallization.

**Compositional analysis: whole rocks.** Whole rock major-element compositions were analyzed by ICP-OES at Centre de Recherches Pétrographiques et Géochimiques (CRPG) in Nancy<sup>109</sup>.

**Compositional analysis: point measurements.** Orthopyroxenes crystals were analysed for major elements (Si, Al, Ca, Mg, Na, K, Ti, Fe, Mn, and P) by electron microprobe micro-analyzers with an acceleration voltage of 15 kV, a beam current of 10 nA and a focused beam of 2 µm (CAMECA SX-Five and SX-100; Service Camparis, Paris). Counting time on peak and background was set at 10 s for all elements, except for Fe and Mg with 80 s. The core to rim compositional profiles in zoned orthopyroxenes had a 2 µm step and an average length of 100 µm (~4 h per profile). Four points were measured across the unzoned crystals. They were
<figure><img src="image_7.png" /></figure>

**Fig. 5** Database showing the comparison between diffusion timescales for studies correlating diffusion timescales and monitoring signals before or during an eruption. Most of them involve olivines or pyroxenes (orthopyroxenes or clinopyroxenes). The range of timescales that can be modelled with diffusion methods go from hours to centuries⁷ but here, they are mainly between months to a decade before the eruption. Studies with well-correlated diffusion timescales and pre-eruptive seismic data are specified. Pre-eruptive timescales before the 2010 eruption of Kizimen studied here are shown (for andesite (and) and dacite (dac); see arrows). The studies used for comparison are: Ruapehu (1969, 1971, 1977, 1995, and 1996)²⁰; El Hierro (2011–2012)²,⁸⁰; Eyjafjallajökull (2010)¹⁴; Bezymianny (2006–2012)²¹; Etna (1991–1993; 2006)¹³,¹⁵; Piton de la Fournaise (2007; 2014–2015)⁷⁸,⁷⁹; Kelud (1990)¹⁸, Llaima (2008)⁸³; Shishaldin (1999)¹⁶; Mt St Helens (1980–1986)¹⁷; Kīlauea (1959)⁷⁷; Merapi (2006; 2010)⁸¹; Irazú (1963–1965)⁸². Symbols correspond to different minerals and the colors to different magma compositions. The symbol position represents the pre-eruptive timescales given by the minerals (median of the dataset, when all data was available) and error bars correspond to the minimum and maximum timescale modelled¹³–¹⁵,¹⁷,²⁰,⁷⁸,⁷⁹,⁸²,⁸³. For some studies, not all the timescales modelled were provided in the articles (or only the timescales ranges were provided)²,²¹,⁷⁷,⁸⁰,⁸¹ so the values taken were the minimum and maximum timescales given in Figures or in the text of these studies. The comparison between these studies is summarized in Supplementary Data 5.

acquired perpendicular to the long axes of the orthopyroxenes and away from the corners to avoid three-dimensional effects such as growth²⁹,⁹⁷,¹⁰⁸,¹¹⁰.

For compositional analyses in magnetite/ilmenite pairs to estimate the temperatures, the counting time was set at 30 s for Fe, Ti, and Al and 10 s for the other elements.

For the measurement of major and minor elements in the residual glasses and melt inclusions of the samples, an acceleration voltage of 15 kV was used, with counting times varying from 10 to 25 s depending on the element and a beam size of 6–10 µm to avoid alkaline migration¹¹¹–¹¹³.

Analytical errors on major and minor elements were calculated to be of 0.7 wt% for Si, ~0.1 wt% for Na, Al, Ca, K, and Ti, 0.07 wt% for Ni, 0.29 wt% for Mg, 0.39 wt % for Fe, 0.14 wt% for Mn, 0.04 wt% for Cr and P and 0.02 wt% for F and Cl (for melt inclusions and residual glasses).

**Timescales of magmatic processes.** To investigate the changes in crystallization conditions in time and date the magmatic processes, interdiffusion timescales were modelled following the shapes of the concentration's profiles in orthopyroxenes and width of the diffusion zone between the plateaus. Using second Fick's law, if the initial conditions, boundary conditions, and diffusion coefficient of the elements are known, fitting of the profiles can be used to obtain the timescales⁶,⁷,¹³,¹⁵,⁵⁰. The diffusion coefficient depends on parameters such as chemical composition ($X_i$; molar fraction of the mineral constituent element), temperature ($T$ in Kelvin), pressure ($P$ in Pa), oxygen fugacity ($fO_2$) and water

fugacity ($fH_2O$)⁶,⁷. Studies have modelled the timescales on different types of crystals (olivines, plagioclase, clinopyroxenes, orthopyroxenes, quartz...) using different diffusion coefficients and modelling⁷. Some programs and softwares exist to model timescales based on EPMA data, for example in olivines (Diffusion Process Analysis; DIPRA¹¹⁴), or on BSE images in clinopyroxenes (on greyscale data using the Non-Isothermal Diffusion Incremental Step model; NIDIS¹¹⁵) or in orthopyroxenes²⁸.

Modelling of the timescales associated to the zonations in the zoned crystals was done by using the method developed in other studies⁹²,¹⁰⁸, using the parametrization of the interdiffusion coefficient $D$ of Fe and Mg elements in orthopyroxenes¹¹⁶. The Fe-Mg interdiffusion profiles were modelled in one dimension, across the c-axis, parallel to the b-axis of the crystals⁹²,¹⁰⁸. $D$ depends on the composition of the orthopyroxenes, temperature, pressure, oxygen fugacity and water fugacity⁶,²⁶. The interdiffusion coefficient we used here is formulated for a temperature between 500 and 800 °C, an oxygen fugacity of $IW + 0.8$ log units above the IW buffer (IW: Iron–Wüstite buffer) and $X_{Fe} \sim 0.10–0.50$¹¹⁶ ($X_{Fe}$ is the molar fraction of the ferrosilite component). The parametrization of $D$ in our model was done with a range of compositions of orthopyroxenes that includes those of this study (En⁴⁴–⁷³; $X_{Fe} = 0.23–0.53$; Supplementary Fig. 4)¹¹⁶, the temperature of the magmas are slightly above the range used in the parametrization of $D$, which is extrapolated at our temperatures. In the experimental study of $D$ we used to model our timescales, an $fO_2$ dependence similar to olivine with an exponent of 1/6 was hypothesized¹¹⁶. This $fO_2$ dependence has been incorporated in the formulation of $D$ in other works¹⁰⁰,¹⁰⁸ but this dependence has never been
verified experimentally for orthopyroxenes with En contents of our dataset¹¹⁷,¹¹⁸. For En-rich orthopyroxenes of the mantle (En⁹⁸ and En⁹¹), between 870 and 1100 °C and oxygen fugacity controlled (fO₂ = 10⁻¹¹–10⁻⁷ Pa), an fO₂ dependence of an exponent between 0 and 1/20 was found to be appropriate from their data¹¹⁸. The estimate for D<sub>Fe–Mg</sub> used here¹¹⁶, when extrapolated to higher temperatures of the study on En-rich orthopyroxenes without considering any fO₂ dependence appear to be consistent with their dataset¹¹⁸. The oxygen fugacity dependence on Fe–Mg interdiffusion coefficient is still under discussion as incorporating it can give longer timescales⁴⁹ but no study has yet examined the appropriate En compositions of the orthopyroxenes and appropriate oxygen fugacity ranges. Therefore, the parametrization of the interdiffusion coefficient D used is the one determined on orthopyroxenes compositions close from our dataset without an fO₂ correction¹¹⁶, as done in previous studies²⁶,⁹⁷. D is defined as follows in Eq. (1)¹¹⁶:

$$ \log D = -5.54 + 2.6X_{\text{Fe}} - \frac{12530}{T} \quad (1) $$

The main assumption of the model is that the initial profile between two zones follows a step function, which is modified by diffusion to form sigmoidal concentration gradients¹¹⁹ (Supplementary Fig. 8). This steep profile becomes broader with time, due to interdiffusion of the different elements. With the eruption, the diffusion is stopped. Maximum timescales are then estimated¹¹⁹.

Profiles with smooth compositional variations of sigmoid shape were studied and those displaying growth dominant profiles were not considered (Supplementary Fig. 7; Supplementary Data 2). These profiles display linear trends or curves and bumps on the profiles and are not consistent with a diffusion sigmoid calculated for the given initial compositional contrast²⁶ (Supplementary Fig. 7). The Al and Ca gradients were compared to the Mg#, as Al and Ca diffuse slower than Fe–Mg⁴⁹,¹⁰⁸,¹²⁰,¹²¹ (Supplementary Fig. 7; Supplementary Data 2).

Diffusion timescales were modelled using a spreadsheet done by D. J. Morgan (now published as a demo version for olivines⁹².

Examples of modelled profiles can be found in Supplementary Fig. 8 and best fit profiles are shown for all modelled profiles in Supplementary Data 2.

Before using the diffusion model, an intercalibration was done between profiles of greyscale values determined from BSE images of zoned orthopyroxenes using the ImageJ software (https://imagej.nih.gov/ij/; version 1.52a) and with EPMA profiles done at the microprobe in the same crystals. We followed the methodology developed in two other studies⁹²,¹⁰⁸.

The greyscale profiles were done on high-resolution north–south oriented images taken with SEM (spatial resolution of 755 nm), and were averaged over a certain width to reduce the noise (between 20 and 50 µm, depending on the crystal rim size)²⁶,⁹⁷. The Mg number is the parameter that seems to best control the greyscale of the images³³,⁹². Best fit diffusion profiles to the calibrated BSE/EPMA profiles were then modelled using the model spreadsheet²⁶,⁵³,⁹²,¹⁰⁸ to obtain the timescales of interdiffusion. As the diffusion coefficient used is composition-dependent, the shape of the curve of the profile was determined by the degree of compositional contrast but the curve width is dependent on time²⁶,⁹². A database of simulated diffusion profiles obeying composition-dependent diffusion under a 1D diffusion geometry was generated using a finite-difference method²⁶,⁵³,⁹²,¹⁰⁸. As all 1D diffusion profiles are self-similar in time for any specific compositional dependence, the curves are identical in their proportional curvature for initial boundary conditions (such as T, fO₂, anisotropy, and compositional contrasts)⁹². Thus, if the boundary conditions and the initial condition are known for a specific measured diffusion profile, the correct diffusion profile shape can be determined from these simulated diffusion curves as the one with the same degree of compositional contrast and asymmetry characteristics⁹². When the appropriate profile has been identified in the database, a linear stretch is applied along the distance x-axis to fit the data⁹². The spreadsheet measures the distance between the 20th and 80th percentile in the measured profile and then scales the appropriate model curve by the correct n factor before overlaying them (centred on the 50th percentile of both curves as a reference point in the centre of the diffusion range)⁹²,¹⁰⁸.

The timescale corresponding to the zonation is then obtained, as the scaling factor is an expression of crystal composition, the diffusivity used for the library of the profiles, the diffusivity in the crystal and the ratio between the timescales of the database and the timescale associated with the zoning in the crystal²⁶ (with the timescale of the crystal being solved as the only unknown)⁹²,¹⁰⁸.

It is possible to write that t<sub>sample</sub> is:

$$ t_{\text{sample}} = \frac{n^2 D_{\text{model}} t_{\text{model}}}{D_{\text{sample}}} \quad (2) $$

Where n is the stretch factor⁹²,¹⁰⁸. The goodness of the fits is done by fitting the profiles by eye²⁶,⁹⁷.

The determination of the temperature, the diffusion coefficient and its measurement are the main sources of uncertainties⁶,¹¹⁹. Uncertainties on the calibration of measurements and point spacing during EPMA analyses, as well as the resolution of SEM images are also sources of uncertainties (pixel size, uncertainty on greyscale values defining a plateau)⁹⁷. Overall, the measurement of the diffusion coefficient and temperature represent the main sources of uncertainty and have a logarithmic effect on final calculated timescales⁶,²⁶,⁵³ (Fig. 2; Supplementary Fig. 6). Here, a Monte Carlo simulation was used to estimate the integrated time uncertainty from the different sources of uncertainty (on T, D₀, Eₐ and resolution of greyscale values of SEM images), giving asymmetric uncertainties, with a larger error bar on the positive error and a smaller error bar on the negative error²⁶.

**Estimation of earthquake magnitudes.** The Kamchatka Branch of Russian Geophysical Survey uses the following magnitude M scale in Eq. (3):

$$ M = \log_{10} U + f(t_{\text{S-p}}) \quad (3) $$

With f(t<sub>S-p</sub>) expressed in Eq. (4):

$$ f(t_{\text{S-p}}) = \begin{cases} 1.42 \log_{10}(t_{\text{S-p}}) + 0.72, & t_{\text{S-p}} < 7 \text{ s} \\ 2.44 \log_{10}(t_{\text{S-p}}) - 0.14, & t_{\text{S-p}} \ge 7 \text{ s} \end{cases} \quad (4) $$

where U is the maximum S-wave ground velocity in micron/sec and t<sub>S-p</sub> is the arrival time difference between S and P waves. The seismic moments M₀ (in N m) are approximately estimated from magnitudes based on the moment magnitude equation in Eq. (5):

$$ M_0 = 10^{(1.5M+9.05)} \quad (5) $$

## Data availability

Geochemical data and all data used in this study are reported in Supplementary Data 1–5. The seismic catalog used for the analysis is available from the website of the Kamchatka Branch of the Geophysical Survey of Russian Academy of Sciences (http://sdis.emsd.ru/info/earthquakes/catalogue.php). Supplementary Data are archived in Zenodo data repository and are available online (https://doi.org/10.5281/zenodo.7307563).

Received: 25 February 2022; Accepted: 11 November 2022; Published online: 19 November 2022

## References

1. Sparks, R. S. J. Forecasting volcanic eruptions. *Earth Planet. Sci. Lett.* **210**, 1–15 (2003).

2. Marti, J. et al. Correlation of magma evolution and geophysical monitoring during the 2011–2012 El Hierro (Canary Islands) submarine eruption. *J. Petrol.* **54**, 1349–1373 (2013).

3. Jay, J. et al. Locating magma reservoirs using InSAR and petrology before and during the 2011–2012 Cordón Caulle silicic eruption. *Earth Planet. Sci. Lett.* **395**, 254–266 (2014).

4. Caricchi, L., Biggs, J., Annen, C. & Ebmeier, S. The influence of cooling, crystallisation and re-melting on the interpretation of geodetic signals in volcanic systems. *Earth Planet. Sci. Lett.* **388**, 166–174 (2014).

5. McNutt, S. R. *Monitoring and Mitigation of Volcano Hazards* 99–146 (Springer, 1996).

6. Costa, F. & Morgan, D. *Time Constraints from Chemical Equilibration in Magmatic Crystals. Timescales of Magmatic Processes: From Core to Atmosphere* (John Wiley & Sons, Ltd, 2010).

7. Costa, F., Shea, T. & Ubide, T. Diffusion chronometry and the timescales of magmatic processes. *Nat. Rev. Earth Environ.* **1**, 201–214 (2020).

8. Costa, F. Clocks in magmatic rocks. *Annu. Rev. Earth Planet. Sci.* **49**, 231–252 (2021).

9. Costa, F., Dohmen, R. & Chakraborty, S. Time scales of magmatic processes from modeling the zoning patterns of crystals. *Rev. Mineral. Geochem.* **69**, 545–594 (2008).

10. Mollo, S. & Hammer, J. E. *Mineral Reaction Kinetics: Microstructures, Textures, Chemical and Isotopic Signatures* 378–418 (Mineralogical Society of Great Britain & Ireland, 2017).

11. Chakraborty, S. Diffusion in solid silicates: A tool to track timescales of processes comes of age. *Annu. Rev. Earth Planet. Sci.* **36**, 153–190 (2008).

12. Morgan, D. J. et al. Magma chamber recharge at Vesuvius in the century prior to the eruption of A.D. 79. *Geology* **34**, 845 (2006).

13. Kahl, M., Chakraborty, S., Costa, F. & Pompilio, M. Dynamic plumbing system beneath volcanoes revealed by kinetic modeling, and the connection to monitoring data: An example from Mt. Etna. *Earth Planet. Sci. Lett.* **308**, 11–22 (2011).

14. Pankhurst, M. J., Morgan, D. J., Thordarson, T. & Loughlin, S. C. Magmatic crystal records in time, space, and process, causatively linked with volcanic unrest. *Earth Planet. Sci. Lett.* **493**, 231–241 (2018).

15. Kahl, M. et al. Compositionally zoned crystals and real-time degassing data reveal changes in magma transfer dynamics during the 2006 summit eruptive episodes of Mt. Etna. *Bull. Volcanol.* **75**, 692 (2003).

16. Rasmussen, D. J. et al. When does eruption run-up begin? Multidisciplinary insight from the 1999 eruption of Shishaldin volcano. *Earth Planet. Sci. Lett.* **486**, 1–14 (2018).

17. Saunders, K., Blundy, J., Dohmen, R. & Cashman, K. Linking petrology and seismology at an active volcano. *Science* **336**, 1023–1027 (2012).

18. Utami, S. B., Costa, F., Lesage, P., Allard, P. & Humaida, H. Fluid fluxing and accumulation drive decadal and short-lived explosive basaltic andesite eruptions preceded by limited volcanic unrest. *J. Petrol.* **62**, 1–29 (2021).
19. Re, G., Corsaro, R. A., D'Oriano, C. & Pompilio, M. Petrological monitoring of active volcanoes: A review of existing procedures to achieve best practices and operative protocols during eruptions. *J. Volcanol. Geotherm. Res.* **419**, 107365 (2021).

20. Kilgour, G. N. et al. Timescales of magmatic processes at Ruapehu volcano from diffusion chronometry and their comparison to monitoring data. *J. Volcanol. Geotherm. Res.* **288**, 62–75 (2014).

21. Davydova, V. O., Shcherbakov, V. D. & Plechov, P. Y. The timescales of magma mixing in the plumbing system of Bezymianny Volcano (Kamchatka): Insights from diffusion chronometry. *Moscow Univ. Geol. Bull.* **73**, 444–450 (2018).

22. Senyukov, S. L. et al. Seismicity of Kizimen volcano. In *Conf. Third Sci. Tech. Conf.'Problems Integr. Geophys. Monit. Russ. Far East* (2011).

23. Dvigalo, V. N., Melekestsev, I. V., Shevchenko, A. V. & Svirid, I. Y. The 2010–2012 eruption of Kizimen Volcano: The greatest output (from the data of remote-sensing observations) for eruptions in Kamchatka in the early 21st century part I. The November 11, 2010 to December 11, 2011 phase. *J. Volcanol. Seismol.* **7**, 345–361 (2013).

24. Ji, L., Lu, Z., Dzurisin, D. & Senyukov, S. Pre-eruption deformation caused by dike intrusion beneath Kizimen volcano, Kamchatka, Russia, observed by InSAR. *J. Volcanol. Geotherm. Res.* **256**, 87–95 (2013).

25. Firstov, P. P. & Shakirova, A. A. Seismicity observed during the precursory process and the actual eruption of Kizimen volcano, Kamchatka in 2009–2013. *J. Volcanol. Seismol.* **8**, 203–217 (2014).

26. Solaro, C. et al. A system dynamics approach to understanding the deep magma plumbing system beneath dominica (Lesser Antilles). *Front. Earth Sci.* **8**. https://doi.org/10.3389/feart.2020.574032 (2020).

27. Martel, C. et al. Effects of $f_{O_2}$ and $H_2O$ on andesite phase relations between 2 and 4 kbar. *J. Geophys. Res. Solid Earth* **104**, 29453–29470 (1999).

28. Metcalfe, A. et al. Magmatic processes at La Soufrière de Guadeloupe: Insights from crystal studies and diffusion timescales for eruption onset. *Front. Earth Sci.* **9**, 1–28 (2021).

29. Fabbro, G. N., Druitt, T. H. & Costa, F. Storage and eruption of silicic magma across the transition from dominantly effusive to caldera-forming states at an arc volcano (Santorini, Greece). *J. Petrol.* **58**, 2429–2464 (2018).

30. Gorbatov, A., Kostoglovod, V., Suárez, G. & Gordeev, E. Seismicity and structure of the Kamchatka subduction zone. *J. Geophys. Res. Solid Earth* **102**, 17883–17898 (1997).

31. Park, J. et al. A Dangling Slab, Amplified Arc Volcanism, Mantle Flow and Seismic Anisotropy in the Kamchatka Plate Corner. In *Plate Boundary Zones* (eds Stein, S. and Freymueller, J.T.), 295–324 (2002).

32. Ponomareva, V. et al. Late Pleistocene-Holocene volcanism on the Kamchatka Peninsula, Northwest Pacific Region. In *Volcanism and Subduction: The Kamchatka Region* (eds Eichelberger, J., Gordeev, E., Izbekov, P., Kasahara, M. and Lees, J.) 165–198 (2007).

33. Churikova, T. G., Dorendorf, F. & Wörner, G. Sources and fluids in the mantle wedge below Kamchatka, evidence from across-arc geochemical variation. *J. Petrol.* **42**, 1567–1593 (2001).

34. Melekestsev, I. V., Ponomareva, V. V. & Volynets, O. N. Kizimen volcano, Kamchatka—A future Mount St. Helens? *J. Volcanol. Geotherm. Res.* **65**, 205–226 (1995).

35. Shapiro, N. M. et al. Deep and shallow long-period volcanic seismicity linked by fluid-pressure transfer. *Nat. Geosci.* **10**, 442–445 (2017).

36. Levin, V., Shapiro, N., Park, J. & Ritzwoller, M. Seismic evidence for catastrophic slab loss beneath Kamchatka. *Nature* **418**, 763–767 (2002).

37. Koulakov, I. et al. Mantle and crustal sources of magmatic activity of Klyuchevskoy and surrounding volcanoes in Kamchatka inferred from earthquake tomography. *J. Geophys. Res. Solid Earth* **125**, 1–29 (2020).

38. Melnik, O., Lyakhovsky, V., Shapiro, N. M., Galina, N. & Bergal-Kuvikas, O. Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas. *Nat. Commun.* **11**, 3918 (2020).

39. Journeau, C. et al. Seismic tremor reveals active trans-crustal magmatic system beneath Kamchatka volcanoes. *Sci. Adv.* **8** https://doi.org/10.1126/sciadv.abj1571 (2022).

40. Koulakov, I. et al. Three different types of plumbing system beneath the neighboring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka. *J. Geophys. Res. Solid Earth* **122**, 3852–3874 (2017).

41. Auer, A., Belousov, A. & Belousova, M. Deposits, petrology and mechanism of the 2010–2013 eruption of Kizimen volcano in Kamchatka, Russia. *Bull. Volcanol.* **80**, 33 (2018).

42. Churikova, T. G., Wörner, G., Eichelberger, J. & Ivanov, B. Minor- and trace element zoning in plagioclase from Kizimen Volcano, Kamchatka: Insights on the magma chamber processes. In *Volcanism and Subduction: The Kamchatka Region* (eds Eichelberger, J., Gordeev, E., Izbekov, P., Kasahara, M. and Lees, J.) 303–323 (2007).

43. Churikova, T. G. et al. Major and trace element zoning in plagioclase from Kizimen Volcano (Kamchatka): Insights into magma-chamber processes. *J. Volcanol. Seismol.* **7**, 112–130 (2013).

44. Shakirova, A. & Firstov, P. Features of the Kizimen Volcano area seismicity prior to and during the 2010–2013 eruption. *J. Volcanol. Geotherm. Res.* **421**, 107420 (2022).

45. Girona, T., Costa, F., Newhall, C. & Taisne, B. On depressurization of volcanic magma reservoirs by passive degassing. *J. Geophys. Res. Solid Earth* **119**, 8667–8687 (2014).

46. Girona, T., Costa, F. & Schubert, G. Degassing during quiescence as a trigger of magma ascent and volcanic eruptions. *Sci. Rep.* **5**, 18212 (2015).

47. Tembrel, I. & Ovsyannikov, A. Activity of Kizimen volcano at summer 2009. *Bull. Kamchatka Assoc Educ. Center Earth Sci.* **14**, 7–9 (2009).

48. Sauerzapf, U., Lattard, D., Burchard, M. & Engelmann, R. The titanomagnetite-ilmenite equilibrium: New experimental data and thermo-oxymetric application to the crystallization of basic to intermediate rocks. *J. Petrol.* **49**, 1161–1185 (2008).

49. Flaherty, T. et al. Multiple timescale constraints for high-flux magma chamber assembly prior to the Late Bronze Age eruption of Santorini (Greece). *Contrib. Mineral. Petrol.* **173**, 75 (2018).

50. Kahl, M., Chakraborty, S., Pompilio, M. & Costa, F. Constraints on the nature and evolution of the magma plumbing system of Mt. Etna volcano (1991–2008) from a combined thermodynamic and kinetic modelling of the compositional record of minerals. *J. Petrol.* **56**, 2025–2068 (2015).

51. Kahl, M., Viccaro, M., Ubide, T., Morgan, D. J. & Dingwell, D. B. A branched magma feeder system during the 1669 eruption of Mt Etna: Evidence from a time-integrated study of zoned olivine phenocryst populations. *J. Petrol.* **58**, 443–472 (2017).

52. Petrone, C. M. & Mangler, M. F. Elemental Diffusion Chronostratigraphy. In *Crustal Magmatic System Evolution* (eds Masotta, M., Beier, C. and Mollo, S.) 177–193 (2021).

53. Solaro-Müller, C. Storage conditions and dynamics of magma reservoirs feeding the major pumiceous eruptions of Dominica (Lesser Antilles Arc). PhD thesis. Paris (France): Université Paris Diderot, 330 (2017).

54. Waters, L. E. & Lange, R. A. An experimental study of $Fe^{2+}$-MgK between orthopyroxene and rhyolite: a strong dependence on $H_2O$ in the melt. *Contrib. Mineral. Petrol.* **172**, 42 (2017).

55. Morgavi, D., Petrelli, M., Vetere, F. P., González-García, D. & Perugini, D. High-temperature apparatus for chaotic mixing of natural silicate melts. *Rev. Sci. Instrum.* **86**, 105108 (2015).

56. Perugini, D. & Poli, G. The mixing of magmas in plutonic and volcanic environments: Analogies and differences. *Lithos* **153**, 261–277 (2012).

57. Morgavi, D., Arienzo, I., Montagna, C., Perugini, D. & Dingwell, D. B. Magma Mixing: History and Dynamics of an Eruption Trigger. In *Volcanic Unrest* (eds Gottsmann, J., Neuberg, J. and Scheu, B.) 123–137 (2017).

58. Putirka, K. D. Thermometers and barometers for volcanic systems. *Rev. Mineral. Geochem.* **69**, 61–120 (2008).

59. Waters, L. E. & Lange, R. A. An updated calibration of the plagioclase-liquid hygrometer-thermometer applicable to basalts through rhyolites. *Am. Mineral.* **100**, 2172–2184 (2015).

60. Ubide, T. et al. The relevance of crystal transfer to magma mixing: a case study in composite dykes from the central Pyrenees. *J. Petrol.* **55**, 1535–1559 (2014).

61. Sparks, R. S. J. & Marshall, L. A. Thermal and mechanical constraints on mixing between mafic and silicic magmas. *J. Volcanol. Geotherm. Res.* **29**, 99–124 (1986).

62. Frost, T. P. & Mahood, G. A. Field, chemical, and physical constraints on mafic-felsic magma interaction in the Lamarck Granodiorite, Sierra Nevada, California. *Geol. Soc. Am. Bull.* **99**, 272 (1987).

63. Grasset, O. & Albarède, F. Hybridization of mingling magmas with different densities. *Earth Planet. Sci. Lett.* **121**, 327–332 (1994).

64. Bateman, R. The interplay between crystallization, replenishment and hybridization in large felsic magma chambers. *Earth-Sci. Rev.* **39**, 91–106 (1995).

65. Perugini, D., Poli, G., Christofides, G. & Eleftheriadis, G. Magma mixing in the Sithonia Plutonic Complex, Greece: evidence from mafic microgranular enclaves. *Mineral. Petrol.* **78**, 173–200 (2003).

66. Perugini, D., De Campos, C. P., Dingwell, D. B., Petrelli, M. & Poli, G. Trace element mobility during magma mixing: preliminary experimental results. *Chem. Geol.* **256**, 146–157 (2008).

67. Perugini, D., Poli, G., Petrelli, M., De Campos, C. P. & Dingwell, D. B. Time-scales of recent Phlegrean Fields eruptions inferred from the application of a ‘diffusive fractionation’ model of trace elements. *Bull. Volcanol.* **72**, 431–447 (2010).

68. Perugini, D., De Campos, C. P., Petrelli, M. & Dingwell, D. B. Concentration variance decay during magma mixing: a volcanic chronometer. *Sci. Rep.* **5**, 14225 (2015).

69. Montagna, C. P., Papale, P. & Longo, A. Timescales of mingling in shallow magmatic reservoirs. *Geol. Soc. London, Spec. Publ.* **422**, 131–140 (2015).

70. Caricchi, L., Sheldrake, T. E. & Blundy, J. Modulation of magmatic processes by $CO_2$ flushing. *Earth Planet. Sci. Lett.* **491**, 160–171 (2018).

71. Ridolfi, F., Renzulli, A. & Puerini, M. Stability and chemical equilibrium of amphibole in calc-alkaline magmas: an overview, new thermobarometric
formulations and application to subduction-related volcanoes. *Contrib. Mineral. Petrol.* **160**, 45–66 (2010).

72. Ridolfi, F. & Renzulli, A. Calcic amphiboles in calc-alkaline and alkaline magmas: thermobarometric and chemometric empirical equations valid up to 1130 °C and 2.2 GPa. *Contrib. Mineral. Petrol.* **163**, 877–895 (2012).

73. d'Augustin, T. Les éléments halogènes dans les magmas, du traçage des conditions de stockage aux flux éruptifs. PhD thesis. Paris (France): Sorbonne Université, 315 (2021).

74. Shinohara, H. Excess degassing from volcanoes and its role on eruptive and intrusive activity. *Rev. Geophys.* **46**, RG4005 (2008).

75. Seropian, G., Kennedy, B. M., Walter, T. R., Ichihara, M. & Jolly, A. D. A review framework of how earthquakes trigger volcanic eruptions. *Nat. Commun.* **12**, 1004 (2019).

76. Cassidy, M. et al. Explosive eruptions with little warning: Experimental petrology and volcano monitoring observations from the 2014 eruption of Kelud, Indonesia. *Geochem., Geophys. Geosyst.* **20**, 4218–4247 (2019).

77. Rae, A. S. P. et al. Time scales of magma transport and mixing at Kilauea Volcano, Hawaii. *Geology* **44**, 463–466 (2016).

78. Albert, H. et al. Magma interactions, crystal mush formation, timescales, and unrest during caldera collapse and lateral eruption at ocean island basaltic volcanoes (Piton de la Fournaise, La Réunion. *Earth Planet. Sci. Lett.* **515**, 187–199 (2019).

79. Sundermeyer, C., Di Muro, A., Gordeychik, B. & Wörner, G. Timescales of magmatic processes during the eruptive cycle 2014–2015 at Piton de la Fournaise, La Réunion, obtained from Mg–Fe diffusion modelling in olivine. *Contrib. Mineral. Petrol.* **175**, 1 (2020).

80. Longpre, M.-A., Klugel, A., Diehl, A. & Stix, J. Mixing in mantle magma reservoirs prior to and during the 2011–2012 eruption at El Hierro, Canary Islands. *Geol. Soc. Am.* **42**, 315–318 (2014).

81. Costa, F., Andreastuti, S., Bouvet de Maisonneuve, C. & Pallister, J. S. Petrological insights into the storage conditions, and magmatic processes that yielded the centennial 2010 Merapi explosive eruption. *J. Volcanol. Geotherm. Res.* **261**, 209–235 (2013).

82. Ruprecht, P. & Plank, T. Feeding andesitic eruptions with a high-speed connection from the mantle. *Nature* **500**, 68–72 (2013).

83. Ruth, D. C. S. et al. Crystal and melt inclusion timescales reveal the evolution of magma migration before eruption. *Nat. Commun.* **9**, 2657 (2018).

84. Albert, H., Costa, F. & Martí, J. Years to weeks of seismic unrest and magmatic intrusions precede monogenetic eruptions. *Geology* **44**, 211–214 (2016).

85. Passarelli, L. & Brodsky, E. E. The correlation between run-up and repose times of volcanic eruptions. *Geophys. J. Int.* **188**, 1025–1045 (2012).

86. Hartley, M. E., Morgan, D. J., MacLennan, J., Edmonds, M. & Thordarson, T. Tracking timescales of short-term precursors to large basaltic fissure eruptions through Fe–Mg diffusion in olivine. *Earth Planet. Sci. Lett.* **439**, 58–70 (2016).

87. Brenna, M. et al. Olivine xenocryst diffusion reveals rapid monogenetic basaltic magma ascent following complex storage at Pupuke Maar, Auckland Volcanic Field, New Zealand. *Earth Planet. Sci. Lett.* **499**, 13–22 (2018).

88. Morgado, E. et al. Transient shallow reservoirs beneath small eruptive centres: constraints from Mg-Fe interdiffusion in olivine. *J. Volcanol. Geotherm. Res.* **347**, 327–336 (2017).

89. Lynn, K. J., Garcia, M. O., Shea, T., Costa, F. & Swanson, D. A. Timescales of mixing and storage for Keanakäko'i Tephra magmas (1500–1820 C.E.), Kilauea Volcano, Hawaii. *Contrib. Mineral. Petrol.* **172**, 76 (2017).

90. Petrone, C. M., Braschi, E., Francalanci, L., Casalini, M. & Tommasini, S. Rapid mixing and short storage timescale in the magma dynamics of a steady-state volcano. *Earth Planet. Sci. Lett.* **492**, 206–221 (2018).

91. Caracciolo, A. et al. Timescales of crystal mush mobilization in the Bárðarbunga-Veiðivötn volcanic system based on olivine diffusion chronometry. *Am. Mineral.* **106**, 1083–1096 (2021).

92. Couperthwaite, F. K., Thordarson, T., Morgan, D. J., Harvey, J. & Wilson, M. Diffusion timescales of magmatic processes in the Moinui lava eruption at Mauna Loa, Hawaii, as inferred from bimodal olivine populations. *J. Petrol.* **61**, 1–19 (2020).

93. Humphreys, M. C. S. et al. Rapid pre-eruptive mush reorganisation and atmospheric volatile emissions from the 12.9 ka Laacher See eruption, determined using apatite. *Earth Planet. Sci. Lett.* **576**, 117198 (2021).

94. Costa, F. & Chakraborty, S. Decadal time gaps between magmatic intrusion and silicic eruption obtained from chemical zoning patterns in olivine. *Earth Planet. Sci. Lett.* **227**, 517–530 (2004).

95. Sundermeyer, C., Gärtjen, J., Weimann, L. & Wörner, G. Timescales from magma mixing to eruption in alkaline volcanism in the Eifel volcanic fields, western Germany. *Contrib. Mineral. Petrol.* **175**, 77 (2020).

96. Rout, S. S. & Wörner, G. Zoning and exsolution in alkali feldspars from Laacher See volcano (Western Germany): constraints on temperature history prior to eruption. *Contrib. Mineral. Petrol.* **173**, 95 (2018).

97. Ostorero, L. et al. Time - window into the transcrustal plumbing system dynamics of Dominica (Lesser Antilles). *Sci. Rep.* **11**, 1–15 (2021).

98. Weber, G., Arce, J. L., Ulianov, A. & Caricchi, L. A recurrent magmatic pattern on observable timescales prior to plinian eruptions from Nevado de Toluca (Mexico). *J. Geophys. Res. Solid Earth* **124**, 10999–11021 (2019).

99. Barker, S. J., Wilson, C. J. N., Morgan, D. J. & Rowland, J. V. Rapid priming, accumulation, and recharge of magma driving recent eruptions at a hyperactive caldera volcano. *Geology* **44**, 323–326 (2016).

100. Allan, A. S. R. et al. A cascade of magmatic events during the assembly and eruption of a super-sized magma body. *Contrib. Mineral. Petrol.* **172**, 49 (2017).

101. Druitt, T. H., Costa, F., Deloule, E., Dungan, M. & Scaillet, B. Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. *Nature* **482**, 77–80 (2012).

102. Chamberlain, K. J., Morgan, D. J. & Wilson, C. J. N. Timescales of mixing and mobilisation in the Bishop Tuff magma body: perspectives from diffusion chronometry. *Contrib. to Mineral. Petrol.* **168**, 1034 (2014).

103. Singer, B. S., Costa, F., Herrin, J. S., Hildreth, W. & Fierstein, J. The timing of compositionally-zoned magma reservoirs and mafic 'priming' weeks before the 1912 Novarupta-Katmai rhyolite eruption. *Earth Planet. Sci. Lett.* **451**, 125–137 (2016).

104. Andersen, N. L. et al. Petrochronologic perspective on rhyolite volcano unrest at Laguna del Maule, Chile. *Earth Planet. Sci. Lett.* **493**, 57–70 (2018).

105. Till, C. B., Vazquez, J. A. & Boyce, J. W. Months between rejuvenation and volcanic eruption at Yellowstone caldera, Wyoming. *Geology* **43**, 695–698 (2015).

106. Matthews, N. E., Huber, C., Pyle, D. M. & Smith, V. C. Timescales of Magma Recharge and Reactivation of Large Silicic Systems from Ti Diffusion in Quartz. *J. Petrol.* **53**, 1385–1416 (2012).

107. Boudon, G., Balcone-Boissard, H. & Morgan, D. J. Systematic pre-eruptive dynamic of the magma plumbing system leading to Plinian eruption at Montagne Pelée Martinique (Lesser Antilles). In *Abstract Cities on Volcanoes* **10** (2018).

108. Allan, A. S. R., Morgan, D. J., Wilson, C. J. N. & Millet, M.-A. From mush to eruption in centuries: assembly of the super-sized Oruanui magma body. *Contrib. to Mineral. Petrol.* **166**, 143–164 (2013).

109. Carignan, J., Hild, P., Mevelle, G., Morel, J. & Yeghicheyan, D. Routine Analyses of Trace Elements in Geological Samples using Flow Injection and Low Pressure On-Line Liquid Chromatography Coupled to ICP-MS: A Study of Geochemical Reference Materials BR, DR-N, UB-N, AN-G and GH. *Geostand. Geoanalytical Res.* **25**, 187–198 (2001).

110. Krimer, D. & Costa, F. Evaluation of the effects of 3D diffusion, crystal geometry, and initial conditions on retrieved time-scales from Fe–Mg zoning in natural oriented orthopyroxene crystals. *Geochim. Cosmochim. Acta* **196**, 271–288 (2017).

111. Balcone-Boissard, H. et al. Deep pre-eruptive storage of silicic magmas feeding Plinian and dome-forming eruptions of central and northern Dominica (Lesser Antilles) inferred from volatile contents of melt inclusions. *Contrib. to Mineral. Petrol.* **173**, 101 (2018).

112. Solaro, C. et al. Petrological and experimental constraints on magma storage for large pumiceous eruptions in Dominica island (Lesser Antilles). **81**, 55 (2019).

113. d'Augustin, T. et al. Evidence for an active, transcrustal magma system in the last 60 ka and eruptive degassing budget ($H_2O, CO_2$, S, F, Cl, Br): The case of Dominica. *Geochem., Geophys. Geosyst.* **21**, https://doi.org/10.1029/2020GC009050 (2020).

114. Girona, T. & Costa, F. DIPRA: A user-friendly program to model multi-element diffusion in olivine with applications to timescales of magmatic processes. *Geochem., Geophys. Geosyst.* **14**, 422–431 (2013).

115. Petrone, C. M., Bugatti, G., Braschi, E. & Tommasini, S. Pre-eruptive magmatic processes re-timed using a non-isothermal approach to magma chamber dynamics. *Nat. Commun.* **7**, 12946 (2016).

116. Ganguly, J. & Tazzoli, V. Fe²⁺-Mg interdiffusion in orthopyroxene: retrieval from the data on intracrystalline exchange reaction. *Am. Mineral.* **79**, 930–937 (1994).

117. Stimpfl, M., Ganguly, J. & Molin, G. Kinetics of Fe²⁺-Mg order-disorder in orthopyroxene: experimental studies and applications to cooling rates of rocks. *Contrib. Mineral. Petrol.* **150**, 319–334 (2005).

118. Dohmen, R., Ter heege, J. H., Becker, H.-W. & Chakraborty, S. Fe-Mg interdiffusion in orthopyroxene. *Am. Mineral.* **101**, 2210–2221 (2016).

119. Morgan, D. J. et al. Time scales of crystal residence and magma chamber volume from modelling of diffusion profiles in phenocrysts: Vesuvius 1944. *Earth Planet. Sci. Lett.* **222**, 933–946 (2004).

120. Smith, D. & Barron, B. R. Pyroxene-garnet equilibration during cooling in the mantle. *Am. Mineral.* **76**, 1950–1963 (1991).

121. Cherniak, D. J. & Dimanov, A. Diffusion in Pyroxene, Mica, and Amphibole. *Rev. Mineral. Geochem.* **72**, 641–690 (2010).

122. Minster, J. B. & Jordan, T. H. Present-day plate motions. *J. Geophys. Res.* **83**, 5331 (1978).
123. Konstantinovskaia, E. Arc-continent collision and subduction reversal in the Cenozoic evolution of the Northwest Pacific: an example from Kamchatka (NE Russia). *Tectonophysics* **333**, 75–94 (2001).

124. Porter, C. et al. "ArcticDEM" https://doi.org/10.7910/DVN/OHHUKH (2018).

125. Portnyagin, M. & Manea, V. C. Mantle temperature control on composition of arc magmas along the Central Kamchatka Depression. *Geology* **36**, 519 (2008).

126. Viccaro, M., Giuffrida, M., Nicotra, E. & Ozerov, A. Y. Magma storage, ascent, and recharge history prior to the 1991 eruption at Avachinsky Volcano, Kamchatka, Russia: inferences on the plumbing system geometry. *Lithos* **140–141**, 11–24 (2012).

## Acknowledgements

We thank the editors E. Liu and J. Aslin and reviewers T. Ubide as well as an anonymous reviewer for their constructive comments and suggestions that helped to improve the manuscript. We would like to thank S. Hidalgo for her help with sample preparation, M. Fialin and N. Rividi for assistance during EPMA analyses, as well as O. Boudouma and S. Borensztajn for SEM imaging and E. Delairis for the thin sections. We also thank D. J. Morgan for providing the files to model the timescales of the orthopyroxenes. We also acknowledge that the image from Fig. 1a is from a Digital Elevation Model provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736 and Fig. 1b is from NASA Earth Observatory images, courtesy of the U.S. Geological Survey. Thank you also to R. Bajou for the discussion on the cumulative number of earthquakes and to F. Costa and S. B. Utami for the discussions on the Fe²⁺-Mg exchange coefficient K_D. This project was funded by the Université Paris Cité, Institut de physique du globe de Paris (doctoral grant from the French Ministry of Higher Education and Research and Innovation) and ANR V-Care-18-CE03-0010. The fieldwork was conducted with permission of the Kronotsky Natural Reserve, obtained by A. Belousov and M. Belousova. The work of N. Shapiro, S. Senyukov, and S. Droznina has been supported by the Ministry of Science and Education of the Russian Federation under project no.

14.W03.31.0033 "Geophysical research, monitoring and forecasting of catastrophic geodynamic processes in the Russian Far East". The work of N. Shapiro has also been supported by the European Research Council under the European Union Horizon 2020 research and innovation program (Grant Agreement 787399-SEISMAZE). An open access license has been selected.

## Author contributions

L.O., H.B.B., G.B., A.B., M.B., and A.A. participated in the field mission in 2019 to collect the samples in Kamchatka (Russia). L.O. performed all data acquisition, prepared the figures, and wrote the initial version of the manuscript during her PhD. L.O., H.B.B., and

G.B. participated in the interpretation of the data and the writing of the manuscript text. N.S., S.S., and S.D. analyzed the seismic catalog of Kizimen and participated in the interpretation of the seismic data. All authors contributed to the writing of the manuscript text.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00622-3.

**Correspondence and requests for materials should be addressed to Lea Ostorero.**

**Peer review information** *Communications Earth & Environment* thanks Teresa Ubide and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Emma Liu, Joe Aslin. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022