Whole-rock compositions of the MSZK xenoliths (Supplementary Table S11) were determined based on mass balance calculations using the xenolith modal composition and the average mineral compositions.

The Group IA and IB single-lithology harzburgite xenoliths, and harzburgitic domains of the Group II xenoliths have homogenous whole-rock compositions (Supplementary Figure S4). In contrast, the dunite domain of xenolith MSZK1301D has higher MgO and lower SiO₂ and Al₂O₃ contents compared to the harzburgites. The Group IB olivine-orthopyroxenite, and Group II orthopyroxenite and websterite domains have lower MgO, but higher SiO₂, TiO₂, and Al₂O₃ contents with respect to the harzburgites (Supplementary Figure S4). The geochemical trend defined by the harzburgites differs from those of peridotite residue after batch and fractional melting (Supplementary Figure S4).

Whole-rock REY contents of the MSZK xenoliths were estimated by mass balance calculations and shown in Supplementary Table S11. Olivine, spinel and orthopyroxene compositions were estimated from clinopyroxene concentrations and mineral/cpx partition relationships using D^mineral/melt values of Ionov et al. (2002) and references therein. Yttrium results were calculated from Yb data based on their distribution in primitive mantle (Y concentration is ∼9.75 times more than Yb concentration) (McDonough and Sun, 1995).

According to the results, the websterite domains and the amphibole-phlogopite vein from Group II xenoliths have the highest REY concentrations. These domains show convex-upward pattern with a negative Eu anomaly (Supplementary Figure S5). The orthopyroxenite and harzburgite lithologies of Group II xenoliths have lower REY concentrations. The patterns of these domains are similar to those depicted by the websterites and the amphibole-phlogopite vein with the exception of the orthopyroxenite in xenolith MSZK1306B, which shows depletion in LREE (Supplementary Figure S5). The Group I xenoliths show flat or slightly LREE-enriched characters (Supplementary Figure S5). The slightly higher LREE and HREE contents compared to elements from Sm to Tb appear in multiple xenoliths, resulting in U-shaped patterns.

Group IA includes four single-lithology harzburgite xenoliths (MSZK0701; MSZK1301A; MSZK1301B; MSZK1310) and dunite domain of the composite xenolith MSZK1301D (Table 1).

The formation of orthopyroxene-rich mantle volumes is often explained by melt extraction and accompanying mantle depletion (e.g., Dick et al., 1984; Puziewicz et al., 2015). Alternatively, orthopyroxene enrichment can be a result of a melt-rock reaction between the lherzolitic wall rock and silica-rich agents leading to reactive harzburgites (e.g., Kelemen et al., 1992, 1998; Smith et al., 1999).

The orthopyroxene modes can be a powerful tool used to explain the origin of the harzburgites. This is because the residual orthopyroxene enrichment triggered by mantle melting is relative, driven by the consumption of clinopyroxene and spinel to a lesser extent (e.g., Hirschmann et al., 1998). In contrast, newly-formed orthopyroxene grains appear in reactive harzburgites at the expense of clinopyroxene (e.g., Zheng et al., 2001) or olivine (e.g., Arai et al., 2003, 2004). In the studied MSZK harzburgites, orthopyroxene contents are high and range between 12–52 vol% (with an average of 22 vol%) (Table 1.), which rather implies reactive origin. Such orthopyroxene-rich lithologies along with subordinate clinopyroxenes (≤4.5 vol%) require high-degree partial melting (>25%), which is irreconcilable with the proposed Cenozoic evolution of the CPR (Szabó et al., 2004). The reactive origin of harzburgites is further supported by the bulk rock major element trends (e.g., positive correlation between MgO vs. Na₂O and MgO vs. FeO) of the MSZK harzburgites differing significantly from the trends of partial melting of Niu (1997) (Supplementary Figure S4). Similarly, the variably enriched trace element patterns of bulk rocks (Supplementary Figure S5) and clinopyroxenes (Figures 6D–F) are irreconcilable with mantle depletion (Niu, 2004) as depleted mantle rocks are typically strongly diluted of highly incompatible trace elements.

The presence of numerous tiny (10–50 µm), elongated or rounded spinel inclusions hosted by orthopyroxenes in the MSZK harzburgites (Figures 3F,G) also cannot be explained by depletion via partial melting, because it never leads to spinel crystallization. Additionally, these inclusions are unlikely to form during mineral breakdown reaction (e.g., garnet) or exsolution because other phases (e.g., clinopyroxene) or pseudomorphs of the decomposing phase would also be expected to appear in such scenario (Falus et al., 2000). However, the metasomatic effect of silica-enriched fluids may play a key role in forming symplectitic intergrowth of orthopyroxene and spinel as proposed by Johan et al. (2017). Also note that the high amounts of Al₂O₃ (>5 wt%) and Cr₂O₃ (>5 wt%), which is required for subsolidus exsolution of spinel inclusions, is impossible to be incorporated in orthopyroxenes. Spinel inclusions, although less frequently, were also described in olivines of harzburgites (Figure 3G), similarly as in reactive harzburgites from Kamchatka (Arai et al., 2003) and the Philippines (Arai et al., 2004).

All prior considerations suggest that the MSZK harzburgites have a reactive origin. This orthopyroxene enrichment fits the definition of stealth metasomatism (O’Reilly and Griffin, 2013). Since there are no microstructures showing reaction (e.g., secondary orthopyroxene appearing at the expense of olivine or clinopyroxene) in the MSZK xenoliths, it is challenging to determine the precursor rock composition. Nevertheless, the MSZK harzburgites are among the BBHVF xenoliths with the highest olivine contents (dominantly >70 vol%; Figure 2). Thus, the dissolution of clinopyroxene may occur at larger amounts than that of olivine. This potentially implies a lherzolitic precursor, a lithology rare in MSZK xenoliths (Créon et al., 2017b), but frequent in other BBHVF localities (Figure 2). The possible reason why the metasomatism leaves no precursor lherzolites in MSZK locality is its pervasive style.

The only non-harzburgite lithology belonging to the Group IA is the dunite domain of the composite xenolith MSZK1301D (Table 1). It is adjacent to an olivine-orthopyroxenite domain with a well-defined contact zone between them (Figure 3A). Dunite veins in harzburgite are widely interpreted as intensive reactions between olivine-saturated melts and orthopyroxene-rich wall rocks (e.g., Kelemen et al., 1995). This suggests that the formation of the dunite vein may have followed the formation of the olivine-orthopyroxenite domain. Since intensive melting of the pyroxenes during metasomatic reactions requires high melt-rock ratios, dunite veins could represent possible melt migration pathways (Yoshikawa et al., 2018) as it was verified by experiments as well (Morgan and Liang, 2005). This interpretation can be applied to the dunite domain of the xenolith MSZK1301D because of the coarse granular texture dominated by stress-free olivines (Figure 3A). Note, however, that recrystallization after the magmatic event can result in similar petrographic features and also explains the lack of reaction textures in other Group I xenoliths.

Orthopyroxene-enrichment via metasomatism requires silica-rich, orthopyroxene-saturated melts. Such melts appear most commonly in volcanic arc environments. Beneath volcanic arcs, the liberation of fluids from a subducting slab can trigger partial melting in the mantle wedge overlying the slab (e.g., Iwamori, 1998; Grove et al., 2006). These melts formed at high pressure saturate first in olivine only, leading to the dissolution of pyroxenes during melt-rock reactions and the formation of dunites (Kelemen et al., 1992). However, ongoing melt-rock reactions can lead to gradual modification of the melt from basalt to basaltic andesite composition, resulting in orthopyroxene saturation (Morgan and Liang, 2005) at lower pressure in constant temperature circumstances (Kelemen et al., 1992). This process can explain the co-precipitation of orthopyroxene and olivine with simultaneous clinopyroxene dissolution; in other words, the formation of harzburgites (Kelemen et al., 1998), even when the magma mass is decreasing (Kelemen et al., 1992).

The migration of silica-rich (even quartz-oversaturated) melts is not unprecedented in the BBHVF. Bali et al. (2008) described a highly unusual quartz-bearing orthopyroxene-rich websterite xenolith from Szigliget (∼10 km southwest of Mindszentkálla; Figure 1B). The quartz appeared as rounded or needle-shaped crystal inclusions in orthopyroxenes and was interpreted to have co-crystallized with its host. The quartz-bearing xenolith was also rich in silicate melt inclusions (SMIs) hosted by pyroxenes. The metasomatic melt represented by the primary SMIs was found to be rich in SiO₂ (66–72 wt%) -and alkalis, and poor in MgO, FeO and CaO. Furthermore, it was strongly enriched in LREE and large ion lithophile elements (LILE), and displayed negative Nb, Ta and Sr and slightly positive Pb anomalies. These geochemical characteristics led Bali et al. (2008) to conclude that a quartz-saturated silicate melt was released from metasedimentary rocks of a subducted oceanic slab. Other studies also suggested migrating subduction-related melts beneath the BBHVF (e.g., Downes et al., 1992; Bali et al., 2002, 2007), occasionally evidenced by silica-rich (>55 wt%) glass veins (Créon et al., 2017b).

We suggest that beneath MSZK, subduction-related extensive metasomatism was caused by a silica-rich, orthopyroxene-saturated melt. In addition to the metasomatism-induced orthopyroxene enrichment, several lines of geochemical evidence point to the supra-subduction zone (SSZ) origin of the MSZK harzburgites. The most important geochemical characteristics are the high Mg# (>0.90; Figures 4B,C) and low Al₂O₃ in orthopyroxenes (mostly <2.5 wt%; Supplementary Figure S3C), as well as the Cr₂O₃-rich spinels (>29 wt%; Figure 8A). The presence of clinopyroxenes displaying U-shaped REY patterns (Figures 6D,E) is also in accordance with SSZ environment (Parkinson and Pearce, 1998). However, the LREE and HREE enrichment compared to Eu, Gd, Tb depletion is not as pronounced (Figures 6D,E) as in other BBHVF xenoliths (Bali et al., 2007; Dobosi et al., 2010; Berkesi et al., 2012; Embey-Isztin et al., 2014). The only Group I xenoliths showing no U-shaped REY patterns are MSZK1301A, MSZK1301B and MSZK1311 samples. The latter two xenoliths are also located outside the field constrained by SSZ xenoliths on the Al₂O₃-CaO diagram of orthopyroxenes (Figure 8B). The reason for this geochemical deviation from the rest of the Group I xenoliths could be due to the influence of a basaltic melt (probably the same as the one responsible for the formation of Group II xenoliths; see chapter 6.3), based on enrichment of these two xenoliths in some basaltic elements (Ti, Al, Fe in clinopyroxene; Ti, Al and Ca in orthopyroxene; Al in spinel). Note that one single-lithology harzburgite xenolith (MSZK1309) in Group I shows a positive Eu anomaly (Figure 6E). Such a positive anomaly is commonly explained by the dissolution of plagioclase in the subducting oceanic crust (Dai and Zheng, 2019), which is in agreement with the proposed evolution of that xenolith.

Group II xenoliths consist of variable lithologies (harzburgite, orthopyroxenite, websterite) of two composite xenoliths (MSZK1305; MSZK1306B) (Table 1). Generally, these lithologies have higher Ti, Mn, Ca and lower Si, Cr, Ni, and Mg# compared to those of Group I (Figure 4; Supplementary Figure S3). The enrichment in basaltic elements is weak in the harzburgite, moderate in the orthopyroxenite and significant in the websterite domains (Figure 4; Supplementary Figure S3). These characteristics, along with the convex-upward REY patterns of clinopyroxenes (Figure 6F), point to an interaction with mafic melts. However, the lack of melt infiltrations, reaction coronas and compositional zoning at the margins of the MSZK xenoliths all suggest that the host basalt is likely not the metasomatic agent. This is further supported by the fact that xenoliths were collected from pyroclastic sequences (Supplementary Figure S1B), which can preserve the original physical and chemical properties of the xenoliths better, due to their fast cooling on the surface.

A general geochemical feature in all lithologies of Group II xenoliths is the negative Eu anomaly (Figure 6F). It may either be related to the crystallization of plagioclase during magma evolution, or to redox state changes affecting the distribution of Eu compared to other REE (Griffin and O’Reilly, 2007). The alkali basaltic melts in the BBHVF originate from the asthenosphere far deeper than the stability field of plagioclase (e.g., Embey-Isztin et al., 1993), thus residual plagioclase retaining Eu in the source region is out of option. In the BBHVF, only one upper mantle xenolith is reported containing plagioclase (Embey-Isztin, 1976), hence the mass accumulation of plagioclase in the lithospheric mantle is also unlikely. Therefore, it is suggested that the negative Eu anomaly results from changing redox conditions of the melt during its upward migration.

The harzburgite domains in Group II composite xenoliths likely represent a mantle volume originally belonging to Group I, subsequently affected by a younger event. This is based on the observation that petrographic features, including modal composition and textural characteristics, do not differ in the single-lithology harzburgites and the harzburgite domains of composite xenoliths (Figure 3; Supplementary Figure S4; Table 1), implying a common origin for them. However, they have several different geochemical features as follows. The harzburgite domains of the composite xenoliths are out of the OSMA (olivine-spinel mantle array; Arai, 1994) range (<87 forsterite content in olivine along with 0.47–0.48 Cr# in spinel; Figure 4A). Furthermore, they are among the harzburgites with the lowest Mg# in their pyroxenes (∼0.88 and ∼0.91 Mg# for orthopyroxene and clinopyroxene, respectively; Figures 4B–E). The subordinate clinopyroxenes in harzburgite xenoliths have lower REY contents (6–24 ppm ΣREY; Figures 6D,E; Supplementary Table S9) compared to the harzburgite domains of composite xenoliths (∼148 ppm ΣREY; Figure 6F; Supplementary Table S9). The other difference is the REY pattern showing flat and enriched distributions in the former (Figures 6D,E) and convex-upward distribution in the latter group (Figure 6F). All these geochemical features suggest that the single-lithology harzburgites and the harzburgite domains of composite xenoliths have different evolution.

The convex-upward pattern of clinopyroxenes in harzburgite domains of the composite xenoliths is also characteristic for clinopyroxenes of the orthopyroxenite and websterite domains in composite xenoliths (Figure 6F). This similarity suggests that the younger geochemical event, presumably linked to mafic melts, that likely added the orthopyroxenite and websterite lithologies to the mantle volume, overprinted the trace element characteristics of the neighboring parts of the presumably older harzburgites via chromatographic processes (Navon and Stolper, 1987), leading to cryptic metasomatism (Dawson, 1984). In other words, the harzburgite host rock of the composite xenoliths likely failed to preserve its original trace element attributes during the formation of orthopyroxenite and websterite lithologies.

The orthopyroxenite lithology, albeit in a subordinate amount, is present in various settings below the Moho, including ophiolitic mantle (e.g., Le Roux et al., 2014) orogenic peridotites (e.g., Malaspina et al., 2006) and upper mantle xenoliths (e.g., Chen and Zhou, 2005; Wang et al., 2008). Most of the studies explain the formation of orthopyroxenites with melt/fluid and peridotite interactions based on resorbed olivines. Several studies suggest that siliceous melts and fluids derived from an oceanic slab could be the most probable metasomatic agents responsible for orthopyroxenite formation (Ertan and Leeman, 1996; Chen and Zhou, 2005). Indeed, orthopyroxene-rich reactive boundary layers were described to form between basaltic andesite melt and peridotite at high pressure (≥2.5 GPa) and temperature (≥1,375 °C) (Mallik and Dasgupta, 2012). However, the crystallization of orthopyroxene-oversaturated reacting melt prevents further melt-rock reactions (Wang et al., 2016). Alternatively, reactions between hydrous ferro-basalt and peridotite in upper mantle conditions (0.8–2 GPa; 1,250–1,385°C) can also lead to the formation of orthopyroxenites as suggested by experimental runs of Wang et al. (2016). Based on their model, the diffusion of water from the hydrous melt into the host peridotite can induce melting in the peridotite, resulting in olivine-rich residues. The mixing of the migrating melt with the melt produced in the host rock results in a third, hybrid melt. This hybrid melt is responsible for the olivine dissolution and orthopyroxene precipitation along the melt-rock interface. Other studies on natural samples also emphasize the possible role of basaltic melt (Arai et al., 2006) and water (Smith et al., 1999; Tamura and Arai, 2006) in the formation of orthopyroxenites. The significance of water is further supported by the common appearance of hydrous phases such as amphibole or phlogopite in orthopyroxenites (Ertan and Leeman, 1996; Chen and Zhou, 2005). The reacted melts assumed to contribute to orthopyroxenite formation have similar water contents in experimental runs (0.61–0.66 wt%; Wang et al., 2016) and natural samples (0.57–1.24 wt%; Arai et al., 2006). Note that orthopyroxenite lithology was only described from the Little Hungarian Plain Volcanic Field (LHPVF; Figure 1A) in the CPR so far (Szabó et al., 1995).

In the MSZK composite xenoliths, the growth of orthopyroxene of the orthopyroxenite domain at the expense of olivine of the neighboring harzburgite domain (Figures 3B,C) suggest that the orthopyroxenite domain was formed later than the harzburgite host rock. It also means that the websterite domains and the amphibole-phlogopite vein, which are in contact with the orthopyroxenite domains only (Figures 3B,C), must have formed after the harzburgite parts.

As it was indicated before, we also suggest that mafic melts are responsible for the orthopyroxenite formation beneath the MSZK locality, based on their basaltic element enrichment (Figure 4; Supplementary Figure S3) and convex-upward REY patterns of clinopyroxenes (Figure 6F). Analyses concluded high water contents for the BBHVF basalts (2.0–2.5 wt%; Kovács et al., 2020), which is capable promoting orthopyroxenite vein formation in the mantle according to the experiments of Wang et al. (2016). Furthermore, the upward-moving basaltic melts likely contacted harzburgite volumes beneath the MSZK locality, which were produced by orthopyroxene enrichment during a previous metasomatism by subduction-related melts and fluids. The dissolution of the orthopyroxenes resulted in more evolved melt compositions, and possible orthopyroxene-saturation and even oversaturation after interacting with large volumes of harzburgite. These processes resulted in the orthopyroxenite formation. The role of harzburgite lithology, as a source of melts responsible for the formation of orthopyroxenite veins, was also highlighted by Bénard et al. (2022). Without orthopyroxene dissolution, basaltic melts would remain orthopyroxene undersaturated with their movement restricted to the olivine-clinopyroxene cotectic line on schematic olivine-clinopyroxene-SiO₂ phase diagram (Kelemen, 1990).

We estimated the equilibrium melt composition for orthopyroxenites as well, using clinopyroxene-basaltic melt and orthopyroxene-basaltic melt distribution coefficients (Ionov et al., 2002; Adam and Green, 2006) (Figure 10A). The estimated equilibrium melt compositions (Supplementary Table S12) show the best fit with SMI reported by Szabó et al. (2009) (Figure 10A). These SMI are hosted by amphibole-bearing spinel lherzolite xenoliths from Szigliget (Figure 1B). The trapped SMI likely represent a melt with a trace element content similar to that of the host alkali basalts, but more enriched in silica (Szabó et al., 2009). According to Szabó et al. (2009), such melt may have formed by the interaction between an alkali mafic melt rising from the asthenosphere and a mantle portion previously metasomatized by a slab-derived melt. This evolution history is very similar to what we suggest for the melt responsible for the orthopyroxenite formation in the upper mantle beneath the MSZK locality. Note that the trace element composition of several studied SMI shows enrichment in U and Pb (Szabó et al., 2009), a feature also characteristic of the Group IB harzburgites of the MSZK xenoliths (Figure 6B). These properties were likely inherited from the metasomatized mantle wedge harzburgites representing the wall rocks next to the upward migrating mafic melts.

The appearance of websterite followed the orthopyroxenite formation, as indicated by resorbed orthopyroxenes in the latter along the contact zone (Figure 3D). The large (3,000–5,000 µm), tabular, stress-free clinopyroxenes of the websterites, with enclosed orthopyroxenes, amphiboles and apatites (Figure 3D), are petrographic features typical of cumulates representing crystallized products (e.g., Perinelli et al., 2017). The enrichment of basaltic elements in the rock-forming minerals of the websterites is even higher compared to the orthopyroxenites (Figure 4; Supplementary Figure S3). Together with the convex-upward REY patterns of clinopyroxenes (Figure 6F), this suggests that websterite formation is related to mafic melts. To further test this hypothesis, the equilibrium melt of websterites was estimated with the use of distribution coefficients for clinopyroxene-basaltic melt and orthopyroxene-basaltic melt (Ionov et al., 2002; Adam and Green, 2006) (Figure 10B). The resulting equilibrium melt composition shows the best overlap with basaltic melts defined by the SMIs in the study of Szabó et al. (2009) (Figure 10B), similarly to what was found for the MSZK orthopyroxenites described above (Figure 10A).

Clinopyroxene is often the liquidus phase in basalts beneath the crust, leading to the formation of clinopyroxene-rich lithologies upon crystallization (e.g., Downes et al., 2004). Thus, the dominance of clinopyroxene in the websterites (∼70 vol%) suggests that they could represent crystallized melt migration paths. The significant modal abundance of orthopyroxenes (∼30 vol%) in the MSZK websterites can be explained by the incomplete alteration of orthopyroxenites or by orthopyroxene dissolution and subsequent orthopyroxene saturation in the basaltic melt. The higher equilibration temperatures of websterites (T_NG>1,000°C; T_REE>1,100°C) compared to orthopyroxenites (T_NG<1,000°C; T_REE<1,100°C) (Table 1) further indicate that melt crystallization is a more important process in websterite formation than melt-rock reactions, which often result in lower temperatures due to the mixing of melt with the colder wall rocks. All these suggest that the MSZK websterites likely represent mafic melts crystallized along melt channels. This is consistent with other studies interpreting websterite veins as products appearing along recrystallization fronts of pervasive melt migration channels (Garrido and Bodinier, 1999; Dantas et al., 2009).

The amphibole-phlogopite vein crosscutting the orthopyroxenite domain of xenolith MSZK1305 (Figure 3C) is likely a product of the last event in the upper mantle of the study area. It is possibly derived from the mafic melt that went through high-degree fractional crystallization and became highly enriched in incompatible elements (Figure 7). Amphibole-rich veins with similar proposed origin were found in composite xenoliths from Szigliget (Embey-Isztin, 1976; Bali et al., 2018) (Figure 1B).