Tracing Volatiles, Halogens, and Chalcophile Metals during Melt Evolution at the Tolbachik Monogenetic Field, Kamchatka

Melt storage and supply beneath arc volcanoes may be distributed between a central stratovolcano and wider fields of monogenetic cones, indicating complex shallow plumbing systems. However, the impact of such spatially variable magma storage conditions on volatile degassing and trace element geochemistry is unclear. This study explores magma generation and storage processes beneath the Tolbachik volcanic field, Kamchatka, Russia, in order to investigate the evolution of the magmatic volatile phase and, specifically, the strong enrichment of chalcophile metals (in particular, Cu) in this system. We present new geochemical data for a large suite of olivine-and clinopyroxene-hosted melt inclusions (and host phenocrysts) from five separate monogenetic cones within the Tolbachik volcanic field. These high-Al composition magmas likely reflect the homogenised fractionation products of primitive intermediate-Mg melt compositions, stored at shallow depths after significant fractional crystallisation. Boron isotope compositions and incompatible trace element ratios of the melt inclusions suggest a deeper plumbing system that is dominated by extensive fractional crystallisation and fed by melts derived from an isotopically homogeneous parental magma composition. Volatile components (H 2 O, CO 2 , S, Cl, F) show that magmas feeding different monogenetic cones had variable initial volatile contents and subsequently experienced different fluid-saturated storage conditions and degassing histories. We also show that melts supplying the Tolbachik volcanic field are strongly enriched in Cu compared with almost all other Kamchatka rocks, including samples from the Tolbachik central stratocones, and other volcanoes situated in close proximity in the Central Kamchatka Depression. The melt inclusions record Cu concentrations ≥ 450 μ g/g at ca . 4–5 wt.% MgO, which can only be explained by bulk incompatible partitioning behaviour of Cu, i.e. evolution under sulphide-undersaturated conditions. We suggest that initial mantle melting in this region exhausted mantle sulphides, leading to sulphide undersaturated primitive melts. This sulphide-free model for the high-Al cone melts is further supported by S/Se and Cu/Ag values that overlap those of the primitive mantle and MORB array, with bulk rock Cu/Ag ratios also overlapping other with other global arc datasets for magma evolution prior to fractionation of a monosulfide solid solution. We therefore demonstrate that the combination of novel chalcophile metal analyses with trace element, isotopic, and volatile data is a powerful tool for deciphering complex magmatic evolution conditions across the entire volcanic field.


INTRODUCTION
The processes of magma generation, migration, and closed-vs open-system storage and evolution beneath large monogenetic cone fields are being continually refined (e.g.Smith & Németh, 2017).Geochemical trends recorded by some monogenetic eruptions within the same volcanic field are commonly interpreted to result from the evolution of a single injection of parental magma.Conversely, other more complex compositional variations have been shown to result from the mixing of temporally related but chemically discrete magma sources (e.g.Brenna et al., 2011).The complexity of monogenetic basaltic eruptions is largely governed by the geometry and depth of the plumbing system feeding the volcanic cones, and characteristically rapid magma ascent rates help to prevent the assimilation of country rock, thus preserving original geochemical signatures (Smith & Németh, 2017, and the references therein).However, investigation of whole rocks alone can be problematic and obscure key geochemical trends.Instead, targeted analysis of melt inclusions offers a direct record of evolving magma compositions prior to eruption and shallow degassing and the associated changes in volatile and metal contents (e.g.Rose-Koga et al., 2021, and the references therein).
Here, we collect new data for volatile components, including halogens (H 2 O, CO 2 , S, Cl, F), B isotopes, incompatible trace elements, and chalcophile metals (Se, As, Cu, Ag) for a suite of >150 melt inclusions (and corresponding bulk rocks) originating from several Holocene monogenetic scoria (± lava) cones across the Tolbachik monogenetic field, Kamchatka.The high density of large, glassy, naturally quenched melt inclusions in these samples render them an ideal candidate to decipher the complete history of how volatiles, halogens, and chalcophile metals evolve in magmas erupted from spatially proximal basaltic monogenetic cones, in a complex magma storage region (Portnyagin et al., 2015;Churikova et al., 2015b;Kugaenko & Volynets, 2019).
We interpret the new data within the existing framework of major element, trace element, and isotopic data for different eruptive products of the Tolbachik monogenetic field (e.g.Supplementary material of Churikova et al., 2015a).Our new data reveal how melts from these monogenetic cones differ from those erupted at the older stratocone edifice and other more primitive cones (e.g.Portnyagin et al., 2015;Churikova et al., 2015aChurikova et al., , 2015b)).We also show that our new δ 11 B data fingerprint an isotopically homogenised melt reservoir that was likely a common source to all the monogenetic cones, prior to smaller batches derived from this common source undergoing their own fractional crystallisation ± mixing histories before individual eruptions.
Additionally, melts from the Tolbachik volcanic massif are anomalously Cu-rich (Vergasova & Filatov, 2012;Portnyagin et al., 2015;Mironov & Portnyagin, 2018), with magma evolution driving melts to higher Cu contents despite evidence for anomalous early and extensive silicate-sulphide immiscibility (Kamenetsky et al., 2017;Zelenski et al., 2018Zelenski et al., , 2022)).This is in contrast to the decrease in Cu recorded in other sulphide-saturated systems globally (Jenner et al., 2010;Chiaradia, 2014).We use the chalcophile metal systematics of these melt inclusions to illustrate the contrasting behaviour of Cu, Se, and Ag during the evolution of these arc magmas.The Cu-rich nature of the Tolbachik system requires evolution under sulphide-free conditions and illustrates how key trends in chalcophile metal ratios (S/Se, Cu/Se, Cu/Ag) can be used to further trace the mobility of S during melt degassing, the recycling of early formed immiscible magmatic sulphides and their metal budgets (e.g.Wieser et al., 2020) and may also reveal information about the processes of mantle melting (e.g. de Hoog et al., 2001;Zhang et al., 2021;Deng et al., 2022).

GEOLOGIC BACKGROUND AND SAMPLE LOCALITIES
The late Pleistocene-Holocene Tolbachik volcanic massif is located in the south of the Klyuchevskaya volcanic group, situated in the north of the Central Kamchatka Depression, a zone interpreted to be undergoing intra-arc extension (Fig. 1A; Churikova et al., 2015b).The massif comprises two stratovolcanoes, Ostry Tolbachik and Plosky Tolbachik, as well as scoria cones and related lava f lows, and is one of the largest and most active volcanic regions in Kamchatka.Since the beginning of the Holocene eruptive activity has been focussed along a SW-NE fissure system that trends through the summit of the Plosky Tolbachik volcano (Fig. 1B).The length of the entire fissure zone is ca.60 km and it is ca.13-15 km wide.Numerous monogenetic cones, including those of recent historical eruptions (1941, 1975-1976, and 2012-2013), are all related to this fissure system.An in-depth review of the eruption history of the massif, including petrology and physical volcanology, prior to the 2012-2013 eruption is available in Churikova et al. (2015a, and the references therein).Further details about the most recent 2012-2013 fissure eruption are presented by Belousov et al. (2015) and Volynets et al. (2015).Phenocrysts and melt inclusions from some of the historical eruptions have also been studied in greater detail (e.g.1941 eruption (Kamenetsky et al., 2017;Gurenko, 2021;Zelenski et al., 2022Zelenski et al., ), 1975Zelenski et al., -1976 eruption (Mironov & Portnyagin, 2018), 2012-2013 eruption (Plechov et al., 2015)), alongside bulk rock measurements (Portnyagin et al., 2015), and therefore provide comparative context to the Holocene monogenetic cones of this study, as discussed in later sections.
Broadly, the composition of erupted melts from the Tolbachik volcanic massif varies between high-Mg basalt and high-Al basalt, the latter comprising ∼90% of all erupted material (Braitseva et al., 1984), and with a suite of intermediate hybrids.General liquidus phenocryst assemblages of the Tolbachik rocks are dominated by olivine ± clinopyroxene ± plagioclase ± Cr-spinel (Churikova et al., 2015a).The whole suite can also be separated into distinct middle-K and high-K rock series that straddle the tholeiitic-calc-alkaline boundary and that show major and trace element differences (e.g.Churikova et al., 2015b).Some modelling has suggested that the compositional variation produced by the Tolbachik plumbing system requires long-term open-system fractionation, accompanied by periodic replenishment and mixing between variably evolved melt batches derived from a primary, middle-K mantlederived picrite melt (Portnyagin et al., 2015).Other modelling by Churikova et al. (2015b) suggests that regional intra-arc rifting has allowed the establishment of the Holocene fissure system, in turn facilitating the upwelling and degassing of melts at shallow depths.This geodynamic change, relative to the periods of older stratovolcano magmatism, has shifted the conditions from waterrich (middle-K rocks) to relatively anhydrous (high-K rocks) and where extensive fractional crystallisation exerts the primary control on melt evolution.Following the convention of Churikova et al. (2015b), throughout the manuscript we refer to 'Trend 1' rocks as those related to the central stratovolcanoes and associated feeder dykes; whereas 'Trend 2' rocks are those erupted from the active monogenetic cone and fissure system.
The samples in this study are taken from five separate high-Al monogenetic cones (i.e.'Trend 2') erupted along the SW-NE fissure system (Fig. 1B).Four samples are located along the sparsely studied NE-limb of the fissure, extending to the northern end of the system between Zimina and Bezymianny volcanoes, whereas the other sample is taken from the Peschanye Gorki cone along the SW-limb.The exact latitudes and longitudes of each cone are available in Supplementary Table 1.Cone ZIM-16-07 is coarsely age constrained to the Holocene, whereas cones TOL-16-08, TOL-16-10A, and TOL-16-14 have been dated at 1.4-7.8ka, and TOL-16-43 was erupted 0.8-1.4ka.The reader is referred to Flerov et al. (1984) and Churikova et al. (2015a, and the references therein) for a comprehensive map of the Tolbachik volcanic field that shows both ages and compositions of individual vents and lava f lows, including the ones in this study.This study also reports previously unpublished trace element concentrations of melt inclusions from samples of the 2012-2013 eruption (Naboko cone; Fig. 1B), reported by Plechov et al. (2015), with these data available in the Supplementary Data A1.
While some cones investigated in this study produce eruptive products of both scoria and lava, only scoriaceous material was sampled to target naturally glassy melt inclusions.The scoria samples are dominated by olivine (± clinopyroxene) phenocrysts, with variably glassy groundmasses that contain plagioclase + olivine ± clinopyroxene ± magnetite/Ti-magnetite (± minor Crspinel) microphenocrysts (up to ca.50 μm length).Olivine phenocrysts show little evidence for resorbed crystal faces or disequilibrium morphologies.Olivine crystals (and clinopyroxene; from sample TOL-16-43 only) were separated from the bulk scoria by hand crushing and picking.Melt inclusion-bearing phenocrysts were extracted, mounted in epoxy resin, and inclusions exposed by non-selective polishing.Most phenocrysts contained multiple, large (>50 μm) naturally glassy melt inclusions, readily visible in transmitted light, and equally distributed throughout core to rim zones (Fig. 2B).Only a minority of the inclusions contained small shrinkage bubbles and/or spinel inclusions (Fig. 2C), and inclusions were visually free of immiscible sulphide droplets.Remnants of groundmass glasses adhered to the phenocryst rims were also free of sulphides.The inclusions were imaged with backscattered electrons prior to any analysis and subsequently reimaged under ref lected light following each analytical technique to confirm the exact location of the analysis spot, as outlined in the next sections.The reader is referred to the accompanying ZIP archive file for ref lected and transmitted light images of all melt inclusions in this study (pre-and post-analysis), along with representative groundmass BSE images from each scoria sample.

Electron microprobe
The major, minor, and volatile (S, Cl) compositions of the melt inclusion glasses and host phenocryst phases were determined using a JEOL JXA-8230 electron microprobe (EPMA) at the University of Leeds Electron Microscopy and Spectroscopy Centre.An accelerating voltage of 15 kV, a beam current of 30 nA, and a focussed beam of 2 μm diameter was used for crystal analyses.A two-step routine was employed for the melt inclusion glasses, where major and minor element concentrations were determined first using a defocussed (10 μm) beam at 10 nA, followed by determination of Cl and S concentrations in the same analysis spot location using a 20 nA beam current with extended counting times and aggregated spectrometer intensities to improve analytical precision and detection limits.Full analytical details, data calibration, and standard reproducibility information are available in the Supplementary Methods of Iveson et al. (2021).

Secondary ion mass spectrometry (SIMS)
Following EPMA analysis, the melt inclusions were analysed by secondary ion mass spectrometry (SIMS) across two sessions at the NERC Edinburgh Ion Microprobe Facility for H 2 O, CO 2 , F, Cl, and trace elements (Cameca IMS-4f) and subsequently for 11 B/ 10 B (Cameca IMS-1270).Full details of the analytical procedures and calibration standard reproducibility data are available in the Supplementary Methods of Iveson et al. (2021).As a confirmation of analytical consistency, a comparison of Cl concentrations in the melt inclusions determined by both EPMA and SIMS methods gives excellent agreement between the two techniques (Supplementary Fig. 1A).

Melt inclusion and host crystal trace elements
Trace element analyses of the melt inclusions and selected host phenocrysts were conducted at The School of Environment, Earth and Ecosystem Sciences, The Open University, using a Photons Machines Analyte Excimer 193 nm laser coupled to an Agilent 8800 triple quadrupole ICP-MS (ICP-QQQ) and following the analytical protocols established in Jenner & O'Neill (2012).The LA-ICP-MS analyses were conducted on the same mounts prepared for, and analysed by, EPMA and SIMS (Fig. 2) that were re-polished with 0.3 μm Al 2 O 3 slurry following the SIMS B-isotope analyses to remove any remaining Au coating.
A laser rep rate of 10 Hz, a f luence of 3.54 J/cm 2 and 30-s background collection followed by 30 s of signal were used for all samples analysed.Host phenocryst data were acquired with a 110-μm beam diameter, whereas melt inclusions were ablated with either a 25-or 50-μm beam depending on the size and shape of the inclusions (e.g.Fig. 2D).Laser aerosols were transported away from the ablation site by He (0.9 L/min) and downstream mixed with 0.77 L/min of Ar in a mixing bulb before introduction into the plasma.A total of 5 mL/min of N 2 was added to the carrier gas for analyses conducted using a 25-μm spot size to increase the ionisation efficiency (Li data were not collected for analyses using a 25-μm spot because addition of N 2 compromises the accuracy of Li analyses).All melt inclusions were also subject to four preablation laser pulses (with the same laser energy and spot size as the sample spots) to remove any residual surface material prior to collection of the spectra.
External calibration of data was conducted using NIST SRM-612.Reference materials BCR-2G and KL2-G were periodically analysed to determine the accuracy and precision of the data (full reproducibility data are given in the Supplementary Data A2).EPMA SiO 2 abundances for both melt inclusions and host crystals were used for internal calibration of LA-ICP-MS data.The collected spectra were processed using the Iolite software package (v. 3.71;Paton et al., 2011) and manually screened to remove data that clearly showed any signal indicative of an unintentional analysis of microphenocrysts or melt inclusions (e.g.strongly increased Mg or Ni intensity in melt inclusion spectra or resolvable Rb or Cs signal in host crystal spectra).Consistent with the visual absence of sulphide droplets, all Cu analyses showed a f lat signal intensity during the ablation period, indicative of a homogeneous Cu distribution in the melt inclusion glass itself.As with EPMA vs SIMS, a comparison of various trace element concentrations measured by both SIMS and LA-ICP-MS methods also gives excellent analytical agreement (Supplementary Figs 1B-D).

Melt inclusion Se and As analysis
( 80 Se + → 96 SeO + ; 75 As + → 91 AsO + ), as per the methods outlined in Wieser et al. (2020).Melt inclusions were ablated in linear transects using a 50-μm laser spot, 50 Hz rep.rate, 15-μm/sec scan speed, and variable transect lengths, beginning and ending in the host olivine.The resultant signal intensities were screened to include only the portion of the analysis transect wholly within the melt inclusion glass (Supplementary Fig. 2).NIST SRM-612 was used as for external calibration of data, and accuracy and precision were monitored using in-house Se and As standards, with repeat analysis yielding values within 1%-7% of the preferred values described in Jenner et al. (2015) and Wieser et al. (2020) (also see Supplementary Data A2).

Whole-rock x-ray fluorescence, solution ICP-MS, and multi-collector-ICP-MS
The major and minor oxide concentrations of the bulk scoria were measured by Bureau Veritas Mineral Laboratories, at the Poland-Krakow facility.Fresh, unaltered sample was mechanically crushed, split, and pulverised and then fused into glass discs with lithium tetraborate/metaborate f lux.The glass discs were then analysed by standard x-ray f luorescence methods, and standard reproducibility data are available in Supplementary Data A3.
For solution ICP-MS, a 100-mg mass of fresh sample was subject to a multi-stage HF-HNO 3 digestion process, with low dry down temperatures (65 • C) between each dissolution step to minimise volatile element loss.This dissolution process was undertaken in the isotope clean laboratory of the Durham University Earth Science Geochemistry Centre.Solutions for analysis were made up to 1000-fold dilutions of the original powder weight in 2% HNO 3 , including method and acid blanks.The trace element concentrations of the final solutions were then analysed using the Agilent 8800 ICP-QQQ at The Open University.Further details of the analytical technique are available in Cox et al. (2019), including the specific operating conditions for Se and Ag analyses.An unknown sample and the BHVO-2 reference material were measured every five or six unknown samples as a monitor block, to measure and correct for instrument drift.Accuracy and reproducibility data for BVHO-2 standard are available in Supplementary Data A2.
The B isotopic ratios of the bulk rocks were measured at the CNR-IGG Pisa, Italy, using a Thermo Neptune Plus multi-collector-ICP-MS.Powdered rock sample was fused with ultrapure K 2 CO 3 and B was then leached into solution.The resulting solutions were purified using Amberlite IRA-743 boron-specific anion exchange resin and an AG 50W-X8 cation exchange column.The full analytical details are available in the Supplementary Methods of Iveson et al. (2021).
Three additional whole rock samples (TOL-16-01, TOL-16-04, and TOL-16-19) were analysed for major element, trace element, and B isotope data only (with no corresponding melt inclusions) and are included for reference in Supplementary Table 1.

Olivine phenocrysts
Olivine phenocrysts across all five samples span a relatively restricted range in major element composition.The most primitive olivine are from sample TOL-16-14 and have an average composition of Fo 76.8 ± 0.2 mol.%, whereas TOL-16-08 has the most evolved olivine populations, with an average composition of Fo 74.6 ± 0.4 mol.% (Table 1).Only two individual analyses from a single olivine in sample  show more magnesian compositions (≥Fo 80 ), potentially indicating a xenocrystic core (Supplementary Data A4).Most olivine are unzoned, or very weakly normally zoned, with core vs rim analyses generally within 1 mol.%Fo content within each crystal.
Consistent with the trends expected from fractional crystallisation of olivine ± clinopyroxene ± plagioclase ± Cr-spinel (e.g.Churikova et al., 2015a), there are broad decreases in Ni and Cr with decreasing Fo content (Supplementary Fig. 3A), whereas there are moderate increases in Sc, Zn, and Y. Lithium, Cu, and Co contents show little overall correlation with Fo mol.% (Supplementary Fig. 3B) but as with other trace elements, each cone sample tends to form its own cluster of data, despite the restricted range in major elements.Full trace element data for the host crystals are available in Supplementary Data A5.
Only a small number of trace element analyses were performed for the clinopyroxene (n = 4; Table 1) and all four crystals show similar trace element concentrations (Supplementary Data A5).

Post-entrapment corrections
The measured olivine-hosted melt inclusion compositions were corrected for post-entrapment crystallisation (PEC) using the Petrolog3 software (Danyushevsky & Plechov, 2011).Glasses were corrected for Fe-loss (Danyushevsky et al., 2000) and recalculated to equilibrium with the host phenocryst, using the Danyushevsky (2001) model calibrated for hydrous melt compositions.Corrections were performed at moderately oxidising conditions of QFM + 1.5 log units using the model of Borisov & Shapkin (1990), and the initial FeO concentration in the melt inclusions was specified as either 10.0 or 10.5 wt.%-close to the uncorrected maximum measured in the different melt inclusion suites and within 1 wt.% of the corresponding bulk-rock composition (Supplementary Table 1).The choice of the 'true' FeO concentration broadly affects the magnitude of the correction, i.e. the amount of olivine subtracted from, or added to, the melt inclusion composition, and given the limited petrographic evidence for PEC (e.g.fully glassy, spherical inclusions, with little/no faceting at the inclusion-host interface; Fig. 2) the choice of FeO content is consistent with this.The full uncorrected and PEC-corrected melt inclusion data are available in the Supplementary Data A6.On average, PEC corrections were small and required only 2.9 ± 2.1% olivine addition to the melt inclusion composition to achieve equilibrium with the host phenocryst.These generally small PEC corrections lead to similarly minor corrections to the melt inclusion trace element concentrations, via a 'correction factor' calculated as the corrected glass K 2 O concentration divided by the raw glass K 2 O concentration.
Conversely to the olivine-hosted melt inclusions, PEC corrections were not applied to the clinopyroxene-hosted MIs.Apparent Kd cpx/liq Fe-Mg values (Kd = (XFe/XMg) cpx / (XFe/XMg) melt ) calculated from cpx-melt inclusion pairs for each phenocryst are variable, ranging from 0.128 to 0.246, with an average of 0.188 ± 0.03 (Supplementary Fig. 4).However, 17 of the 42 melt inclusions show apparent Kd melts from Seguam and Fuego volcanoes, with very similar melt and phenocryst major element and H 2 O contents.The positive correlation between apparent Kd cpx/liq Fe-Mg and melt inclusion CaO concentration, along with the less spherical and more faceted melt inclusion shapes (Supplementary Fig. 4A), suggests that the most evolved cpx-hosted melt inclusions may have been subjected to moderate PEC (up to ca. 35%).However, since we do not use the cpx-melt inclusion pairs to explicitly reconstruct changes in equilibrium crystallisation conditions, in the following sections all cpx-hosted melt inclusions are discussed as uncorrected compositions, and in all figures, only cpx-hosted melt inclusions with Kd cpx/liq Fe-Mg > 0.20 are plotted.

Olivine-hosted melt inclusions
After PEC corrections, the olivine-hosted melt inclusion compositions straddle the basalt-basaltic andesite to trachybasaltbasaltic trachyandesite fields (Fig. 3A).SiO 2 and alkali concentrations for ZIM-16-07 are similar to basalts of the 'Naboko' event produced during the 2012-2013 eruption (Plechov et al., 2015), whereas the most primitive TOL-16-43 inclusions overlap with the more evolved groundmass glasses erupted from the 1975-1976 and other 0.8-7.8-kacones (Mironov & Portnyagin, 2018).Generally, as with the olivine phenocryst compositions, each cone forms its own cluster of melt compositions, though with some  overlap (e.g.Fig. 3B, inset).A plot of MgO/Al 2 O 3 vs K 2 O can be used to differentiate the high-Mg and high-Al suites of Tolbachik rocks (e.g.Flerov et al., 1984), confirming that the new melt inclusion samples belong to the high-Al group of melts, unlike high-Mg melt inclusions sampled from rare eruptions in the southern part of the massif, including some 0.8-7.8-kacones, and the 1941 and 1975-1976 eruptions (Fig. 3B).The middle-K rocks of the Plosky Tolbachik and Ostry Tolbachik stratocones (i.e.'Trend 1') fall on a separate lower-alkali trend compared with the melt inclusions from this study that are higher-K (Churikova et al., 2015b), and which also have FeO/MgO ratios that show a more tholeiitic affinity (FeO/MgO between 2-2.5 at SiO 2 = 50-53 wt.%).
Prior data for the high-K basaltic rocks at Tolbachik show that they have typical arc-like trace element signatures, with strong but variable enrichment in large ion lithophile elements and light rare earth elements (Churikova et al., 2015b).Relative to the high-Mg series, they are also strongly enriched in Rb, Ba, Li, B, Be, and in most incompatible trace elements (except Sr) including HFSE (Zr, Nb, Ta, Hf), Y, and REE (Churikova et al., 2015a).The trace element signatures of the melt inclusions from the cone samples in this study are consistent with these observations.Trace elements such as Li and B behave conservatively and correlate positively with indices of fractionation (i.e.K 2 O) and increase by approximately threefold relative to the higher-Mg melt inclusions from the 1941 eruption (Fig. 4).Consistent with the higher alkali contents, melt inclusions from the 2012-2013 eruption are generally more fractionated than the studied cone samples, and average 63.7 ± 2.3 μg/g B, relative to an average of 44.7 ± 5.1 μg/g B across the new samples in this study.
Trace element ratio binary diagrams indicate that the melt inclusions match the 'Trend 2' rocks of Churikova et al. (2015b).In particular, the new cone samples have trace element ratios (e.g.lower Sr/Sm and Sr/Y; Supplementary Fig. 5) that indicate a greater inf luence of plagioclase crystallisation, relative to 'Trend 1' rocks (comprising the lower parts of stratovolcanoes, most of the upper stratovolcano, and some dykes; Churikova et al., 2015b).Incompatible trace element ratios (Fig. 5) again show tight clusters for each cone, and Nb/Zr ratios that are similar across the range in K 2 O contents of the new cone samples, most of the stratocone samples, the 2012-2013 eruption melt inclusions, and the high-Mg 1941 eruption melt inclusions (Fig. 5B and C).Conversely, Rb/Ba shows a broad positive correlation with increasing K 2 O (Fig. 5F) and strongly overlaps with the 2012-2013 whole rock data, and other 'Trend 2' rocks (Fig. 5E).These data are discussed in greater detail in later sections.The full trace element data are available in Supplementary A6.
Melt inclusion CO 2 contents were not reconstructed for the formation of shrinkage bubbles (e.g.Moore et al., 2015;Wallace et al., 2015), and thus H 2 O-CO 2 concentrations measured in the glasses represent minimum values (e.g.Gaetani et al., 2012).However, the majority of melt inclusions lack shrinkage bubbles (Fig. 2B), and there are consistent differences between the volatile inventories of the different cone samples (Fig. 6A).The highest water contents (ca. 3 wt.%)are recorded by the most evolved sample ZIM-16-07, whereas the highest CO 2 contents (ca.1250 μg/g) are preserved in the most primitive TOL-16-43 melt inclusions.The highest calculated saturation pressures calculated from H 2 O-CO 2 solubility modelling (VESIcal; Iacovino et al., 2021 using MagmaSat;Ghiorso & Gualda, 2015) are ca.220 MPa, in melt inclusions from TOL-16-43.The three closely situated cones (Fig. 1B) all record very similar maximum H 2 O contents (ca. 2 wt.%) and CO 2 contents (100-400 μg/g) giving a cluster of minimum entrapment estimates between 50 and 110 MPa, similar to the volatile concentrations from the 2012-2013 eruption (Fig. 6A).In other Tolbachik magmas, the highest H 2 O contents (up to 5.3 wt.%) and CO 2 contents (>7000 μg/g) were measured in higher-Mg 1941 samples (Kamenetsky et al., 2017;Zelenski et al., 2022), yielding minimum entrapment pressures between 300 and 400 MPa, and up to Fig. 3. Olivine (PEC corrected)-and clinopyroxene (uncorrected, anhydrous basis)-hosted melt inclusion compositions and bulk rock compositions from the new cone samples in this study, along with literature data for several other Tolbachik eruptions (both bulk rock and melt inclusion samples).[A] Total alkali-silica plot showing the melt inclusions lie mostly in the trachybasalt and basaltic trachyandesite fields and generally with distinct clusters for each cone sample.[B] The MgO/Al 2 O 3 ratio vs K 2 O concentration for the melt inclusions and bulk rocks, with an expanded inset for the melt inclusions in this study and the 2012-2013 eruption inclusions.Compositional delineations after Flerov et al. (1984) and Portnyagin et al. (2015).Note the two grey squares with black crosses are bulk rock data from Portnyagin et al. (2015) for the same Peschanye Gorki cone (one high-Mg, one high-Al) as analysed in this study (TOL-16-43; grey square white cross).Literature data sources for melt inclusions and bulk rocks included for comparison are consistent in all subsequent figures and are as follows: Plosky and Ostry stratocones and feeder dykes (Churikova et al., 2015b); 0.8-7.8ka cones (Portnyagin et al., 2015;Mironov & Portnyagin, 2018); 1941 eruption (Churikova et al., 2015b;Kamenetsky et al., 2017);1975-1976eruption (Portnyagin et al., 2015;;Mironov and Portnyagin, 2018); 2012-2013 eruption (Churikova et al., 2015b;Plechov et al., 2015;Portnyagin et al., 2015;Volynets et al., 2015).As in all subsequent figures, only cpx-hosted melt inclusions with Kd cpx/liq Fe-Mg > 0.20 are plotted.
≥700 MPa if corrected for potential CO 2 -loss to post-entrapment shrinkage bubbles.
Sulphur contents of the melt inclusions correlate positively with MgO/Al 2 O 3 ratios and negatively with K 2 O (Fig. 6B), for both the high-Mg melt inclusions reported in the literature as well as four of the samples in this study, again with each cone tending to form a cluster.Maximum S contents up to ca. 1000 μg/g, however, are found in the ZIM-16-07 cone and do not correlate with the MgO/Al 2 O 3 ratio of the inclusions.These inclusions also show the largest range in S contents, from ca. 200 to >1000 μg/g, relative to the more clustered S contents of the other cones.In contrast to H 2 O, CO 2 , and S, both Cl and F are essentially immune from diffusive re-equilibration and shrinkage bubble formation during melt inclusion post-entrapment modification at magmatic temperatures (Bucholz et al., 2013).Average Cl contents across all melt inclusions are 585 ± 85 μg/g, and average F concentrations are 360 ± 33 μg/g.The Cl/Nb and F/Zr ratios of the inclusions correlate positively, with tightly clustered data for each cone (Fig. 6C).Application of the Webster et al. (2015) Cl-solubility model to the melt inclusion compositions at 1000 • C and 200 MPa storage conditions predicts maximum Cl solubilities (when in equilibrium with hydrosaline brine) ca. 2 wt.%Cl, significantly higher than actual measured Cl concentrations.Overall, melt inclusions from ZIM-16-07 record the highest H 2 O, S, and Cl/Nb and F/Zr values of the samples in this study.The more primitive 1941 eruption melt inclusions (Kamenetsky et al., 2017;Zelenski et al., 2022) have average Cl and F contents that are both higher than the high-Al melts in this study, 1247 ± 99 μg/g and 554 ± 58 μg/g, respectively, that testify to apparent volatile loss during melt evolution.The primitive 1941 eruption melt inclusions also record anomalously high S contents, up to ca. 11 500 μg/g, and evidence for co-existing magmatic sulphide globules with anhydrite in the quenched inclusion glasses (Zelenski et al., 2022).See later Discussion sections for further consideration of these data.

Clinopyroxene-hosted melt inclusions
The cpx-hosted melt inclusions from sample TOL-16-43 are more strongly fractionated than the olivine-hosted inclusions from the same sample, and the inclusions with Kd cpx/liq Fe-Mg > 0.20 are basaltic trachyandesites (Fig. 3A).Relative to the ol-hosted inclusions, they have approximately twice the K 2 O contents with only a moderate increase in SiO 2 (Fig. 3B).
Trace elements such as B (Fig. 4B) and Rb (Supplementary Fig. 6) show perfectly incompatible behaviour, increasing proportionally with K 2 O.However, Li appears to have been affected by postentrapment loss (Fig. 4A; Audétat et al., 2018).A comparison of the most primitive cpx-hosted inclusion (based on K 2 O concentration) relative to the average of the ol-hosted inclusions from TOL-16-43 (Supplementary Fig. 6) suggests ∼50% fractional crystallisation is required to explain the twofold increase in the most incompatible trace elements, assuming a bulk D mineral/melt ≈ 0 for both B and Rb.Incompatible trace element ratios, e.g.Nb/Zr and Rb/Ba, are broadly the same as the ol-hosted inclusions (Fig. 5C and F).
Relative to the ol-hosted inclusions, volatile species (H 2 O, CO 2 , and S) show lower concentrations (Fig. 6).H 2 O-CO 2 calculated saturation pressures cluster at ≤50 MPa, and S contents are strongly degassed, averaging ca.200 μg/g.Fluorine and Cl contents are slightly enriched (Supplementary Fig. 6) but not to the extent expected for ∼50% f luid-absent fractional crystallisation of an anhydrous mineral assemblage.Cl/Nb and F/Zr values are therefore the lowest of all samples measured (Fig. 6C) and lower than the ol-hosted inclusions from the same sample.(Kamenetsky et al., 2017).Note that modelling by Churikova et al. (2015b) has demonstrated that the composition of the 1941 eruption can be considered close to parental to the 'Trend 2' rock suite, to which these new melt inclusions belong.Both Li and B show fully incompatible mineral-melt partitioning behaviour, increasing proportionally with K 2 O from the high-Mg inclusions to the high-Al ones.The dashed lines show a simple fractional crystallisation model starting with an average 1941 eruption melt composition and where bulk D Li mineral/melt , D B mineral/melt , and D K2O mineral/melt = 0, with the stars representing 10% increments of fractionation.Lithium and B concentrations both suggest ∼50%-60% crystallisation is responsible for the evolution from the 1941 melt composition to the monogenetic cones of this study, increasing to ca. 70% for the 2012-2013 eruption.Note that clinopyroxene-hosted melt inclusions may be susceptible to post-entrapment Li modification (e.g.Audétat et al., 2018).three melt inclusions from ZIM-16-07 are isotopically lighter than δ 11 B = +1 .On average, clinopyroxene-hosted inclusions from TOL-16-43 are moderately heavier than olivine-hosted inclusions from the same sample (δ 11 B = +3.3± 0.2 vs δ 11 B = +2.1 ± 0.6 , respectively).These very tightly clustered inclusions are in stark contrast to the isotopic data for melt inclusions of the 1941 Tolbachik eruption (Gurenko, 2021) that show an apparently large range in δ 11 B, from −9.7 ± 2.3 to +5.4 ± 2.0 , and which are found in different eruptive sample types (lava vs lapilli and scoria).

B isotopes
The small range in δ 11 B values of the melt inclusions do not correlate strongly with any major or trace element signature, such as MgO/Al 2 O 3 , Ce/B, and Nb/Zr (Fig. 7), or volatile concentrations or ratios, including Cl/Nb and S/K 2 O. Consistent with other melt inclusion suites in Kamchatka (Iveson et al., 2021), the whole rock δ 11 B values generally ref lect the average δ 11 B of the corresponding melt inclusions from that sample, and range from δ 11 B = +1.65 to +2.44 in these whole rocks (Supplementary Table 1).

Chalcophile elements: Cu, Se, As, and Ag
Previous work has identified high Cu concentrations in various eruptive products of the Tolbachik volcanic fields (Vergasova et al., 2007;Vergasova & Filatov, 2012;Portnyagin et al., 2015;Mironov & Portnyagin, 2018), and rocks from the most recent 2012-2013 eruption were some of the most strongly Cu-enriched (up to ca. 450 μg/g; Volynets et al., 2015) of any measured samples.The new melt inclusions in this study have average Cu concentrations >220 μg/g, also reaching a maximum of ca.450 μg/g in two inclusions from TOL-16-08.Across the whole sample suite, only six of the inclusions from TOL-16-43 have Cu <100 μg/g, and which correlates positively with Ni for those inclusions, potentially indicating a role for PEC modification and/or diffusive re-equilibration of original Cu contents (e.g.Audétat et al., 2018).When data from all 'Trend 2' Tolbachik samples, excluding 'Trend 1' stratocones and feeder dykes, are taken together there is a generally negative relationship with MgO; high-Mg basalts average ca.125 μg/g Cu at >10 wt.% MgO, increasing to >400 μg/g for the high-Al suite with ∼4 wt.% MgO (Fig. 8).This contrasts with the generally positive relationship between Cu and MgO for other Kamchatka volcanic rocks, such as Bakening and Shiveluch, that show pseudo-linear downward trajectories with decreasing MgO concentrations, similarly to the MORB array (Fig. 8).Copper data for Gorely volcano (Gavrilenko et al., 2016) show a distinct downwards kink ca. 3 wt.%MgO, which coincides with the onset of magnetite crystallisation and concurrent decreases in V, FeO, and TiO 2 in those rocks (e.g.Jenner et al., 2010).The older 'Trend 1' Plosky and Ostry stratocones and associated rocks range from <40 to ca. >200 μg/g Cu, with no strong correlation with melt MgO.In summary, the rocks of 'Trend 2' (including the new monogenetic cone samples) from Churikova et al. (2015b) are consistently higher Cu than those from 'Trend 1' (stratovolcanoes and feeding dykes).
Selenium concentrations in the Tolbachik melt inclusions cluster around an average of 0.178 ± 0.026 μg/g and are therefore greater than the corresponding bulk-rock concentrations for all samples except TOL-16-08 (Supplementary Table 1).The Tolbachik melts span to higher Cu and lower Se compared with global MORB and produce a range in Cu/Se that shows the highest values currently measured in any glassy arc sample (Cu/Se > 1500; Fig. 9A).The only other comparable continental arc data for Se are from Antuco volcano, Chile, reported by Cox et al. (2019).Very high bulk-rock Cu/Se (spanning to >7000) in the Antuco samples results from bulk Se concentrations that are an order of magnitude lower than those measured in these Tolbachik melt inclusions (Se = 0.029 ± 0.016 μg/g; Cox et al., 2019).The Cu/Se values of the Tolbachik melt inclusions are higher than Hawaiian glass samples (Wieser et al., 2020) and backarc basin glasses prior to sulphide saturation (Jenner et al., 2010(Jenner et al., , 2015)), and result from generally elevated Cu contents relative to these other datasets.Conversely, the Tolbachik samples have average S/Se values of ca.3000 ± 600 that are comparable to oceanic plateau basalts and overlap primitive mantle values (Fig. 9B).The Tolbachik samples have not only high Cu/Se but also significantly higher As (2.42-3.03μg/g) compared with global MORB samples, Hawaiian glasses, and backarc basin glasses.Arsenic contents are comparable to the Antuco stratocone samples at equivalent MgO contents (ca.2-3 μg/g at ca. 4 wt.%MgO; Cox et al., 2019).
Downloaded from https://academic.oup.com/petrology/article/63/9/egac087/6677107 by guest on 10 February 2023  Volynets et al. (2013Volynets et al. ( , 2015) ) and Portnyagin et al. (2015) define two group of samples-one with low K 2 O and low Nb/Zr and another with high K 2 O and high Nb/Zr (green circled areas).These groups overlap some literature measurements (Tatsumi et al., 1995;Churikova et al., 2001;Portnyagin et al., 2007a;panel [A]) and were used as end-member compositions for binary mixing trajectories (e.g.Portnyagin et al., 2015; dotted curves with stars = 10% mixing increments) in the Tolbachik compositional suite (panels [A] and [D]).However, exactly the same samples of the 2012-2013 eruption analysed using ICP-MS by Volynets et al. (2013Volynets et al. ( , 2015) ) show the absence of a high Nb/Zr signature (Volynets et al., 2015), which is consistent with most ICP whole rock measurements by other authors (Turner et al., 1998;Churikova et al., 2001;Kamenetsky et al., 2017;this study;panel [B]).Moreover, neither ICP-MS measurements of whole rocks nor the measurements of melt inclusions from this and other studies (Portnyagin et al., 2007b, Nekrylov, 2015;Kamenetsky et al., 2017;this study;panel [C]) can corroborate the existence of the high Nb/Zr end-member from Portnyagin et al. (2015).Although ICP-MS data for whole rocks (panel [B]) and measurements of melt inclusions (panel [C]) show some variance, systematic changes in the Nb/Zr ratio vs K 2 O are not observed.Instead, the constant ratio of incompatible elements across a range in K 2 O contents indicates that fractional crystallisation is the strongest control on the evolution of the Tolbachik rock suite.The plot of Rb/Ba vs K 2 O (panel [E]) does not ref lect the mixing of magma but demonstrates their division into two trends in accordance with (Churikova et al., 2015b).The incorporation of Ba during abundant plagioclase crystallisation leads to an increase in the Rb/Ba ratio for rocks of 'Trend 2'.Note that, as can be seen in panel [F], melt inclusions of 'Trend 1' rocks were not measured in this study.Also note that the arrows indicating the general trend for fractional crystallisation vs melt mixing process are schematic.In panels [C] and [F], melt inclusions from the cones analysed in this study are colour-coded by sample and are consistent with all other figures.The outline colour differentiates the analytical technique; black = SIMS, orange = LA-ICP-MS.

Olivine-melt partitioning of Cu
An additional method to investigate the variability in melt Cu concentrations is to apply equilibrium crystal-melt partition coefficients to the host phenocryst compositions.Portnyagin et al. (2017) specifically analysed olivine-melt inclusion pairs from Tolbachik cone 'Mt.1004', and two other primitive Kamchatkan volcanoes.They report a very consistent D Cu ol/melt of 0.028 ± 0.009 (Supplementary Fig. 7A), across all three volcanoes, and across a

Fig. 6.
[A] H 2 O-CO 2 concentrations measured in the melt inclusion glasses in this study, and melt inclusions of the 2012-2013 eruption (Plechov et al., 2015).Saturation isobars for all melt inclusions from samples TOL-16-08 (red curves), TOL-16-43 (grey curves; dotted = cpx-hosted inclusions), and ZIM-16-07 (purple curves) were computed at the specified pressure using VESIcal (v.0.1.8;Iacovino et al., 2021;Wieser et al., 2021) with the MagmaSat model (Ghiorso & Gualda, 2015) at 1000 • C for olivine-hosted inclusions and 900 data for melt inclusions and host olivine from the 1941 Tolbachik eruption is also consistent with this olivine-melt partition coefficient (Supplementary Fig. 7A).Application of D Cu ol/melt = 0.028 to the Cu contents of the moderately more evolved olivine from the high-Al, high-K cones in this study (Supplementary Fig. 3B) would therefore yield equilibrium melt Cu contents ranging from ca. 320 μg/g to a maximum of ca.455 μg/g.This highest calculated melt Cu concentration is in accordance with two TOL-16-08 melt inclusions (Fig. 8) and some 2012-2013 bulk rocks.However, the majority of the melt inclusions cluster between ca.200 to 250 μg/g Cu, and if we calculate an apparent crystal-melt partition coefficient solely from our olivine-melt inclusion pairs, we derive a D Cu ol/melt of 0.048 ± 0.007 (Supplementary Fig. 7A).This value is approximately twice the value of Portnyagin et al. (2017)

Petrogenetic context within the Tolbachik volcanic massif
The new melt inclusions reported in this study are the first integrated dataset with volatile, trace element, chalcophile, and B isotope data for the high-Al, high-K suite of Tolbachik 'Trend 2' monogenetic scoria cones.These data also include the first geochemical results for the northern branch of the Tolbachik monogenetic fault zone, including the compositions of samples ZIM-16-07 and TOL-16-01 from the two northerly cones that were previously chemically undescribed (Churikova et al., 2015a).These data allow us to better constrain the petrogenetic context of these melts and their relationship to the other rock types produced by the Tolbachik volcanic massif.
Prior modelling has indicated that the systematic major and trace element differences between the high-Mg and high-Al suites, and middle-K ('Trend 1') vs high-K ('Trend 2') rocks, cannot be explained solely by variable extents of fractional crystallisation of the same parental magma at the same P-T-H 2 O-fO 2 conditions (Portnyagin et al., 2015;Churikova et al., 2015b).Sr-Nd isotope systematics have also shown that crustal contamination during magma ascent and storage plays only a very minor role in the composition of the 'Trend 1' rock suite (Churikova et al., 2015b).Thus, the rocks of the Tolbachik massif originate from the same, variably depleted mantle source, with the primary difference between the two trends ascribed to fractional crystallisation of relatively hydrous ('Trend 1' (stratovolcano samples)) vs relatively anhydrous ('Trend 2' (dykes and monogenetic cones)) parental Fig. 7. Variations in the melt inclusion and bulk rock B isotope ratio (δ 11 B) with Ce/B [A], and Nb/Zr ratios [B].Ce/B ratios are generally indicative of slab contributions to the mantle melting region, whereas Nb/Zr ratios are affected by source mantle enrichment/depletion, and/or degree of melt extraction.The limited variation in melt inclusion δ 11 B (1.9 ± 0.7 ) across all the samples analysed supports their derivation from a common, isotopically homogeneous parental melt.Overall, heavy isotope ratios are consistent with the addition of hydrous altered oceanic crust material (potentially in the form of down-dragged mantle wedge serpentinite) to an N-MORB depleted mantle composition.melts (Churikova et al., 2015b).This fractional crystallisationdominated regime can also explain the relationship between the melts of the 1975-1976 and 2012-2013 eruptions, where geochemical evidence suggests a common parent melt without requiring inputs from multiple sources and/or significant opensystem contamination (Volynets et al., 2015).
However, Portnyagin et al. (2015) suggest that the wider Tolbachik plumbing system is best described by an open-system model, whereby fractionation of the primary high-Mg moderate-K melts is accompanied by periodic mixing and hybridisation with fresh recharge melt.This open-system behaviour and extensive magma mixing is ref lected in the eruption of disequilibrium phenocryst assemblages (Flerov et al., 2015) and is consistent with the geophysical interpretations of Kugaenko & Volynets (2019).We find that the binary mixing models of Portnyagin et al. (2015) cannot satisfactorily explain the trend in Nb/Zr ratios vs K 2 O for our new melt inclusions or prior whole rock samples.Specifically, when only higher precision SIMS and solution/laser ablation ICP-MS data for Tolbachik whole rocks and melt inclusions are considered (Figs.5B, C, E, and F), the high Nb/Zr end-member purportedly erupted in 2012-2013 rocks is not apparent.Therefore, the lack of a systematic change in the Nb/Zr ratio vs K 2 O contents shown by the remaining filtered data does not support their derivation through binary mixing between a low K 2 O and low Nb/Zr end-member and a high K 2 O and high Nb/Zr end-member (Figs.5B, C).The ca. linear trend of Nb/Zr with K 2 O is better interpreted as incremental fractional crystallisation.
We have demonstrated that some trace elements, e.g.Li, B, behave conservatively, and changes in concentration of these elements between primitive and more evolved Tolbachik rocks can be reasonably modelled by simple fractionation (Fig. 4).Furthermore, the trends in Rb/Ba vs K 2 O (Fig. 5E and F) can also be broadly explained through extensive plagioclase fractionation, typifying 'Trend 2' rocks relative to 'Trend 1' rocks.Taking Rb and Ba contents similar to the average of the 1941 eruption melt inclusions (Kamenetsky et al., 2017)  The trace element data for the 2012-2013 eruption (Plechov et al., 2015; Supplementary Data A1) generally overlap with bulk rock data from the same eruption (Volynets et al., 2015), and they have Nb/Zr values effectively identical to the cone samples in this (again confirming that the apparently high Nb/Zr end-member in literature data is likely an analytical artefact) (Fig. 5C).Other incompatible trace element ratios Ba/Nb, Ce/B) in the 2012-2013 melt inclusions also overlap with the Holocene cone samples, and CaO/Al 2 O 3 and MgO/Al 2 O 3 (e.g.Fig. 3B inset) values are indistinguishable from the ZIM-16-07 cone sample specifically.The 2012-2013 inclusions have absolute K 2 O, Li, and B contents consistent with a further ∼25% fractional crystallisation of the average cone sample melt inclusion (Fig. 4), assuming bulk crystal-melt partition coefficients ca.0, for these components.The new cones in this study therefore ref lect a slightly less evolved high-Al, high-K melt composition, similar to that from which the 2012-2013 eruption melt was derived.The lower H 2 O, S (Fig. 6) and Cl/Nb (average ca.72 ± 15; Supplementary Data A1) in the 2012-2013 eruption samples are also evidence for additional f luid-saturated fractional crystallisation, relative to the new monogenetic cone data, potentially in a storage region at similar depths (e.g.Fig. 6A).However, we reiterate that neither the new melt inclusions in this study nor those sampled from the 2012-2013 eruption (Plechov et al., 2015) were corrected for H 2 O-CO 2 -S-loss to shrinkage bubbles (e.g.Venugopal et al., 2020) and therefore must ref lect minimum entrapment pressures.(Mironov & Portnyagin, 2018), one with Cu = 130 μg/g and one with Cu = 166 μg/g (see Section 5.4.1 for further details).The grey region shows the range of compositions produced by using each starting melt inclusion.Along the curves, F = percent melt remaining after fractional crystallisation.Representative D Cu mineral/melt values were specified for olivine (0.03), clinopyroxene (0.1), and plagioclase (0.02) after Liu et al. (2015) and Portnyagin et al. (2017).Note the sharp decrease in melt Cu contents recorded in the Manus backarc basin samples following magnetite-triggered sulphide saturation (Jenner et al., 2010) and a similar drop in Cu contents at similar MgO contents for Gorely volcano, apparently also triggered by magnetite saturation (Gavrilenko et al., 2016).These trends contrast with the gradual decrease in melt Cu contents recorded by MORB glasses and Bakening/Shiveluch volcanoes, indicating a role for sulphide liquid and MSS throughout melt evolution, respectively.Literature data sources as in Fig. 3 and Gorely (Gavrilenko et al., 2016); Bakening (Dorendorf et al., 2000;Iveson et al., 2021); Shiveluch (Ponomareva et al., 2015;Iveson et al., 2021); Manus backarc basin samples (Jenner et al., 2012(Jenner et al., , 2015)); MORB array (Jenner & O'Neill, 2012;Reekie et al., 2019).Note, cpx-hosted inclusions and the group of low-Cu ol-hosted inclusions from TOL-16-43 have not been plotted to reduce confusion, given the likely overprinting by PEC modification.

Boron isotope variability
Our new δ 11 B isotope measurements allow us to trace the involvement of a common magma source supplying the monogenetic cones.The B isotope ratio of a melt should not be fractionated during closed system evolution (consistent with incompatible crystal-melt partitioning; Fig. 4B; Kaliwoda et al., 2014; and incompatible f luid-melt partitioning in basaltic systems; Hervig et al., 2002).Given the limited inf luence of crustal contamination across the Tolbachik rock suites, the B isotope ratio of the melt inclusions should ref lect subducted-slab contributions to the mantle source during partial melting (de Hoog & Savov, 2018).The Ce/B ratios and heavy isotopic signatures (relative to MORB values; δ 11 B = −7.1 ± 0.9 ; Marschall et al., 2017) are consistent with the addition of hydrous altered oceanic crust material (potentially in the form of down-dragged mantle wedge serpentinite) to a N-MORB depleted mantle composition (Iveson et al., 2021) and are similar to the bulk rock value measured at Tolbachik by Ishikawa et al. (2001) (δ 11 B = +3.43 ± 0.32 ; measured using a Cs 2 BO 2 +graphite method and thermal-ionisation mass spectrometry).The data generally overlap with the melt inclusion records at other Kamchatka volcanic centres, e.g.Klyuchevskoy, Kamen, and Shiveluch, at equivalent Ce/B values (Iveson et al., 2021).
Thus, the overall small variations in melt δ 11 B (1.9 ± 0.7 ; Fig. 7) between 35 melt inclusions from distally located cones indicates that they were supplied from an isotopically homogeneous parental magma reservoir.These new data are also in agreement with the average δ 11 B reported by Gurenko (2021) for naturally quenched scoria-derived melt inclusions at Tolbachik, which show δ 11 B = +2.1 ± 2.7 .The overlap between the compositions of the 1941 eruptive products (Gurenko, 2021) and the Holocene cones of this study confirms that the B isotopic ratio of the parental Tolbachik melts appears to not be significantly modified during magma evolution from the middle-K, higher-Mg basalt (1941 eruption; Fig. 3B) to the high-K, higher-Al melts in this study.However, the most northern cone (ZIM-16-07; Fig. 1B) does record three inclusions that are isotopically lighter and more than two standard deviations outside of the average of the whole dataset.While trace element ratios of this cone are consistent with the other samples (Fig. 5), the relative volatile enrichment (S, H 2 O, Cl/Nb; Fig. 6) tentatively suggests that this could indicate contributions from a subsidiary melt source in this region.However, a more in-depth study of the δ 11 B composition of the primary mantle-derived high-Mg melts would help to decipher if any isotopic heterogeneity is present before later homogenisation Fig. 9. Chalcophile element trends as a function of MgO (wt.%) in the melt inclusions and bulk rocks of this study along with available literature data for melt inclusions, groundmass glasses, and bulk rocks from other Kamchatka volcanoes and other tectonic settings.[A] Cu/Se ratios; [B] S/Se ratios; [C] Cu/Ag ratios.Schematic arrows apply to panel [B] only and illustrate the broad trends expected by each process on a given data set.Data sources for literature samples included for comparison are as follows: Kīlauea and L ōi'hi glasses (Wieser et al., 2020); Reykjanes Ridge (Jenner et al., 2012); Oceanic plateau basalts (Reekie et al., 2019); discovery southern Mid-Atlantic ridge (Lissner et al., 2014); Manus and Lau backarc basin samples (Jenner et al., 2012(Jenner et al., , 2015)); MORB array (Jenner & O'Neill, 2012;Reekie et al., 2019); Kamchatka (Iveson et al., 2021); Antuco (Cox et al., 2019).Primitive mantle S/Se value taken from Palme & O'Neill (2014).Note that two bulk rock Tolbachik samples are not shown in panel [B] due to S concentrations below detection limit.
to produce the high-Al monogenetic cone eruptions across the whole volcanic field.
Unlike Gurenko (2021), we do not find a >15 range in B isotopes within the same sample.The extremely light values (δ 11 B ≤ −9 ) reported by Gurenko (2021) for experimentally rehomogenised melt inclusions from the 1941 eruption, and its correlation with δ 34 S only but no other major element, trace element, or volatile ratio appear to most likely have arisen due to sample processing effects: experimental rehomogenisation of melt inclusions in lavas vs naturally glassy scoria.Isotopically light values close to or lighter than MORB have only previously been described in Kamchatka at volcanic centres located in the far back-arc (Sredinny Ridge) and/or from tectonic regions with strongly reduced slab contributions (e.g.Bakening volcano) (Ishikawa et al., 2001;Iveson et al., 2021).We find the occurrence of light B at Tolbachik to be unlikely, given the abundant evidence for high slab-derived f luid f luxes beneath this portion of the Central Kamchatka Depression (Churikova et al., 2001;Portnyagin et al., 2007b), and therefore the potentially anomalous results of Gurenko (2021) will not be considered further.

Isolated magma batches after deeper mixing
The new melt inclusion suites from each individual monogenetic cone form tight clusters in major element, trace element, trace element ratio, and halogen compositional space, as do the trace elements within each olivine phenocryst population (Supplementary Fig. 3).However, they strongly overlap in their δ 11 B signatures, as discussed previously.We therefore interpret these trends to ref lect the independent evolution of each melt batch prior to eruption, in closely situated but relatively isolated storage regions.The spread in the eruption ages of the monogenetic cones also suggest that the physical conditions of melt storage have been relatively invariant over this period.A plumbing system with the requisite geometry beneath the Tolbachik volcanic massif has also been suggested by Kugaenko & Volynets (2019), based on quantitative seismic data.Kugaenko & Volynets (2019) demonstrate that the Tolbachik massif comprises a complicated magma plumbing system, with independent magma conduits and reservoirs, arranged in a subvertical and sublateral geometry.Additional geophysical data of Kugaenko et al. (2018) also demonstrate that there is a shallow magma reservoir at a depth of ca.4-8 km, developed along the rift system from the central to the southern part of the volcanic massif.This reservoir facilitates the migration and storage of magma, after melt homogenisation at deeper levels.This depth estimate derived from geophysical methods is also in good agreement with the minimum calculated volatile saturation isobars, specifically TOL-16-43 melt inclusions (Fig. 6A).
The major element composition of the new melt inclusions, in particular average K 2 O contents and correspondingly incompatible trace element concentrations, demonstrate only slight variations in the degree of fractionation for each cone (e.g.Figs.3B  inset, 4, and 5).Therefore, assuming all cones were fed from a common parental melt composition, the distinct Cl/Nb and F/Zr clusters (Fig. 6C) require variable loss of Cl and F from the melts prior to eruption.Similarly, the 1941 eruption melt inclusions (Kamenetsky et al., 2017;Zelenski et al., 2022) have average H 2 O, CO 2 , S, F, and Cl contents higher than all melt inclusions in this study (Cl/Nb = 912 ± 41, F/Zr = 9.9 ± 1.4), with extremely high CO 2 concentrations (up to ca. 7500 μg/g) measured in some inclusions (Zelenski et al., 2022).Both F and Cl dissolved in hydrous basaltic melts should be stable to very low pressures (0-5 MPa) (Webster et al., 2018;Thomas & Wood, 2021) and halogen concentrations should increase during f luid-absent crystal fractionation due to incompatible crystal-melt partitioning affinities (e.g.Webster et al., 2020, and the references therein).However, a small reduction in the Cl dissolution capacity of the melt would accompany melt evolution (i.e. a decrease in the molar [Al + Na + Ca + Mg / Si] ratio, Webster et al., 2018).The high-Mg melt inclusion data show a strong reduction in measured Cl contents as a function of increasing melt differentiation (Supplementary Fig. 8), supporting the sequestration of Cl to aqueous f luids during fractionation.However, this Cl-reduction trend f lattens (or the melt inclusion suites show a slight increase) with melt differentiation for the new cones in this study, along with the 2012-2013 eruption samples (Supplementary Fig. 8).The modelling results of Webster et al. (2020) indicate that much greater loss of Cl to a co-existing f luid phase would be expected in systems undergoing f luidsaturated fractionation and concurrent decompression, whereas f luid-saturated evolution at fixed pressures would lead to smaller reductions in melt Cl contents.Therefore, the Cl trends in the new melt inclusions may ref lect stalling of the magmas in the final shallow storage region immediately prior to eruption, combined with the effect of strong decreases in D Cl fluid/melt due to lower bulk Cl contents in the system, i.e. a reduction in the affinity for the f luid phase at lower total Cl (Webster et al., 2018(Webster et al., , 2020)).
The Peschanye Gorki cone (sample TOL-16-43) is one of the few monogenetic cones to erupt both high-Mg and high-Al melts (Portnyagin et al., 2015).Based on the total alkali contents and MgO/Al 2 O 3 ratios of the ol-hosted melt inclusions, they are the most primitive samples in this study (Fig. 3).Conversely, the cpxhosted inclusions from the same sample are the most evolved, with strongly incompatible trace element concentrations, e.g.B and Rb, suggesting ∼50% fractional crystallisation relative to the ol-hosted melts (Section 4.3.2;Supplementary Fig. 6).The average mol.Mg/(Mg + Fe TOT ) ratio of the clinopyroxene phenocrysts from TOL-16-43 is 0.74 ± 0.01, which is only slightly lower than the olivine phenocrysts from the same sample (0.76 ± 0.01) (Table 1).
Assuming that clinopyroxene and olivine with approximately the same Mg# will be in equilibrium in hydrous basaltic melts, this likely indicates co-crystallisation of both phenocryst phases.However, this is difficult to reconcile with the major and trace element concentrations, and volatile contents, of the melt inclusions within each phase.While H 2 O concentrations in a melt inclusion may reequilibrate rapidly with a degassed carrier melt across a clinopyroxene crystal (e.g.Lloyd et al., 2016), diffusive re-equilibration of Cl and F should be negligible.Both H 2 O and CO 2 contents of the cpx-hosted inclusions are degassed relative to ol-hosted inclusions, but Cl and F concentrations are moderately elevated.They are not, however, enriched to the concentrations expected from ∼50% f luid-absent crystal fractionation that is indicated by other incompatible trace elements.Therefore, as discussed in the previous section, some loss of Cl (± F) to a co-existing volatile phase seems necessary to explain the Cl/Nb and F/Zr ratios of the clinopyroxene melts.No ol-hosted melt inclusion overlaps with any cpx-hosted melt inclusion in H 2 O-CO 2 content (Fig. 6A) so if the cpx-hosted inclusions are indeed recording late-stage storage at <50 MPa pressure, the volatile contents of the ol-hosted inclusions have not re-equilibrated to these conditions.Thus, as discussed by Lloyd et al. (2016), two scenarios may explain the lack of overlap between compositions of ol-and cpxhosted melt inclusions erupted from the same cone.Either the phenocrysts are derived from unique melts that are crystallising at different depths and subsequently mixed prior to eruption or the phenocrysts are preferentially trapping samples of the same melt that is undergoing decompression-driven crystallisation, leading to a decrease in volatile concentrations (e.g.Fig. 6).On average the δ 11 B ratios of the cpx-hosted inclusions are moderately heavier (∼1 ) than the ol-hosted inclusions (Fig. 7) and may thus support the first scenario, indicating contributions from a different melt composition, albeit one with similar incompatible trace element ratios (Fig. 5C and F).Conversely, an additional consideration arises from the morphology of the melt inclusions within the clinopyroxene crystals.The more elongate and faceted 2-D shapes (e.g.Supplementary Fig. 4) suggest they were preferentially trapped along bands or melt channels within the phenocrysts and potentially experienced greater post-entrapment modification.Given the homogeneity in volatile contents, trace element ratios of the inclusions, and trace element contents of the host clinopyroxene phenocrysts, it is likely they were all trapped at the same stage during melt evolution.Speculatively, in the latter scenario, if a late-stage perturbation to the magma storage region triggered rapid crystallisation prior to eruption, it may also have led to boundary layer enrichment during melt inclusion entrapment, in particular more slowly diffusing major and trace elements (e.g.P 2 O 5 , Ba).Ultimately, if olivine and clinopyroxene were indeed cotectic, we suggest that late-stage syn-eruptive entrapment of more fractionated melts in clinopyroxene only (and followed by more PEC modification) is the most reasonable explanation for the compositional differences outlined in this section.

Chalcophile element systematics
Tolbachik and global data: discriminating sulphide-saturated vs sulphide-undersaturated evolution and mantle melting?
In order to interpret the behaviour of chalcophile elements in the Tolbachik monogenetic cone and fissure system, we compare and contrast our new data with other global magmatic trends.This allows us to distinguish three broad end-member scenarios, as defined by different datasets.
The first scenario is defined by the MORB array and exemplified by some of the Kamchatka volcanic suites (e.g.Bakening and Shiveluch).All three of these datasets demonstrate a downward trajectory of melt Cu contents with decreasing MgO (Fig. 8).This behaviour is consistent with the sulphide-saturated evolution of the magmas, as Cu is highly compatible in sulphide phases (e.g.Ripley et al., 2002;Jenner et al., 2010).In the case of MORB, the melts are relatively reduced, and contain sufficiently high S and Fe concentrations to become saturated with respect to a sulphide liquid during early stages of evolution (e.g.Mathez, 1976;Jenner, 2017).In the case of Bakening, recent modelling by Deng et al. (2022) confirms that low degrees of melting at low fO 2 conditions (ca.QFM ≤ +1) will leave residual mantle sulphides, and thus magmas will be sulphide-saturated at an early stage of magma differentiation in the crust.
The nature of the stable sulphide phases fractionating from MORB and convergent margin magmas can be determined by examining Cu/Ag systematics during magma differentiation (Fig. 9C).This is because sulphide melt shows equal compatibility for both Cu and Ag, whereas monosulfide solid solution (MSS) favours Cu over Ag (Jenner et al., 2010;Li & Audétat, 2012).MORB glasses show a decrease in Cu and constant Cu/Ag with decreasing MgO (Fig. 9c), which is consistent with sulphide melt fractionation.The concomitant decrease in melt Cu contents and Cu/Ag during differentiation of convergent margin magmas is instead indicative of exsolution of MSS (Jenner et al., 2010;Li & Audétat, 2012).This behaviour is shown by Antuco (Chile), and other Kamchatka samples, spanning a range of MgO contents.Indeed, recent modelling and interpretation of global Cu data reveals that most primitive arc magmas likely undergo pervasive sulphide-saturation, facilitated in particular by temperaturepressure-fO 2 conditions that favour Fe-depleting amphibole and occasionally garnet fractionation (Barber et al., 2021), lowering melt S solubility.This commonality is further supported by the ubiquitous preservation of sulphide inclusions (crystalline and polymineralic, with a Cu-rich and Cu-poor phase) in phenocrysts from Ecuadorian volcanic centres, spanning a wide compositional range from basalt to dacite (Georgatou et al., 2018).
The second broad scenario is where magmas are initially sulphide-undersaturated, which is likely because of either the high f O 2 of the source (S solubility increases with increasing f O 2 ) or the high melt fraction during melting (Deng et al., 2022, and the references therein).These sulphide under-saturated convergent margin magmas show an initial increase in Cu with decreasing MgO, until the melts eventually reach sulphide saturation (Jenner et al., 2010(Jenner et al., , 2015)).For example, backarc basin magmas from the Eastern Manus basin show an initial increase in Cu with decreasing MgO and constant Cu/Se until sulphide saturation is reached at ∼3 wt.% MgO.Sulphide saturation in these melts is triggered by the appearance of magnetite on the liquidus, which leads to a decrease in melt FeO content and a reduction in the sulphur species from sulphate to sulphide (Jenner et al., 2010).Consequently, Eastern Manus Basin magmas with <3 wt.% MgO show a decrease in Cu, Cu/Ag and Cu/Se with decreasing MgO (Fig. 9C), which is attributable to MSS fractionation.The effect of magnetite-triggered sulphide saturation is also displayed by the data of Gavrilenko et al. (2016) for Gorely volcano, Kamchatka (Fig. 8).
The final scenario is applicable to rocks of the Tolbachik monogenetic cone field (i.e.'Trend 2' rocks).Primitive high-Mg melts produced, for example during the 1975-1976 eruption and at other Holocene monogenetic cones (Portnyagin et al., 2015;Mironov & Portnyagin, 2018), show high initial Cu concentrations (ranging from 106 to 166 μg/g, with an average of 135 μg/g) for rocks with MgO >9 wt.% (Fig. 8).Concentrations spanning to 166 μg/g Cu are notably higher than average primitive MORB melts, which have only ca.115 μg/g Cu (Jenner, 2017).Modelling of mantle melting conducted by Mironov & Portnyagin (2018) showed that the elevated Cu (and S) concentrations in these primary Tolbachik melts may be achieved by melting of the same mantle source, but under more oxidising conditions for Tolbachik ( QFM +1.1-1.3 log units) relative to MORB ( QFM ≤0), and in both cases in the presence of residual sulphide.However, if melting was in the presence of sulphide, it is unclear why the 1975-1976 eruption and other Tobachik magmas were sulphide-undersaturated and show a pronounced increase in Cu with decreasing MgO (Fig. 8).In contrast, Deng et al. (2022) recently showed that such high Cu in primitive magmas requires mantle sulphide to be exhausted during melting.Importantly, Deng et al. (2022), used S and Cu concentrations that were applicable to melting of primitive (200 μg/g S and 20 μg/g Cu) through to depleted mantle (120 μg/g S and 17.6 μg/g Cu), whereas Mironov & Portnyagin (2018) used initial starting compositions that had significantly higher S (200-300 μg/g) and Cu (30 μg/g) than even the primitive mantle.Application of the Deng et al. (2022) model to the arc magmas of the Tolbachik massif, at representative conditions of fO 2 = QFM +1.1-1.5 and melting degrees of ca.8%-11% (e.g.Mironov & Portnyagin, 2018;Iveson et al., 2021), demonstrates that mantle melting would exhaust residual sulphide, generating sulphide-undersaturated Cu-rich melts.
Following melt generation at Tolbachik, Cu contents increase as melts become increasingly fractionated and reach maximum Cu contents of ca.450 μg/g at MgO ≈ 4 wt.% (Fig. 8).This unusual Cu-enrichment behaviour requires that Tolbachik melts do not saturate in a sulphide phase at any point during magmatic evolution (likely because the magmas do not reach magnetite saturation prior to eruption) and that bulk Cu partitioning behaviour remains incompatible throughout crystallisation.Notably, most of the whole rock and melt inclusion data for monogenetic cones that were analysed in this study show similar behaviour (Fig. 8).Importantly, unlike other Kamchatka data that show a decrease in Cu/Ag that is consistent with MSS fractionation, our new Tolbachik whole rock analyses overlap the range in Cu/Ag shown by the MORB array (Fig. 9C).This trend is consistent with our interpretation that the magmas were sulphide-(MSS) undersaturated during differentiation.Unfortunately, given the general lack of melts erupted from the Tolbachik monogenetic field with less than ca. 4 wt.%MgO, there is no information about the later evolutionary trends of these Cu-rich melts, as compared with other systems where post-magnetite-triggered sulphide saturation trends are discernible, e.g. the Manus basin (e.g.Fig. 8; Jenner et al., 2010).
Previous studies have demonstrated that increases in Cu with decreasing MgO can be attributed to degassing-related dissolution of sulphides during magma evolution (i.e.dissolution of the sulphides that were in contact with the melt, but not those that were hosted as inclusions in minerals) and consequent partitioning of Cu, Ag, and Se from the sulphides back into the melt (Reekie et al., 2019).In this scenario, Cu is not be expected to mimic perfectly incompatible behaviour because the proportion of sulphide that was fractionated from the melts (i.e. the proportion that remains in an underlying mush zone) and the proportion of sulphide that was trapped as inclusions in minerals would not be in contact with the degassing melt to be able to redissolve.As a hypothetical exercise, we have used Petrolog3 (Danyushevsky & Plechov, 2011) to simulate sulphide-free fractional crystallisation (ol + cpx + plag; 75% fractionation) and concurrent decompression (700 MPa to 300 MPa, δP/δT = 20 bars per 1 • C) at fO 2 = QFM + 1.5 log units (representative of the CKD in general; Nekrylov et al., 2018), for two starting melt compositions.We use two melt inclusions from Mironov & Portnyagin (2018); one with 130 μg/g Cu and one with 166 μg/g Cu, with similar MgO contents (ca.9.5 wt.%)In this model, bulk D Cu min/melt is close to 0, and Fig. 8 (thick dotted curves) shows that between ca.40 to 65% fractional crystallisation of these starting composition is required to explain the high Cu contents recorded by some 2012-2013 eruption samples and melt inclusions in this study.Depending on the initial Cu concentration for each sample, the simple model can broadly replicate the overall shape of the Cu vs MgO trend for the whole Tolbachik suite, and the spread in Cu contents in lower-Mg rocks can largely be enveloped by the two end-member model trends.Hence, the simplest explanation for the trend to extremely high Cu of the Tolbachik samples, together with their high initial Cu, is that the melts were sulphideundersaturated during differentiation and that loss of S during degassing simply made the magmas even more unlikely to reach sulphide saturation during further differentiation.Finally, we emphasise that the new data here cannot further discriminate the main processes controlling the difference in chalcophile metal behaviour between our 'Trend 2' rocks (monogenetic scoria-lava cones and fissures) and the 'Trend 1' rocks (Tolbachik central stratocone).The 'Trend 1' rocks show little evidence for anomalous Cu-enrichment (Fig. 8; Churikova et al., 2015b) and are similar to other closely situated volcanoes in the Central Kamchatka Depression (e.g.Klyuchevskoy and Kamen) and must indicate sulphide-saturated fractionation.The exact reason for the change in Cu systematics during tectonic evolution from older stratocone volcanism to more recent monogenetic cone and fissure eruptions remains enigmatic and is ultimately beyond the scope of this study.Speculatively, the increase in Cu concentrations between these time intervals, together with the shift from differentiation of sulphide-saturated (decreasing Cu with decreasing MgO) to sulphide-undersaturated melts (increase in Cu with decreasing MgO) might be suggestive of an increase in mantle f O 2 or a decrease in mantle fertility (required to exhaust sulphide during the most recent episodes of magmatism).

Tracing selenium behaviour
To our knowledge, there are currently no comparable Se data for arc-derived melt inclusion glasses and Se data for arc rocks in general are currently scarce.Our new data show that for all but one of the monogenetic cones (TOL-16-08), the wholerock Se concentrations are lower than those measured in the corresponding melt inclusions (Supplementary Table 1).Similarly, the whole rock data show higher Cu/Se compared with the melt inclusions (Fig. 9a).These data indicate the overall loss of Se from groundmass glasses, relative to melt inclusions, because the melt inclusions were shielded from external degassing.This trend is consistent with the other comparable datasets for arc bulk rocks measured by Cox et al. (2019) from Antuco volcano, Chile, and Iveson et al. (2021) for other Kamchatka volcanoes.These two comparable datasets both show very low S and Se concentrations relative to the MORB array, and S/Se ratios generally much lower than ca.1000.These trends are interpreted to ref lect low-pressure degassing of both these elements during eruption, and the preferential loss of S relative to Se.However, at Tolbachik it also appears that S and Se are degassing at the same rate during magmatic storage, as indicated by the high-Al melt inclusions: the S/Se of the inclusions overlap with primitive mantle values (Fig. 9B), which is consistent with the exhaustion of mantle sulphide during mantle melting (e.g.Reekie et al., 2019;Deng et al., 2022; this section).We caution that only further targeted Se analyses of higher-Mg and higher-S melt inclusions and groundmass glasses from primitive Tolbachik rocks will be able to fully decipher the behaviour of Se during melt evolution at the Tolbachik monogenetic cone field.

Sulphide immiscibility in the 1941 eruption and at cone 'Mt. 1004'
One final outstanding question is how the extremely sulphiderich olivine samples from the 1941 eruption and cone 'Mt.1004' are consistent with this petrogenetic model for sulphideabsent evolution of the Tolbachik massif.There is abundant evidence that some primitive melt batches in the Tolbachik monogenetic cone plumbing system are not sulphide-free.For example, silicate-sulphide liquid immiscibility in primitive melts, even under oxidising conditions (fO2 ≈ QFM + 1 to +1.5), has been investigated by Zelenski et al. (2018Zelenski et al. ( , 2022) ) and Kamenetsky et al. (2017).These authors describe densely concentrated immiscible sulphide droplets from separate monogenetic cones (the 1941 eruption and cone 'Mt.1004'), and which occur exclusively as inclusions in olivine phenocrysts, with none found in the coexisting groundmass (Zelenski et al., 2017).Such petrographic textures imply that any sulphides that were not trapped as isolated inclusions must therefore have dissolved during either the ascent of magmas or as a result of wider melt degassing of sulphur.
Several uncertainties surrounding their origins remain, for example why these samples specifically undergo localised, heterogeneous saturation in an immiscible sulphide liquid, despite the oxidising conditions (e.g.Kamenetsky et al., 2017;Zelenski et al., 2018Zelenski et al., , 2022)), and the origins of the 'excess' sulphur required to trigger such pervasive sulphide saturation (Zelenski et al., 2022).Furthermore, it is counterintuitive that the Cu contents of the olivine-hosted melt inclusions co-entrapped with sulphide inclusions from these samples are comparable to Cu contents in glasses from other sulphide-free Tolbachik melts at equivalent MgO contents (e.g.Fig. 8).We therefore conclude that while several models have explored the importance of sulphide melt resorption and chalcophile metal release back into an accompanying silicate melt (e.g.Kerr & Leitch, 2005;Reekie et al., 2019;Wieser et al., 2020, Chen et al., 2021;Holwell et al., 2022), we currently do not have the necessary data to decipher any contributions from this process to the compositions of these unusual Tolbachik samples.Thus, we suggest that the simpler model for sulphide-undersaturated evolution across the entirety of the compositional range of the new monogenetic cone melts analysed in this study is the most likely explanation for the trends discussed here.

CONCLUSIONS
Our new detailed geochemical investigation of melt inclusions from the Tolbachik monogenetic field contribute to our understanding of melt generation and evolution in this intra-arc rift zone, which is characterised by high melt production and eruption of magmas with contrasting major and trace element signatures.
We have demonstrated the following.
1) The major and trace element compositions of the new high-Al monogenetic cone melt inclusions are most satisfactorily explained through protracted fractional crystallisation processes, following derivation from a more primitive melt composition resembling those erupted during the 1941 or 1975-1976 eruptions.The 2012-2013 eruption melts are some of the most fractionated melts measured in the Tolbachik massif, but Nb/Zr values are indistinguishable from the new monogenetic cone samples.2) Tightly clustered δ 11 B ratios of the melt inclusions suggest that each cone was fed from a common, isotopically homogeneous parental source melt.Furthermore, δ 11 B ratios of naturally quenched melt inclusions from the 1941 eruption overlap the new data for the high-Al cones, supporting their derivation from a similar mantle source.3) Overall, the relatively heavy δ 11 B values (δ 11 B = 1.9 ± 0.7 ) are consistent with the position of the volcanic massif relative to the depth of the subducting slab and with contributions from hydrous altered oceanic crust material to a N-MORB depleted mantle wedge composition (e.g.Liu et al., 2020;Iveson et al., 2021).Bulk rock δ 11 B values also consistently ref lect the average δ 11 B of the corresponding melt inclusion suite.4) H 2 O-CO 2 contents indicate minimum storage pressures ranging from ca.50 to 200 MPa, and the southern Peschanye Gorki cone records some of the highest CO 2 contents measured in the high-Al suite (up to ca. 1200 μg/g).These volatile concentrations are consistent with geophysical evidence for a shallow magma reservoir at a depth of ca.4-8 km (e.g.Kugaenko et al., 2018).
5) Halogen (Cl/Nb and F/Zr) ratios are strongly correlated and clustered for each cone, demonstrating isolated melt evolution beneath each cone prior to eruption.Chlorine, F, and S abundances are lower in the high-Al melts relative to more primitive high-Mg melts, indicating a role for f luid-saturated evolution and volatile degassing during decompression, fractionation, and storage.6) A comparison of olivine-and clinopyroxene-hosted melt inclusions from the same sample show compositional differences requiring ∼50% fractional crystallisation, whereas host phenocrysts are apparently close to equilibrium with each other.This discrepancy likely results from the entrapment of more fractionated melt compositions pre-or syneruption, followed by more extensive post-entrapment modification of the clinopyroxene-hosted melt inclusions.7) The new monogenetic cone melt inclusions are not anomalously S-rich, but contain some of the highest Cu contents measured in any arc melt with 3-4 wt.% MgO.Despite evidence for early sulphide saturation in some anomalous Tolbachik melts, Cu contents continue to increase with melt evolution.We have shown that this trend requires melts to evolve under sulphide-undersaturated conditions across a range of MgO contents, relative to other global datasets and Kamchatka volcanoes where the formation of an immiscible MSS is clearly distinguished by fractionation of Cu/Ag ratios.8) Novel analysis of Se in the melt inclusions demonstrates that during shallow storage and degassing initial S and Se concentrations may be modified.However, S/Se values that overlap with primitive mantle signatures require mantle melting to exhaust residual sulphides, followed by melt degassing of both of these elements at the same rate from the melt.9) Further targeted Se and Ag analysis of the primitive high-Mg and high-S melt inclusions from both sulphide-bearing and sulphide-free monogenetic cone eruptions would be key in unravelling any contributions from sulphide melt resorption and chalcophile metal release in the Tolbachik plumbing system (e.g.Holwell et al., 2022).

Fig. 1 .
Fig. 1. [A] Simplified map of the Kamchatka peninsula showing the three main tectonic delineations of the Kamchatka arc; Eastern Volcanic Front (EVF), Central Kamchatka Depression (CKD), and Sredinny Ridge (SR).The approximate location of the Tolbachik Volcanic Massif and study area is given by the black square.[B] Enlarged view of the study area.The names and locations of main stratocone edifices are shown by the large white stars, with some of the Tolbachik-related eruptive products discussed in the text shown by smaller black stars.The locations and IDs of the monogenetic cones yielding melt inclusions and whole rocks analysed in this study are shown by the red circles and labels.Note that the samples from the 2012-2013 eruption with new data in this study (from Plechov et al., 2015) originate from the Naboko cone (lower star).For this figure and all following figure legends, the [−16] designation has been removed from the sample name for brevity (e.g.TOL-16-08 = TOL-08).A comprehensive map of the Tolbachik monogenetic lava field with ages and compositions of individual vents and lava f lows, modified after Flerov et al. (1984), is available in Churikova et al. (2015a).The base map image is an SRTM elevation model available from NASA/JPL/NIMA (https://photojournal.jpl.nasa.gov/catalog/PIA03374).

Fig. 2 .
Fig. 2. Ref lected light ([A], [C], [D]) and transmitted light ([B]) images of a single representative olivine phenocryst (TOL-16-10A-olivine21) before and after analysis by SIMS and LA-ICP-MS techniques.[A] Polished phenocryst prior to any analysis, showing two large, glassy and pseudo-spherical inclusions exposed at the surface.[B] Transmitted light image of [A] showing multiple unexposed inclusions of varying sizes, homogeneously distributed throughout the phenocryst and mostly free of post-entrapment shrinkage bubbles.[C] Phenocryst following analysis of the two large inclusions by SIMS for B isotopes and additional inclusions now exposed after further polishing.[D] Phenocryst following LA-ICP-MS analysis of both inclusions (one 50-μm spot, one 25-μm spot) and the host (110-μm spot).

Fig. 4 .
Fig. 4. Melt inclusion Li[A]  and B [B] concentrations vs K 2 O for the new monogenetic cones in this study, the 2012-2013 eruption, and the higher-Mg 1941 eruption(Kamenetsky et al., 2017).Note that modelling byChurikova et al. (2015b) has demonstrated that the composition of the 1941 eruption can be considered close to parental to the 'Trend 2' rock suite, to which these new melt inclusions belong.Both Li and B show fully incompatible mineral-melt partitioning behaviour, increasing proportionally with K 2 O from the high-Mg inclusions to the high-Al ones.The dashed lines show a simple fractional crystallisation model starting with an average 1941 eruption melt composition and where bulk D Li mineral/melt , D B mineral/melt , and D K2O mineral/melt = 0, with the stars representing 10% increments of fractionation.Lithium and B concentrations both suggest ∼50%-60% crystallisation is responsible for the evolution from the 1941 melt composition to the monogenetic cones of this study, increasing to ca. 70% for the 2012-2013 eruption.Note that clinopyroxene-hosted melt inclusions may be susceptible to post-entrapment Li modification (e.g.Audétat et al., 2018).

Fig. 5 .
Fig. 5. Variability in incompatible trace element ratio (Nb/Zr and Rb/Ba) vs K 2 O in Tolbachik rocks.Panels [A] and [D] are x-ray f luorescence (XRF) measurements of whole rocks, panels [B] and [E] are ICP-MS measurements of whole rocks, and panels [C] and [F] are measurements of melt inclusions by SIMS and/or LA-ICP-MS.Whole rock XRF data fromVolynets et al. (2013Volynets et al. ( , 2015) ) andPortnyagin et al. (2015) define two group of samples-one with low K 2 O and low Nb/Zr and another with high K 2 O and high Nb/Zr (green circled areas).These groups overlap some literature measurements(Tatsumi et al., 1995;Churikova et al., 2001;Portnyagin et al., 2007a; panel [A]) and were used as end-member compositions for binary mixing trajectories (e.g.Portnyagin et al., 2015;  dotted curves with stars = 10% mixing increments) in the Tolbachik compositional suite (panels [A] and[D]).However, exactly the same samples of the 2012-2013 eruption analysed using ICP-MS byVolynets et al. (2013Volynets et al. ( , 2015) ) show the absence of a high Nb/Zr signature(Volynets et al., 2015), which is consistent with most ICP whole rock measurements by other authors(Turner et al., 1998;Churikova et al., 2001;Kamenetsky et al., 2017; this study; panel [B]).Moreover, neither ICP-MS measurements of whole rocks nor the measurements of melt inclusions from this and other studies(Portnyagin et al., 2007b, Nekrylov, 2015;Kamenetsky et al., 2017; this study; panel [C]) can corroborate the existence of the high Nb/Zr end-member fromPortnyagin et al. (2015).Although ICP-MS data for whole rocks (panel[B]) and measurements of melt inclusions (panel [C]) show some variance, systematic changes in the Nb/Zr ratio vs K 2 O are not observed.Instead, the constant ratio of incompatible elements across a range in K 2 O contents indicates that fractional crystallisation is the strongest control on the evolution of the Tolbachik rock suite.The plot of Rb/Ba vs K 2 O (panel [E]) does not ref lect the mixing of magma but demonstrates their division into two trends in accordance with(Churikova et al., 2015b).The incorporation of Ba during abundant plagioclase crystallisation leads to an increase in the Rb/Ba ratio for rocks of 'Trend 2'.Note that, as can be seen in panel[F], melt inclusions of 'Trend 1' rocks were not measured in this study.Also note that the arrows indicating the general trend for fractional crystallisation vs melt mixing process are schematic.In panels [C] and [F], melt inclusions from the cones analysed in this study are colour-coded by sample and are consistent with all other figures.The outline colour differentiates the analytical technique; black = SIMS, orange = LA-ICP-MS.

Fig. 8 .
Fig.8.Variations in concentration vs MgO (wt.%) for a range of samples from the Tolbachik volcanic massif, along with literature data for other volcanic centres of the Kamchatka arc, Manus backarc basin, and the global MORB array.Almost all 'Trend 2' Tolbachik rocks are strongly enriched in Cu compared with other global datasets and show a generally negative relationship with MgO.However, 'Trend 1' rocks of the Plosky and Ostry stratocones are generally low Cu, similar to other Kamchatka rock suites.The thin dotted curves show best-fit evolutionary trends to the specified data series.The thick dotted curves show the results of Petrolog3 modelling of decompression fractional crystallisation for two starting melt inclusion compositions(Mironov & Portnyagin, 2018), one with Cu = 130 μg/g and one with Cu = 166 μg/g (see Section 5.4.1 for further details).The grey region shows the range of compositions produced by using each starting melt inclusion.Along the curves, F = percent melt remaining after fractional crystallisation.Representative D Cu mineral/melt values were specified for olivine (0.03), clinopyroxene (0.1), and plagioclase (0.02) afterLiu et al. (2015) andPortnyagin et al. (2017).Note the sharp decrease in melt Cu contents recorded in the Manus backarc basin samples following magnetite-triggered sulphide saturation(Jenner et al., 2010) and a similar drop in Cu contents at similar MgO contents for Gorely volcano, apparently also triggered by magnetite saturation(Gavrilenko et al., 2016).These trends contrast with the gradual decrease in melt Cu contents recorded by MORB glasses and Bakening/Shiveluch volcanoes, indicating a role for sulphide liquid and MSS throughout melt evolution, respectively.Literature data sources as in Fig.3and Gorely(Gavrilenko et al., 2016); Bakening(Dorendorf et al., 2000;Iveson et al., 2021); Shiveluch(Ponomareva et al., 2015;Iveson et al., 2021); Manus backarc basin samples(Jenner et al., 2012(Jenner et al., , 2015)); MORB array(Jenner & O'Neill, 2012;Reekie et al., 2019).Note, cpx-hosted inclusions and the group of low-Cu ol-hosted inclusions from TOL-16-43 have not been plotted to reduce confusion, given the likely overprinting by PEC modification.

Table 1 :
Average major element, trace element, and calculated end-member mineral compositions of the host phenocrysts to the melt inclusions from each monogenetic cone sample

Table 2 :
Average (and standard deviations)for the major, minor, trace, volatile, and B isotopic compositions of the melt inclusions from each monogenetic cone sample

Table 2 :
Continued Average compositions of the host phenocryst phase also given.Data are uncorrected for the effect of PEC on melt inclusion composition (see Results section for further details).