Earth and Planetary Science Letters Deciphering variable mantle sources and hydrous inputs to arc magmas in Kamchatka

The chemistry of primitive arc rocks provides a window into compositional variability in the mantle wedge, as well as slab-derived inputs to subduction-related magmatism. However, in the long-term cycling of elements between Earth’s internal and external reservoirs, a key unknown is the importance of retaining mobile elements within the subduction system, through subduction-related metasomatism of the mantle. To address these questions, we have analysed olivine-hosted melt inclusions and corresponding bulk rocks from the Kamchatka arc. Suites of melt inclusions record evidence for entrapment along melt mixing arrays during assembly of diverse parental magma compositions. Systematic variations in parental magma B/Zr, Nb/Zr, Ce/B, and δ 11 B are also apparent among the different eruptive centres studied. These element ratios constrain the nature of subduction-related metasomatism and provide evidence for ambient mantle heterogeneity and variable degrees of mantle melting. High Nb/Zr and low B/Zr in back-arc rocks indicate smaller degree melts, lower slab-derived inputs, but relatively enriched mantle compositions. Similarly, small monogenetic eruptive centres located away from the main stratocones also tend to erupt magmas with relatively lower slab contribution and overall smaller melting degrees. Conversely, arc-front compositions reﬂect greater slab contributions and larger degree melts of a more depleted ambient mantle. Across-arc variations in δ 11 B (ranging from ca . − 6 (cid:2) in the rear-arc and Sredinny Ridge to + 7 (cid:2) in the Central Kamchatka Depression) are generally consistent with variable addition of an isotopically heavy slab-derived component to a depleted MORB mantle composition. However, individual volcanic centres (e.g. Bakening volcano) show correlations between melt inclusion δ 11 B and other geochemical indicators (e.g. Cl/K 2 O, Ce/B) that require mixing between isotopically distinct melt batches that have undergone different extents of crustal evolution and degassing processes. Our results show that while melt inclusion volatile inventories are largely overprinted during shallower melt storage and aggregation, incompatible trace element ratios and B isotope compositions more faithfully trace initial mantle compositions and subduction inputs. Furthermore, we suggest that the signals of compositional heterogeneity generated in the sub-arc mantle by protracted metasomatism


Introduction
Magmas erupted at subduction zones represent a surface expression of the cycling of elements between Earth's internal and external reservoirs. Their compositions reflect ambient mantle chemistry, slab contributions, and later crustal processing (including polybaric fractionation, storage, degassing, and mixing). However, a key unknown is the importance of prior subductionrelated metasomatism in retaining volatile and fluid-mobile elements (FME) within the lithospheric mantle. Understanding the timescales and mechanisms of volatile and FME processing, transfer, retention, and recycling in subduction systems requires quantitative constraints on potential mantle and slab reservoirs.
Constraints on subduction zone slab and mantle processes are largely derived from whole-rock chemistry. These studies have demonstrated that the ambient sub-arc mantle chemistry (prior to subduction modification) can be more depleted than the depleted MORB mantle (DMM) or reflect mixing with enriched mantle similar to the sources of ocean-island basalts (OIB) (e.g. Pearce and Peate, 1995;Turner et al., 2017;and references therein). Trace element patterns observed in whole-rock samples are consistent with a varying ambient mantle that is overprinted by slab melts and altered oceanic crust-derived aqueous fluids (e.g. Wieser et al., 2019). However, while whole-rock data is a valuable tool, it cannot provide reliable information about the behaviour of volatile elements, and it obscures the extent to which samples of arc lavas reflect disequilibria. Rapidly quenched melt inclusions, on the other hand, can preserve volatile contents and sometimes provide a snapshot of real melt compositions from deeper in the magmatic system. Melt inclusions can therefore provide insight into subduction zone mantle processes that are otherwise inaccessible. In this study we will also demonstrate that melt inclusions can reveal deep crustal mixing processes that are not evident from wholerock data alone.
Boron concentrations and isotopes are a sensitive proxy for both overall fluid fluxes and the relative contributions of distinct hydrous slab sources to melting. The heavy 11 B isotope is selectively mobilised relative to 10 B during prograde metamorphism in the down-going slab (e.g. de Hoog and Savov, 2018;and references therein), such that the δ 11 B of slab-derived contributions becomes progressively lighter with increasing depths to the slab and distances from the trench. This generates arc-normal trends towards lower B and δ 11 B as observed in the Andean Central Volcanic Zone (Rosner et al., 2003) and Washington Cascades (Leeman et al., 2004). Boron isotopes also have the potential to discriminate fluid sources, as some studies have suggested that serpentinitederived fluids impart strongly elevated B and δ 11 B whereas altered oceanic crust and sediment inputs generally yield more modest enrichments due to prior forearc dehydration (Savov et al., 2007;Tonarini et al., 2011).
While whole-rock δ 11 B data consistently indicate the transfer of subduction zone fluids into the sub-arc mantle (e.g. Tonarini et al., 2001;Leeman et al., 2017), the whole-rock scale also conceals information about how small melt aliquots are generated and assembled into primitive arc magmas. In contrast, δ 11 B in melt inclusions has revealed compositionally variable subduction inputs (Gurenko et al., 2018); variations in fluid fluxes related to subduction geometry (Jones et al., 2014), source mantle compositions for OIBs (Walowski et al., 2019), and distinct subducted water sources (Cooper et al., 2020).
The Kamchatka arc represents a case study to examine the processes of subduction-related metasomatism across its long history of older subduction and 'cool' thermal regime, resulting from the subduction of old (>100 Ma) Pacific oceanic crust (e.g. Ruppert et al., 2007). Here, we combine trace element, volatile, and δ 11 B compositions of suites of melt inclusions, and their corresponding whole-rocks, to explore the controls on geochemical variability across the Kamchatka arc. We evaluate the relative importance of heterogeneity in the mantle and subduction component, and degree of mantle melting, for generating the compositional diversity of arc magmas. We show that melt inclusion arrays can preserve δ 11 B signatures indicating variable subduction inputs to a geochemically heterogeneous ambient mantle, despite later crustal magma assembly and mixing.

Geological background and sample settings
The Kamchatka arc lies at the junction of the Pacific plate and the Okhotsk and Bering plates (Fig. 1). The margin is tectonically complex, with along-and across-arc changes in slab dip and Moho depth; evidence for slab-tearing and a relict slab; and intra-arc rifting associated with anomalously high magma production (e.g. Avdeiko et al., 2007;Portnyagin and Manea, 2008;Hayes et al., 2018). Previous studies have identified a melting regime dominated by hydrous flux melting, with contributions from a variably depleted MORB mantle wedge, an OIB-type enriched mantle, slab-derived aqueous fluids, slab-derived melts, and sediment melt (Churikova et al., , 2007Portnyagin et al., 2007a;Simon et al., 2014). Integrated thermodynamic-geochemical modelling has also attempted to link variations in whole-rock δ 11 B (Ishikawa et al., 2001) to the thermal structure and punctuated dehydration processes in the subducting Pacific slab (Konrad-Schmolke et al., 2016). Our new melt inclusion and whole-rock samples derive from seven separate volcanic centres across the Kamchatka arc, and represent distinct mantle source domains and melting regimes (Fig. 1B).

Eastern Volcanic Front (EVF) and Rear Eastern Volcanic Front (REVF)
The EVF (Fig. 1A) represents a typical transitional arc, constituting a chain of large stratovolcanoes situated along a constant slab-depth contour. According to the model of Hayes et al. (2018), sub-arc slab depths along the EVF are 80-100 km, somewhat shallower than the global average. Magmatism at the EVF generally originates from moderate degrees of partial melting of a DMM source, fluxed by water-rich slab-derived fluids (Churikova et al., , 2007Portnyagin et al., 2007b). The EVF sample in this study is from the 'Medvezhka' cone; a monogenetic centre situated ca. 10 km from the edifices of both Avachinsky and Koryaksky (Supp. Fig.  DR1). Recent data on Holocene eruptions from EVF volcanoes suggests the Medvezhka cone may be compositionally more related to Koryaksky volcano, rather than Avachinsky (Krasheninnikov et al., 2020). The sample will be henceforth referred to as Medvezhka to avoid confusion regarding its petrogenesis.
The REVF consists of intermittent volcanism occurring up to ca. 100 km behind the EVF, and was once the active volcanic front (late Miocene, 7-5 Ma), before accretion of terranes to the eastern side of Kamchatka caused cessation and migration of the active -the names and locations of the volcanic centres with new melt inclusions and whole rocks analysed in this study (larger white dots), and those from literature studies with pre-existing melt inclusion and/or whole rock data included for comparison (smaller black dots); [C]depth contours to the top of the subducting slab (kilometres), taken from the subduction zone geometry model (Slab2) of Hayes et al. (2018). Base map image is an SRTM elevation model available from NASA/JPL/NIMA (https://photojournal .jpl .nasa .gov /catalog /PIA03374). trench to its current location (e.g. Avdeiko et al., 2007). REVF samples in this study include scoria from monogenetic cones related to Bakening stratovolcano. There has been little work on primitive tephra from Bakening, although studies of lavas and xenoliths have invoked mixing between variably depleted mantle sources, influenced by slab components (Dorendorf et al., 2000) and/or carbonate-rich slab-melts (Widom et al., 2003). Data from Farafonova Pad', a monogenetic cone located ca. 40 km northeast of Bakening (Fig. 1B), suggest that magmas in this region are generated by relatively low degrees of mantle partial melting and smaller slab contributions than the EVF (Portnyagin et al., 2007a).

Central Kamchatka Depression (CKD) and Northern Central Kamchatka Depression (NCKD)
The CKD (Fig. 1A) represents an active intra-arc rift zone with sub-arc slab depths ca. 20 to 40 km deeper than the EVF. Geochemical trends are attributed to gradual MORB mantle source depletion through protracted magmatism with the addition of distinct slab-derived fluids and/or hydrous slab-melts to the mantle source (Yogodzinski et al., 2001;Churikova et al., 2007Churikova et al., , 2013Portnyagin et al., 2007b). Mixing and/or assimilation of metasomatised lithospheric mantle by ascending primary melts has also been invoked under the CKD to explain the variable δ 18 O of wholerocks and olivine phenocrysts (Auer et al., 2009). Klyuchevskoy and Kamen volcanic centres are generally grouped together based on their close proximity (Supp. Fig. DR1) and evidence that some Kamen magmas share a common source with Klyuchevskoy volcano (Churikova et al., 2013). Samples from two satellite cones of Klyuchevskoy and a scoria deposit associated with Kamen volcano were analysed in this study (Supp. Fig. DR1).
The northern end of active volcanism in Kamchatka is marked by Shiveluch volcano (Fig. 1B). Volcanism here is thought to be related to the Aleutian transform fault zone (e.g. Yogodzinski et al., 2001;Portnyagin et al., 2005a) that marks the northerly termination of the actively subducting Pacific plate and shallowing of the slab dip (Fig. 1C). Samples from a monogenetic cone ca. 16 km from the edifice of Shiveluch and an unusual phlogopitebearing scoria deposit are included in this study. The Shisheisky volcanic complex (Fig. 1B) represents a field of Quaternary-age monogenetic cones ca. 50 km north of Shiveluch (Portnyagin et al., 2007b;Bryant et al., 2011).

Sredinny Ridge (SR)
The SR lies ca. 400 km behind the currently active trench ( Fig. 1A) with slab depths ca. 300-400 km (Fig. 1C), but once represented the active volcanic front around Pliocene times. Current volcanism results from melting of a fluid-metasomatised mantle aided by strong upwelling (Volynets et al., 2010) and includes both island-arc basalt (IAB) melts with typical subduction-like FMEenriched compositions, as well as within-plate type basalts (WPT) relatively enriched in HFSEs and REEs Volynets et al., 2010). We analysed samples from Ichinsky and Achtang volcanoes, previously studied by Churikova et al. (2001).

Experimental and analytical methods
Olivine crystals were hand-picked from 8 rapidly quenched tephra and 5 lava flows (Supp. Table 1). The tephra samples contained naturally glassy melt inclusions, whereas melt inclusions from the lavas were experimentally rehomogenised (Supplementary Methods). Olivine were mounted and polished, expos-ing melt inclusions that were then visualised by reflected light and backscattered electron imaging. Selected inclusions were analysed first for major and minor elements and S contents using the JEOL JXA-8230 electron microprobe (EPMA) at the University of Leeds. Analytical reproducibility was generally better than 3% relative for all major elements (Supp. Methods). Many inclusions were >50 μm width, so EPMA analyses were positioned to avoid electron beam-induced damage for later secondary ion mass spectrometry (SIMS) measurements (i.e. for H and Li). The inclusions were then analysed by SIMS at the NERC Edinburgh Ion Microprobe Facility for volatile elements, F, Cl, and trace elements (Cameca IMS-4f; with reproducibility generally better than 5% relative, Supp. Methods), and subsequently for 11 B/ 10 B (Cameca IMS-1270) following similar methods to Walowski et al. (2019). Full details of the melt inclusion preparation procedures, high-temperature homogenisation, analytical methods, and reproducibility/standard data are available in the Supplementary Methods.
Data were screened for contamination by the host phenocryst and imaged after SIMS analysis to confirm the exact analysis locations. We have not recalculated the melt inclusion compositions to correct for post-entrapment modification (PEC; see the Supplementary Methods for detailed discussion), but MI data interpretations are based on trace element ratios that are incompatible during olivine crystallisation (e.g., Nb, Zr, B, Ce) and 11 B/ 10 B ratios, which should be unaffected by post-entrapment crystallisation during ascent/storage.
Whole-rock major, minor, and trace element compositions were acquired through fused-bead XRF and solution ICP-MS analyses, conducted on fresh tephra and lava chips from the same samples as the melt inclusions. Whole-rock powders were also analysed for B concentrations and 11 B/ 10 B ratios by flux fusion and multicollector ICP-MS at CNR-IGG Pisa, Italy, following the procedures outlined in the Supplementary Methods.

Melt inclusions
A summary of the major, minor, trace element, volatile, and isotopic composition of melt inclusions for all volcanic centres is presented in Table 1, with the full dataset for all analysed inclusions available in Supplementary Table 1. Melt inclusions are hosted by olivine phenocrysts ranging in composition from Fo 89 to Fo 71 . Uncorrected melt-inclusions compositions of rear-arc volcanics from the REVF and SR range fall within or near the trachybasalt field, while inclusions from the Medvezhka, Kluchevskoy, and Kamen EVF/CKD samples are mostly basaltic andesite. Though most eruptions from Shiveluch are andesitic in composition, the inclusions in our samples range from trachybasalt to trachyandesite.
All melt inclusion samples show typical arc-like trace element patterns relative to N-MORB, with enrichment in B and LILEs (K, Rb, Ba), negative Nb anomalies, and decreasing abundances of increasingly compatible trace elements (Supp. Fig. DR2). Samples from the CKD are generally the most fractionated (Fig. 2) and have the highest B abundances, followed by the EVF, REVF, and then SR, but with some overlap. The single WPT basalt melt inclusion from Ichinsky (sample ICH-96-10) shows the most enriched trace element composition, and lowest B abundance of all samples analysed (Supp. Fig. DR2).
Measured H 2 O and CO 2 concentrations (Supp. Fig. DR3) show that, even when uncorrected for possible post-entrapment volatile loss, melt inclusions from Bakening are the most volatile-rich, and preserve H 2 O > 4 wt% and CO 2 > 2250 ppm, giving minimum entrapment pressures ≥5 kbar (Newman and Lowenstern, 2002). Other unhomogenised inclusions also preserve a range in volatile contents, yielding minimum pressure estimates up to ca. 5 kbar (Supp. Fig. DR3). These H 2 O and CO 2 contents represent minima due to the possibility of rehomogenisation during shallower magma storage, and likely do not preserve direct information about hydrous slab contributions to the mantle melting region (see later Discussion).
Melt inclusions show variable trace element and geochemical trends with SiO 2 (as an index of melt fractionation, Fig. 2). For some suites, variations in trace element concentration are consistent with up to ∼50% fractional crystallisation of olivine ± clinopyroxene + plagioclase (e.g. Klyuchevskoy-Kamen, see below). However, for Bakening, excluding more strongly fractionated inclusions with ≥ ∼53 wt% SiO 2 , the melt inclusion suite shows increasing K 2 O, B and δ 11 B with increasing SiO 2 , but decreasing S, Cl, Ca and Ti. Incompatible trace element concentrations are mostly flat or slightly decreasing with increasing SiO 2 , and Ce/B and Cl/K 2 O decrease systematically (Fig. 2). These features cannot be explained by magmatic fractionation and must therefore reflect mixing of different melt batches prior to eruption.
The Shiveluch melt inclusions are also not obviously related to one another by fractional crystallisation, showing decreasing Ce/B and Cl/K 2 O with increasing SiO 2 and an absence of systematic trace element trends. The inclusions are also not in equilibrium with their highly forsteritic host olivines (Table 1). Most inclusions have low Ba, LREE and LILE but a subset of three inclusions has elevated concentrations of most trace elements and LILE. Within the limits of data availability, this sub-group is associated with low δ 11 B (see also below).
For Klyuchevskoy-Kamen, the variations in trace element concentration can largely be explained by up to ∼50% fractional crystallisation of olivine ± clinopyroxene + plagioclase; incompatible elements are progressively enriched with decreasing MgO and Fo%. A subsidiary high-Al group of melt inclusions could be related to Al-rich magmas (Auer et al., 2009) but is not significantly different in trace element composition from the main population. There is also a subset of five inclusions in high-Fo olivines that typically have higher MgO contents and slightly lower Ce/B, K 2 O, and Rb (Supp. Table 1).
Melt inclusion δ 11 B varies widely across the arc, from ca. −7 in Bakening and SR samples, to ca. +7 in KGV samples. Both Bakening and the CKD also show large intra-sample ranges of ∼5 ; at Bakening this occurs across a small range in B concentrations (ca. 7.5-11 ppm, Table 1). Melt inclusion suites for Bakening, Shiveluch and Klyuchevskoy-Kamen generally show negative correlations of δ 11 B with Ce/B. For Bakening and Shiveluch, δ 11 B also correlates with major and trace element indices (see above), whereas at Klyuchevskoy-Kamen volcanoes, δ 11 B increases strongly with decreasing host Fo content (Supp. Fig. DR4).

Whole rocks
In general, the whole rocks are basaltic compositions with ∼48 to 54 wt.% SiO 2 and ∼5 to 9.6 wt.% MgO. Trace element concentrations are typically close to the average of the melt inclusions from that sample and δ 11 B values typically lie within the range defined by the melt inclusions. The full major, minor, trace, and δ 11 B compositions of the complementary whole rock samples are available in Supplementary Table 2.

Discussion
Our new dataset allows us to address three broad objectives: (i) interrogate the compositional variability of the sub-arc mantle in different regions of the Kamchatka arc, using key trace element ratios; (ii) decipher the nature and source of the hydrous subduction component(s), as recorded by the δ 11 B signatures of melt inclusions and whole rocks; and (iii) derive information about shallower magma aggregation and melt mixing processes at individ- Table 1 Average major, minor, and selected trace element compositions of the melt inclusions, along with volatile concentrations and boron isotope ratios.  ual volcanic centres. Along with our new data we have compiled published melt inclusion and/or whole rock analyses for 17 total volcanic centres across the entire Kamchatka arc that provide additional context (Fig. 1B).

Variable ambient mantle source composition and variable melting extents
Variation in the depletion/enrichment of the ambient mantle and variable mantle melting extents (F) can lead to similar trace element fractionations. In Kamchatka, variation in F and varia-tion in ambient mantle compositions have both been proposed (e.g. Churikova et al., 2001) and so must be considered concurrently. Middle-to heavy-REEs, Nb, and Zr are not readily mobilised from the slab (e.g. Carter et al., 2015), thus serve as useful tracers of ambient mantle and F. Because Nb is more incompatible than Zr during mantle melting (Kamber and Collerson, 2000), either low F or prior mantle enrichment can lead to elevated Nb/Zr. Ratios of middle to heavy rare-earth elements, such as Dy/Yb, are only strongly fractionated by F in the presence of residual garnet and are not strongly affected by ambient enrichment/depletion processes. To disentangle ambient mantle and melting effects for  Gale et al., 2013;and references therein). The models indicate that compositional trends in SR volcanoes, in particular Ichinsky, can be matched through variable, but low degrees of melting, and require residual garnet to explain the high Dy/Yb. Bakening whole-rock samples also record a similar petrogenesis but with a less enriched mantle. There is also strong evidence for two different mantle sources and variable melting degrees beneath the SR (cf. within-plate type and island-arc type basalt fields; Volynets et al., 2010). In contrast, the tight clustering in the CKD samples indicates consistent and larger degrees of partial melting. EVF volcanoes typically have similar moderate extents of melting to the CKD, except literature data for Gamchen and Avachinsky which show some evidence of slightly increased melting degree (see text for further discussion). Whole-rock data literature sources as follows: Mutnovsky  Kamchatka, we use a batch mantle melting model that accounts for varying mantle mineralogy (i.e. depleted vs. enriched; garnetbearing vs. garnet-free) under hydrous conditions at specified P-T conditions (Wieser et al., 2019).
Variations in both mantle composition and degree of melting are required to fully account for changes in melt composition across the arc. Fig. 3 shows whole-rock data from Kamchatka samples and mid-ocean ridge basalts alongside the models of mantle melting for enriched and depleted ambient mantle. Both Bakening (REVF) and Ichinsky (SR) datasets are consistent with variable but low degrees of melting (≤10 %), from a variably enriched mantle source, and both deviate from the MORB array towards higher Dy/Yb at a given Nb/Zr, consistent with residual garnet in the source (to be consistent with the data, the melting models must incorporate subsolidus garnet which is exhausted at low to moderate F). Furthermore, the Ichinsky WPT-basalts require very low melt fractions (∼2 %), from a strongly enriched OIB source, relative to the IAB-type basalts which require ≥5-10% melting.
Globally, Dy/Yb of arc volcanics correlate with local Moho depth (Turner and Langmuir, 2015), suggestive of deeper and lower degrees of melting where the overriding plate is thickest. This also appears to be true within Kamchatka, where both back-arc and Shisheisky Complex areas, which have deeper Mohos (Iwasaki et al., 2013), have elevated Dy/Yb. The lower Nb/Zr of the Shisheisky Complex is more consistent with a depleted ambient mantle, while the rear-arc volcanics require a more enriched source. Conversely, EVF and CKD generally all overlap each other, and the MORB array, requiring higher degrees of melting (≥10 %) of a DMM-like mantle source (Fig. 3). This is consistent with either higher bulk water contents or hotter and lower P melting conditions. The modelled melting trend also agrees well with the conclusions of Portnyagin et al. (2007a) and Bergal-Kuvikas et al. (2017) that magmatism at Klyuchevskoy requires 10-13% melting of a moderately depleted, garnet-free mantle source.
Thus, EVF and CKD volcanoes both require moderate degrees of melting of relatively depleted mantle, while all rear-arc data indicate low degrees of melting of more enriched mantle sources. Throughout the arc, variations in Nb/Zr likely reflect these com- bined effects. With this caveat in mind, we examine Nb/Zr and B/Zr in the melt inclusions to trace the subduction component added to the mantle source. B/Zr is not fractionated during fluid-absent mantle melting and/or olivine crystallisation (e.g. Marschall et al., 2017), and as B is highly enriched in serpentinised rocks (10 s to 100 s ppm B; Pabst et al., 2012) it is an excellent tracer of the addition of slab-derived material to the mantle. Thus, broad decreases in B/Zr coupled with increasing Nb/Zr (Fig. 4) indicate reduced subduction inputs and an increasingly enriched mantle source (coupled with lower melt fractions) towards the back-arc. Relatively higher B/Zr in EVF and CKD volcanoes supports previous interpretations of significantly higher subduction fluxes here (Churikova et al., 2007;Portnyagin et al., 2007a).
Both Shiveluch (sample SHIV-16-01; Table 1) and Medvezhka have a subset of melt inclusions with both low Nb/Zr and B/Zr similar to those of Farafonova Pad', indicating a subsidiary, depleted source with limited slab influence. Our Medvezhka melt inclusions are also significantly different from the high-Mg 'avachite' melts described by Portnyagin et al. (2005b) (Fig. 4). This finding suggests that small volume eruptive centres (e.g. the monogenetic Medvezhka cone sampled in this study, and the 'avachite' sample with uncertain provenance) preferentially sample melts from heterogeneous compositional domains of the mantle (Rawson et al., 2016). We suggest that such melts will otherwise tend to be homogenised during the repeated episodes of magmatism respon-sible for construction of the main stratocones. This highlights that, while the EVF as a whole would be expected to show higher B/Zr than recorded by AVX-16-02 (as seen in other literature studies of EVF volcanoes, e.g. Fig. 4), these samples provide finer details about more 'exotic' melts that may otherwise not observed in rocks of the main stratocones.

How is mantle source heterogeneity and evidence for variable melt fraction preserved?
We have established that primary melt composition is a product of both source heterogeneity and degree of mantle melting, especially at low degrees, and there is evidence for mixing between ambient mantle source compositions on a sample scale. Examining the distribution of Nb/Zr vs. Nb for each volcanic centre individually therefore allows us to infer the degree of mixing and geochemical homogenisation that occurs between melt generation and eruption.

EVF and REVF -Avachinsky and Bakening
Medvezhka melt inclusions are bimodal in Nb/Zr, with a group extending towards enriched Nb/Zr, and another overlapping N-MORB (Fig. 5A). Existing primitive whole-rocks for Avachinsky (and Koryaksky; Castellana, 1998) are highly consistent with the more depleted melt inclusions, as are the high-Mg avachite melt in-  Fig. 3. The preservation of significant Nb/Zr variability within individual suites of melt inclusions and whole rocks from the same volcano, and between volcanoes located in close proximity, suggests small length scales of source mantle compositional heterogeneity, superimposed on variations in melt fraction for some suites. While geochemical variation in e.g. Sredinny Ridge volcanoes [B], can be attributed to variable melting, deviation from the model lines in the suites of data for EVF, REVF, CKD, and NCKD volcanoes shows that ambient mantle variability most strongly controls Nb/Zr in these regions. Note that some panels have different axes limits to improve figure clarity. The two melt inclusions highlighted with squares in [D] are from the SHIV-16-01 high-K 2 O phlogopite-bearing tephra (see Section 5.2.4.). The high-Si Shiveluch melt inclusions have not been screened for the effects of late saturation in HFSE-bearing accessory phases (Tolstykh et al., 2000;Humphreys et al., 2008). Literature data sources as in Figs clusions (Portnyagin et al., 2005b). However, the more enriched end-member component recorded by the melt inclusions and corresponding bulk rock appears to be previously unrecognised at Avachinsky/Koryaksky.
Melts erupted from the Bakening monogenetic cones clearly display both a variably enriched source composition and evidence for variations in melt fraction (e.g. Figs. 3, 5A). The melt inclusions have Nb/Zr (Fig. 5A) similar to the more enriched mafic lavas (<55 wt.% SiO 2 , Dorendorf et al., 2000). In general, the average of each melt inclusion suite is similar to its whole rock composition, indicating that the whole rocks are generated by mixing multiple aliquots of variable source material. Conversely to Bakening, melt inclusions from Farafonova Pad' overlap the depleted Avachinsky and Koryaksky compositions ( Fig. 4; Fig. 5A). Taken together, the data suggest that magmas erupting in the EVF and REVF are formed from a mantle source region with contributions from both enriched and depleted domains, and monogenetic cones are more likely to sample a mantle source with lower overall slab contributions, relative to rocks erupted from the main central vents.

SR -Achtang and Ichinsky
The Ichinsky WPT basalts record the most enriched Nb/Zr values of the arc (Figs. 3; 4, Fig. 5B). Together, the Achtang and Ichinsky data record a continuum from moderately to very strongly enriched, with the overall trend closely matched by the melting model (Section 5.1; Fig. 5B). The presence of enriched mantle under the SR is consistent with the modelling results of Volynets et al. (2010) that require ∼30 to ∼55 wt.% contribution from a garnet lherzolite OIB-type source to account for the trace-element and REE signatures of mafic melts in this region.

CKD -Klyuchevskoy and Kamen
Almost all reported melt inclusions from the KGV overlap with MORB Nb/Zr compositions, with several samples extending towards the highly depleted S. MAR MORB field (Fig. 5C). This is consistent with melt extraction from a DMM source (Section 2.2). The horizontal spread towards higher Nb contents in our inclusions is likely attributable to melt fractionation (± small degrees of melt inclusion PEC), whereas the trend towards lower Nb in some literature data potentially results from 'over-heating' (i.e. re-melting the olivine host along the walls of the inclusion and thus 'diluting' the melt composition) during experimental homogenisation.
The combined melt inclusion dataset (this study and literature sources) covers ≥100% of the whole-rock Nb/Zr variability ( Fig. 5C), however the overall variation is much less compared with the EVF, SR, and NCKD samples. Given the large degrees of melting inferred for Klyuchevskoy magmas (Section 5.1.), the relatively reduced Nb/Zr variability likely indicates a locally more homogeneous mantle source region under this portion of the CKD.

NCKD -Shiveluch and Shisheisky complex
Both Shiveluch volcano and the neighbouring Shisheisky complex show large ranges in Nb/Zr (Fig. 5D). Our melt inclusions from Shiveluch (sample SHIV-16-09) are distinctly more enriched than  Fig. 3), but confirms that variations in both are required to explain the Kamchatka volcanic rocks. The Th/Yb ratio is a definitive tracer of sediment contributions to the mantle melting region, and the offset of the Kamchatka data towards much higher Th/Yb than the MORB array indicates that the majority of the arc requires a large transfer of slab-derived material. The whole-arc trend also shows a general concave-downwards curvature at high Nb/Yb and Th/Yb, almost overlapping the most enriched MORB glasses. This curvature and convergence of the arc samples with MORB glasses indicates that the WPT basalts of Ichinsky volcano require increasingly less slab inputs to the mantle source in this extreme back-arc region. Literature data sources as in Fig. 3. recent and historic andesitic tephra (Ponomareva et al., 2015), but other basaltic to andesitic tephra are similarly enriched (Tolstykh et al., 2015). The two SHIV-16-01 melt inclusions (Table 1) overlap with the main tephra suite. This unusual sample is taken from the high-K 2 O SHsp tephra layer (Volynets et al., 1997;Ponomareva et al., 2015) where phlogopite appears to be a primary liquidus phase (Krawczynski et al., 2019). The petrogenesis of this sample is not well understood, but our data suggest a strongly contrasting source relative to the other, typically middle-K 2 O, melts erupted at Shiveluch. As expected, the fractionated, high-Si rhyolitic melt inclusions plot towards higher Nb concentrations, but parallel to the Nb/Zr trend defined by the more primitive samples, and encompass a similar range in Nb/Zr (Fig. 5D).
The Shisheisky complex whole rocks also show a very wide range in Nb/Zr with most samples at the enriched end of the array, overlapping with the Shiveluch melt inclusions in this study. As with the EVF-REVF volcanoes, these trace element characteristics point to a highly heterogeneous mantle source region, coupled with some variation in degree of mantle melting (Figs. 3, 5D). Almost all Kamchatka whole-rocks have elevated Th/Yb relative to MORB (Fig. 6), indicating transfer of sediment and slab melts to the mantle melting region (e.g. Johnson and Plank, 1999) due to the extremely low partition coefficient of Th between aqueous fluids and eclogitic residue (e.g. Rustioni et al., 2019). The magnitude of the slab input decreases as the depth to the slab increases towards the back arc, reflected in the convergence towards the MORB array at high Nb/Yb for the most enriched WPT basalts of Ichinsky (Section 2.3, Fig. 6).
The incoming oceanic sediment pile along the Kurile-Kamchatka trench generally comprises diatomaceous silty clay and ooze, interlayered with volcanic ash (Plank, 2014). Siliceous oozes have a large range of [B] and δ 11 B (e.g. 86 to 148 ppm and −11.8 to +5.4 , Tonarini et al., 2011) but with no δ 11 B data available for Kamchatka trench sediments, we cannot quantify the relative δ 11 B contributions from siliceous vs. carbonate material. Given the generally light δ 11 B values of siliceous materials and their dominance in the Kamchatka sediments (Plank, 2014) we suggest that carbonates likely do not influence the overall transfer of 11 B to the sub-arc mantle in Kamchatka. Hence, the major controls on magma δ 11 B arise from ambient mantle isotopic heterogeneity and other hydrous slab contributions as outlined below.

Two component δ 11 B mixing model
Neither Ce/B nor δ 11 B are fractionated during fluid-absent partial melting and olivine crystallisation, reflected by identical Ce/B ratios and δ 11 B values for MORB and DMM (Marschall et al., 2017). The relationship between Ce/B and δ 11 B therefore provides a robust mechanism to assess the contributions of slab-derived components to the mantle melting region (de Hoog and Savov, 2018). To illustrate this subduction process, we calculate a simple, nonunique two-component mixing model (Fig. 7).
Model end-member compositions are the DMM (δ 11 B = −7.1 ± 0.9 and Ce/B = 10 ± 1; Marschall et al., 2017), an OIB-like mantle (La Réunion; δ 11 B ranges from −5.8 to −10.6 and Ce/B = 15.9 ± 1.6; Walowski et al., 2019), and a hypothetical 'subduction component' (where the combined contributions from AOC-derived fluids and melts, and sediment melts yield a δ 11 B range from +7 to −1 , and Ce/B = 0.1). The teal and grey fields in Fig. 7 illustrate the possible range in melt δ 11 B and Ce/B that could be produced through incremental addition of this variable subduction component to an average DMM and La Réunion mantle source.
Overall, the co-variation between increasingly heavy melt inclusion δ 11 B and decreasing Ce/B (Fig. 7) indicates that Kamchatka magmatism can be reproduced by addition of an isotopically heavy subduction component to a mantle source with both DMM and OIB-like compositions. However, greater detail is seen within individual melt inclusion suites, pointing to additional complexity in the interaction of slab materials with the ambient mantle wedge, and later melt mixing, as discussed below.

EVF and REVF
Melt inclusions (Fig. 7) and whole rocks (Fig. 8) from Medvezhka yield generally positive δ 11 B values, broadly consistent with isotopically heavy contributions to a mixture of typical DMM and a more enriched mantle composition. However, the relatively low B/Zr (Fig. 4) and high Ce/B in these samples are anomalous for arc front volcanism. We attribute this to the higher likelihood of satellite monogenetic cones sampling mantle regions which are less affected by slab-related metasomatism (Rawson et al., 2016) and caution that we do not consider this sample to be representative of the whole EVF (e.g. Krasheninnikov et al., 2020). positions suggest that the subduction component responsible for magmatism in this region has experienced significant prior B depletion. Such compositions could be generated by melting of a dehydrated slab (Tonarini et al., 2011), but this is inconsistent with the trace element and radiogenic isotope characteristics of the Bakening magmas (Dorendorf et al., 2000). Furthermore, the strong linear melt inclusion trend of increasing δ 11 B with decreasing Ce/B ( Fig. 7) and other geochemical parameters such as Cl/K 2 O (but not with host Fo content, Supp Fig. DR4), is not easily explained by the simple binary mixing model above. Combined with the flat trace element trends (Section 4) and the inability of crystallisation or degassing to modify δ 11 B (Hervig et al., 2002), we conclude that the isotopic variation in these Bakening samples is caused by mixing between end-member melts, generated from either isotopically distinct mantle sources (Fig. 7 inset), or from melts with variable slab contributions. We therefore infer that both B/Zr and δ 11 B signatures record information about the initial subduction input (slab contributions) and are preserved despite melt mixing in the crust (see Section 5.3.7 for further discussion).

SR
Melt inclusions from both Achtang and Ichinsky also show negative δ 11 B (Fig. 7) and lower overall slab contributions (Fig. 6). At whole-rocks (this study and literature data) for volcanoes of the Kamchatka arc.
Also shown are the average δ 11 B ± 2SD of MORB (Marschall et al., 2017) and the range in La Palma OIB data (Walowski et al., 2019). The whole rock δ 11 B values generally fall in the middle of the range of values recorded by the melt inclusions from that sample. This supports a petrogenetic model whereby the erupted melts are formed through the amalgamation of heterogeneous aliquots of liquid, that are derived from magma mixing processes and/or variable initial source regions.
Ichinsky, there is also no significant difference in δ 11 B between the WPT melt inclusion (Ce/B >10) and the more typical IAB samples (Ce/B ≤6), and little correlation between δ 11 B and Ce/B, Cl/K 2 O, or host olivine forsterite contents (Supp. Fig. DR4). The negative δ 11 B of these melts therefore may reflect the ambient mantle under the SR following protracted periods of earlier subduction metasomatism while the region was on the active volcanic front (Section 2.3.). Active volcanism in the SR is interpreted to result largely from decompression re-melting of fluid-metasomatised mantle (Volynets et al., 2010), and any current slab-derived B sources are expected to be strongly 11 B depleted. Recent data (Sugden et al., 2020) . 7; Ishikawa et al., 2001) and a typical DMM source.
However, we suggest that the significant variability in δ 11 B of melt inclusions from Klyuchevskoy-Kamen volcanoes again reflects melt mixing processes under the CKD, potentially resulting from interaction of ascending primitive melts with intensely metasomatised lithospheric mantle (e.g. Auer et al., 2009). The CKD melt inclusions show a negative correlation between δ 11 B and Cl/K 2 O, but unlike Bakening, higher δ 11 B also correlates strongly with increasingly evolved olivine compositions at Fo ≤80 mol.% (Fig. Supp. Fig.   DR4). We interpret this to reflect melt mixing and crystal recycling between variably evolved melt bodies under the KGV, as illustrated by the bi-modal olivine population, presence of a subsidiary high-Al melt, and variable δ 18 O signatures of olivine phenocrysts (Auer et al., 2009).
Given the evidence for high B fluxes in this region (Fig. 4), it is unexpected that the Shiveluch melt inclusions are isotopically light relative to inclusions from the other CKD volcanic centres. Tentatively, this could be related to the shallowing of the slab dip under the CKD, from ∼45 • under Tolbachik (Fig. 1C) to ∼30 • under Shiveluch (e.g. Portnyagin and Manea, 2008), leading to earlier fractionation of 11 B into the liberated fluid phase prior to reaching the melting region under Shiveluch (de Hoog and Savov, 2018). Although previous studies have suggested a component of true slab melting to explain NCKD magma geochemistry (Yogodzinski et al., 2001), comparative δ 11 B data for adakites (sensu stricto) is needed to test this hypothesis.
A final complexity for Shiveluch is that the predominantly andesitic melt inclusions in our primitive basalt sample are hosted by high Fo olivine crystals (Table 1; Supp. Fig. DR4). The melt inclusions could represent mixing of a small amount of andesitic material into a more primitive basalt; the interaction of primitive basalt with harzburgitic mantle (Volynets et al., 1997); or more evolved melts trapped in mantle-derived antecrysts during rapid ascent and interaction with shallowly stored magma (Gordeychik et al., 2018). Regardless of the exact mechanism of entrapment, the correlations between δ 11 B, Ce/B (Fig. 7), and Cl/K 2 O ratios (Supp. Fig. DR4) at Shiveluch volcano likely reflect mixing between isotopically heterogeneous melt batches that are preserved during later shallow magma storage and/or mixing events, as described below.

Melt inclusions preserve melt mixing arrays
The presence of diverse melt compositions trapped within olivine crystals from a single hand sample attests to the assembly of arc magmas through the aggregation of multiple small aliquots of melts at depth within complex plumbing systems. The range in melt inclusion trace element and isotope compositions, e.g. particularly for Bakening, Shiveluch and the CKD, thus reflect melt mixing arrays. Both Bakening and Shiveluch record a more primitive end-member, as well as a more evolved end-member that may have been affected by crustal assimilation and fractional crystallisation during magma storage. There are two possible scenarios to account for the generation of the different melt end-members: i) an isotopically heterogeneous ambient mantle + an homogeneous subduction input, or ii) homogeneous ambient mantle + an isotopically variable subduction input, as outlined in the previous sections and illustrated by the model (Fig. 7).
Unfortunately, the trace element and B isotopic composition of the Kamchatka deep crust is almost entirely unconstrained, and it is not possible to differentiate the effects of addition of an evolved, degassed, and stronger slab signature melt to a more primitive melt from the effects of pure crustal assimilation. We infer that mixing took place during ascent, between later-produced more primitive melts with reduced slab input (lower B/Zr, lighter δ 11 B) and earlier, shallowly stored degassed melts. Thus, the range in volatile contents preserved in melt inclusion arrays from the same sample are best interpreted as resulting from the pressure at which more evolved high-Si magmas are stored, and not the initial input from the subducting slab. Such a scenario is also supported by the fact that both Bakening and Shiveluch central cones erupt dacite-rhyolite compositions, and have been shown to evolve in open-systems with repeated episodes of mixing and homogenisation (e.g. Dorendorf et al., 2000;Humphreys et al., 2006b).
Regardless of the exact petrogenesis of the different melt endmembers, mixing of isotopically distinct batches is likely a process common to other regions of the Kamchatka arc, and other subduction zones globally. Nb/Zr must record the ambient mantle heterogeneity and initial degree of melting (Section 5.1), whereas the B/Zr and B isotope ratios reflect the initial slab contributions.

Conclusions
Using key trace element and isotopic ratios of melt inclusion suites alongside corresponding bulk rocks we have been able to identify the combined effects of ambient mantle composition, mantle melting degree, the nature of slab-derived contributions, and later melt mixing processes, on the composition of Kamchatka magmas. Broadly, higher slab-derived contributions to depleted mantle sources, coupled with higher degrees of melting in the EVF and CKD are reflected in high B/Zr, low to moderate Nb/Zr, and heavy δ 11 B signatures. Conversely, 11 B-depleted slab-derived contributions to lower degree melts of more enriched mantle compositions in the REVF are recorded by larger ranges in Nb/Zr and lighter δ 11 B. Current volcanism in the extreme back-arc SR results from decompression melting of a fluid-metasomatised source. Therefore, the δ 11 B of the SR volcanoes, and likely other global volcanic systems lacking direct fluid input (i.e. collision zones; Sugden et al., 2020), may reflect older geodynamic regimes and earlier metasomatism of the mantle when the region was on the active arc front (e.g. Leeman, 2020).
Our modelling of the combined δ 11 B and trace element dataset across this new suite of melt inclusions indicates that isotopically distinct melt batches could be produced by both distinct hydrous inputs and a variable mantle wedge, and whole rock δ 11 B values typically lie within the average defined by the melt inclusions. These melts are stored at different depths in the respective plumbing systems and undergo typical magmatic differentiation and degassing (± crustal contamination). The erupted magmas undoubtedly record melt mixing arrays, but Nb/Zr ratios retain information about source mantle enrichment/depletion and melting degree, while B/Zr and δ 11 B signatures reflect the hydrous subduction contribution. Our approach highlights the advantages of combining melt inclusion and whole-rock data to investigate small length scale geochemical heterogeneity, which may otherwise be obscured by the accumulative nature of bulk rock samples. Despite later crustal evolution and degassing overprints on the melt inclusion volatile contents, we show that the preservation of varied trace element and isotopic compositions at a single volcanic centre demonstrates that melt inclusions and whole-rock samples record the assembly of diverse melt batches from proximal mantle regions (often sampled by satellite monogenetic cones) in complex transcrustal plumbing systems, a process likely common to global arc volcanism.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.