Tephrostratigraphy and Magma Evolution Based on Combined Zircon Trace Element and U-Pb Age Data: Fingerprinting Miocene Silicic Pyroclastic Rocks in the Pannonian Basin

1 Introduction

Eruptions of large volume silicic magma are usually accompanied with caldera collapse, pyroclastic flows (ignimbrites) traveling up to 100 km distance and pyroclastic falls dispersed over 100’s of kilometre away from the source (e.g., de Silva, 1989; de Silva and Francis, 1989; Wilson et al., 1995; Wörner et al., 2000; Ellis et al., 2012; Lipman and Bachmann, 2015; Barker et al., 2020). Therefore, their volcanic products act as important marker beds in stratigraphy, in addition to providing key information about the origin of the magmas. In the past, several geochemical tools, such as composition of glass and the main mineral phases have been used to correlate scattered deposits of an eruption event (de Silva and Francis, 1989; Shane and Smith, 2000; Pearce et al., 2002; Harangi et al., 2005; Pearce et al., 2007; Lowe, 2011; Aydar et al., 2012; Gisbert and Gimeno, 2016; Lowe et al., 2017; Petrelli et al., 2017). However, syn- and post-eruption alteration often obscures the original compositional features of pyroclastic deposits. A possible solution to correlate scattered, variable altered silicic volcanic deposit is to use zircon chemistry, as zircon is both a common accessory mineral in silicic magmas and it is very resistant to mechanical and chemical erosion. Since U and Th can effectively enter its crystal lattice, it is widely used in geochronology to date eruption events as well as to constrain the timescale of magma storage evolution. Zircon contains significant amount of other trace elements, such as Y, Ti, middle and heavy rare earth elements and their concentrations depend on the nature of the magma as well as the crystallization conditions, which can then be used for tectonomagmatic provenance studies (Hoskin and Ireland, 2000; Belousova et al., 2002; Hoskin and Schaltegger, 2003; Grimes et al., 2015). Therefore, this mineral can yield a time-stamped record of the magma types and it can be a key-mineral in correlating eruption products (Aydar et al., 2012) even in case of strong alteration.

The Pannonian Basin (eastern-central Europe; Figure 1) produced among the largest silicic volcanic eruptions in Europe for the last 20 Myr (Lukács et al., 2018a). Volcanic ash material was deposited throughout central and southern Europe (e.g., Rocholl et al., 2018; Lukács et al., 2018a; Rybár et al., 2019; Sant et al., 2020) and therefore plays an important role in constraining the regional Paratethys chronostratigraphy in this area (Steininger et al., 1988; Sant et al., 2017). Based on a thorough zircon-based study, Lukács et al. (2018a) could identify at least five major eruption events from 18.2 Ma to 14.4 Ma mostly sourced around the southern foreland of the Bükk Mts., called Bükkalja Volcanic Field (BVF; Northern Hungary, Figure 1). Miocene silicic pyroclastic deposits are widespread in the Pannonian Basin although most of them are buried by young sediments due to post-volcanic subsidence. Nevertheless, there are many scattered occurrences in Northern Hungary, tens of kilometre from the BVF. In this study, we explore how combined zircon U-Pb dates and trace element content allow fingerprinting the main eruption units of the BVF, and whether such chemical fingerprints of zircons can be used to correlate distal deposits lacking precise age data and/or stratigraphic position. We argue that such a zircon-based correlation technique could be extended to volcanic ash deposits found in different Paratethys basins of central and southern Europe, several hundreds of kilometre from the BVF where volcanic glasses are altered and they yield no more information about the erupted magmas (Lukács et al., 2018a; Rocholl et al., 2018). Zircons are also known to be a useful monitor of magma evolution (Miller et al., 2007; Claiborne et al., 2010; Wotzlaw et al., 2013; Deering et al., 2016; Buret et al., 2017; Szymanowski et al., 2017; Claiborne et al., 2018; Yan et al., 2018; Yan et al., 2020) and can be used to unravel the pre-eruption processes in the crustal reservoirs emplaced in continental crust. Obtaining reliable mineral-melt distribution coefficients is a key for it as well as to infer the composition of the melt from which those zircons crystallized, particularly when no direct information for the melt phase (e.g., strongly altered glass) is available. Thus, our aim with this study is to present a novel approach to use zircons as a correlation tool and magma monitor for a large volume silicic volcanic system.

FIGURE 1. Occurrences of Miocene silicic pyroclastic rocks in Northern Hungary. Studied samples in the Bükkalja volcanic field presented this paper are denoted by numbers (Supplementary Table S1; after Lukács et al., 2018a), whereas other samples used for correlation are shown with sample names (IT = Ipolytarnóc; MSZ = Mátraszele; NEMTI = Nemti; SZVG = Szentkút; Tar = Tar; Tar-Gh = Tar Gömör hill; SZIV = Szilvásvárad). Europe is an ESA satellite picture from 2009.

2 Geological Setting

The Pannonian Basin was formed by major lithospheric extension in the hinterland of the Carpathian-Dinarides orogens during the Miocene (Royden et al., 1982; Tari et al., 1999; Horváth et al., 2006; Horváth et al., 2015). The continental rifting was diachronous, started probably already at around 28 Ma at the southwestern margin, although the major stage occurred during the early Miocene (19–16 Ma) in the western areas and it took place only at about 10–9 Ma in the central-eastern subbasins (Fodor et al., 1999; Tari et al., 1999; Ustaszewski et al., 2010; Toljić et al., 2013; Horváth et al., 2015; Balázs et al., 2016).

The continental extension was accompanied with widespread volcanism involving eruption of various magmas since the early Miocene (Szabó et al., 1992; Harangi, 2001; Seghedi et al., 2004, Seghedi et al., 2005; Harangi and Lenkey, 2007; Seghedi and Downes, 2011; Lukács et al., 2018a; Harangi and Lukács, 2019). The last eruptions occurred in the Late Pleistocene at the northern part of the Pannonian Basin (100 ka; Šimon and Halouzka, 1996; Šimon and Maglay, 2005) and in the southeastern Carpathians (30 ka; Harangi et al., 2015; Molnár et al., 2019). This long-lasting volcanism began with an ignimbrite flare-up event when more than 4000 km³ silicic magma erupted during about 3 Myr (Lukács et al., 2018a). Lukács et al. (2018a) distinguished at least five major eruptions in the BVF (17.5 ± 0.3 Ma Eger ignimbrite, 17.055 ± 0.024 Ma Mangó ignimbrite, 16.816 ± 0.059 Ma Bogács unit, 14.880 ± 0.014 Ma Demjén ignimbrite and 14.358 ± 0.015 Ma Harsány ignimbrite based on pooling of the youngest zircon population in each unit). These eruptions provided significant volume of erupted material covering extensive areas across the Pannonian Basin and elsewhere in Europe (e.g., Lukács et al., 2018a; Rocholl et al., 2018; Rybár et al., 2019; Sant et al., 2020). These eruption products act as important marker layers in the early to mid-Miocene and could help to constrain the regional Paratethys chronostratigraphy used in central Europe.

The volcanic products of three of the main eruption units (Eger, Mangó and Harsány) are rhyolitic and contain crystal-poor pumices (Lukács et al., 2007; Harangi and Lenkey, 2007; Lukács et al., 2015; Lukács et al., 2018a). Phenocrysts are quartz, plagioclase and biotite with accessory minerals of zircon, allanite, apatite and ilmenite. The Demjén ignimbrite is rhyodacitic and a notable feature is the occurrence of amphibole and sporadically clinopyroxene, while quartz is much restricted compared to the rhyolitic units. The most complex eruptive product is the Bogács unit (Czuppon et al., 2012), where a gradual compositional and petrological change is observed upwards. In contrast to the other eruption units, the pyroclastic rocks of the Bogács unit are characteristically crystal-rich (15–40 vol%). Its lower subunit is rhyolitic and contains unwelded and welded ignimbrites (with two types of fiamme), whereas the upper subunit is a mixture of various juvenile clasts (scoria, pumice and hybrid clasts) suggesting pre-eruption mixing and mingling of crystal mushes and melts with dacitic to rhyolitic compositions (Czuppon et al., 2012). Remarkably, the mineral phases and the glasses are usually free of syn-eruption or post-eruption alteration in all units (Harangi et al., 2005) and therefore, they preserve the original compositional features.

This study focuses on finding correlation tools using zircon U-Pb dates and trace element contents to constrain the relationship of scattered silicic pyroclastic occurrences in the BVF and in the northern Pannonian Basin (Figure 1). These pumiceous volcanoclastics were dated previously by K-Ar technique that yielded ages from 20 to 13 Ma (Hámor et al., 1980; Márton and Pécskay, 1998; Póka et al., 2004; Zelenka et al., 2004; Pécskay et al., 2006). They were divided into formal lithostratigraphic units: Gyulakeszi Rhyolite Tuff, 19.6 ± 1.4 Ma; Tar Dacite Tuff, 16.4 ± 0.8 Ma and Galgavölgy Rhyolite Tuff, Harsány Rhyolite Tuff, Cserehát Rhyolite Tuff, Felnémet Rhyolite Tuff with overlapping ages from 14.6 to 9 Ma (Hámor et al., 1980; Márton and Pécskay, 1998; Gyalog and Budai, 2004; Less et al., 2005).

3 Samples and Analytical Methods

More than 950 zircon crystals were analyzed for U-Pb dating in combination for trace element content from 30 samples. Samples used in this study are presented in Figure 1, Table 1 and Supplementary Table S1 with their stratigraphic position, GPS coordinates, lithology and volcanic facies. Geological background of the studied localities is given in the Supplementary Data Sheet S1. The BVF samples represent all the units identified by Lukács et al. (2018a); Supplementary Table S1), while samples out of BVF are from several scattered localities in the northern part of the Pannonian Basin (Figure 1; Supplementary Table S1).

TABLE 1. Localities and main features of the studied samples in Northern Hungary.

Zircon crystals of BVF units presented in Lukács et al. (2015), Lukács et al. (2018a) were targeted for trace element analysis. We used the same sample aliquots and crystal mounts as Lukács et al. (2015), Lukács et al. (2018a) measured for zircon U-Pb dating (see Supplementary Table S1; Lukács et al., 2018b). In addition, we performed combined U-Pb dating and trace element analysis for the samples collected in Northern Hungary using the same methodologies. The list of the studied samples with GPS coordinates of the localities is presented in Table 1.

Zircon U-Pb isotopes and trace element contents were analyzed by LA-ICP-MS at the ETH Zürich. The laser ablation system was a Resonetics Resolution 155 type which was coupled to a Thermo Element XR SF-ICP-MS. We used 30 µm spot size, 5 Hz repetition rate, 2.0 J/cm² energy density (fluence) and 20 s (in session for only geochemistry) or 40 s (in sessions for U-Pb dating and chemistry) ablation time after five cleaning pulses and 30 s of gas blank acquisition.

U-Pb dating of zircons was preferably taken on rim parts of the sectioned crystals and for each sample we conducted at least 45 spot analyses on various zircon grains. We performed the analyses in three sessions. Analytical conditions were similar to what Lukács et al. (2018b) used and summarized in Supplementary Data Sheet S1. GJ-1 reference zircon (Jackson et al., 2004) was used as a primary standard, while Zircon 91500, Plešovice, Temora2, OD-3, AUSZ7-1, AUSZ7-5, AUSZ8-10 were measured as validation reference materials (Wiedenbeck et al., 1995; Black et al., 2004; Sláma et al., 2008; Iwano et al., 2013; Kennedy et al., 2014; von Quadt et al., 2016, respectively). The average precision of the reference materials (RM) ranged from 0.8% to 2.4% (2 SE) in case of secondary RM older than 30 Ma, while the younger have 6.4–10.1% (2 SE) average precisions (further details are in Supplementary Table S2). Validation reference materials, covering an age interval between 1063.5 Ma and 1.431 Ma, were used to correct the matrix-dependent age offsets following the procedure described in Sliwinski et al. (2017). The youngest reference material was AUSZ8-10 of which six crystals were measured by ID-TIMS method to obtain its reference age value. Analytical conditions were the same as in von Quadt et al. (2016). Results of the analyses are presented in Supplementary Table S2, the weighted mean age of the six crystals gives 1.431 ± 0.002 Ma (MSWD = 0.63). For the LA-ICP-MS data reduction, we used IOLITE (Paton et al., 2011) combined with VizualAge (Petrus and Kamber, 2012) software. The Th disequilibrium correction was performed after alpha dose-correction assuming a constant Th/U partition coefficient ratio of 0.33 ± 0.063 (1 σ; Rubatto and Hermann, 2007) and using the equations of Schärer (1984). Zircon U-Pb dates for the 10 samples in Northern Hungary were not corrected for common Pb contents; however, during data reduction, integration intervals were set to exclude the common Pb contaminated signal intervals and data were filtered according to their discordance ([(²⁰⁷Pb/²³⁵U Age)–(²⁰⁶Pb/²³⁸U Age)]/(²⁰⁷Pb/²³⁵U Age) < 10%). Average uncertainty of the individual zircon dates is given as 2 SE and it is between 1.3 and 2.1% (relative 2 SE).

In case of trace element analysis we used NIST610 as primary reference material and zircon 91500 for quality control, and an in-house synthetic reference material (Synthetic Zircon Blank) to correct the Nb concentrations in some sessions. We analyzed the samples in 12 sessions and in case of 9 sessions U-Pb isotopes and trace elements were determined simultaneously (see Supplementary Table S1). With some variations between sessions the target elements were Si, Zr, REE, Y, Hf, P, Nb, Ta, U, Th, Ti (see Supplementary Table S1), and either Al, Rb, Ba, Ca, Fe or Al, Sr were measured for monitoring glass, apatite and iron oxide inclusions. Si (15.2 wt% in Zircon) was used as internal standard for data reduction done by IOLITE. Spot compositions that has larger values than the detection limit for Al, Fe, Ca were discarded from the database. Further data filtering was made using the results of Principal Component Analysis (PCA) and data having anomalously high values in Y, La, HREE, P, Nb, possibly due to incorporation of glass or mineral inclusion in the analyzed compound volume were also omitted. Statistical analyses were performed by the NCSS 2020 Software (NCSS, LLC. Kaysville, Utah, United States, ncss.com/software/ncss).

4 Results

4.1 Zircon Geochemistry of the Bükkalja Volcanic Field Samples

The analyzed zircons (data are presented in Supplementary Table S1) have Hf contents between 7100 and 12,000 ppm, the Zr/Hf ratio is in the range between 38 and 66 (Figure 2 and Supplementary Figure S1 in Supplementary Data Sheet S2). The Th/U ratio varies from 0.2 to 1.1 and U contents from 73 to 2300 ppm, although majority of the U contents is between 73 and 1000 ppm. REE and Y concentrations (330–4300 ppm) positively correlate with U and Th contents (Figure 2). The U/Yb ratio varies between 0.4–2.7 which is typical for continental crustal igneous rocks (Grimes et al., 2007). Zircon Ti concentrations cluster between 2–8 ppm for all units except for the samples of the Bogács unit, where zircon crystals show higher concentration values (4–13.5 ppm). Below, we describe briefly the principal trace element characteristics of the five main eruptive units of the BVF.

FIGURE 2. Zircon trace element variation diagrams for the main BVF eruption units (left panel) and for the scattered distal samples in Northern Hungary (right panel).

The Eger ignimbrite is represented by zircon crystals from one single pumice clast collected at the lower level of the Tihamér quarry. The Mangó ignimbrite has a much wider identified distribution in the BVF. Here, to represent it, we analyzed zircons from 6 different localities including the upper level of Tihamér quarry (Supplementary Table S1). Zircons of Eger ignimbrite and Mangó ignimbrite show roughly similar chemical variation in Nb (<8 ppm), P (145–680 ppm), Ti (2.6–8.5 ppm) and Eu anomaly (0.15–0.37). However, the Mangó ignimbrite has a little bit wider range in zircon Hf (8600–12000 ppm) and Y (540–2600 ppm) concentrations and Th/U (0.2–1.1) and Zr/Hf (45–54) ratios compared to the zircons from the Eger ignimbrite (Hf = 9200–10900 ppm; Y = 386–2470 ppm; Th/U = 0.4–0.9; Zr/Hf = 47–52). Zircons of the Eger and Mangó ignimbrites can be well separated in the Yb/Dy vs. Th/U plot (Figure 2), where the Eger zircons have typically higher Yb/Dy ratio at the same Th/U values.

Zircon crystals of Bogács Unit were analyzed from four types of juvenile clasts (Supplementary Table S1): a large fiamme (A-fiamme type; Td-Fi), a large hybrid (Td-H2N) and a small hybrid (Td-Hk1) clasts and a scoria (Td-S) of Tibolddaróc locality (Czuppon et al., 2012; Lukács et al., 2018a). They all show fairly similar zircon compositions and have distinct chemical features compared to the other BVF zircons. The difference is obvious in the restricted and lower Hf concentrations (generally between 7100 and 9500 ppm) together with higher Zr/Hf ratios (47–66) and also in the lower Yb/Dy ratios (2.4–4.3). These zircon crystals have Th/U ratio between 0.4–1.0, Nb, Y and P contents of 1–9 ppm, 660–3900 ppm and 590–1540 ppm, respectively. Ti content (4–13.5 ppm) of Bogács zircons shows positive correlation with the Eu-anomaly. The small hybrid clast (Td-Hk1) has zircons that have noticeably restricted Y, Nb, Ta, U, P and Hf concentrations within the Bogács zircon population.

Trace element composition of zircon of the Harsány and Demjén ignimbrite units was published in Lukács et al. (2015). Zircons of Demjén ignimbrite show similar Hf (8000–11200 ppm) range as Mangó and Eger ignimbrite, but a usually lower and restricted variation in Nb (mostly <5 ppm), Y (mostly 300–2000 ppm), P (mostly 100–400 ppm) contents. The Th/U and Zr/Hf ratios range between 0.35 and 0.8 and between 44 and 53, respectively. They show a wide range in the Yb/Dy ratio from 3.3–8.0, and give the highest values within the BVF zircons. These high Yb/Dy ratios at a given Th/U values overlap the zircons from the Eger ignimbrite. Noteworthy, a bimodality in zircon composition can be observed, particularly in the Eu/Eu* and the Yb/Dy values and this may imply involvement of two distinct silicic melts and crystal cargo in the eruptive product. On the other hand, zircons from the Harsány ignimbrite give wide compositional variation in most elements and trace element ratios. Hf concentrations span from 8300 to 12000 ppm, while Th/U from 0.2 to 0.85 and Zr/Hf ranges between 38 and 54. The wide variation is typical in case of Y (450–4300 ppm), Nb (up to 17.5 ppm) and P (150–780 ppm) contents. Remarkably, zircons have typically the highest negative Eu-anomaly (Eu/Eu* is typically 0.07–0.24) among the BVF samples. They have also low Yb/Dy values along with the low Th/U ratios (Figure 2). Zircons from all the other samples in the BVF share compositional features of one of the main eruptive units described above.

4.2 U-Pb Zircon Chronology and Geochemistry of Pyroclastic Samples From Localities Outside the Bükkalja Volcanic Field

Majority of the single crystal dates varies between ∼20 and ∼14 Ma, and only a few extremely old xenocrystic cores were measured in some samples (∼610–350 Ma). Each sample has multicomponent age populations according to the calculated MSWD (Mean Squared Weighted Deviation) values of weighted mean ages, which have therefore no geologic meaning. Thus, we calculated the youngest population age (low-MSWD weighted mean) without the younger outlier dates (masked by Pb loss) for each sample along with the propagated uncertainty including the 1.5% external error of LA-ICP-MS measurements. These youngest zircon crystallization ages are interpreted as representing the closest ages to the volcanic eruptions. The results are summarized in Table 2, and are presented in Supplementary Figure S2 in Supplementary Data Sheet S2.

TABLE 2. Summary of zircon U-Pb dates and determined ages of Miocene ignimbrite samples from distal occurrences in Northern Hungary.

The youngest age populations of the studied zircon samples are in the range from 14.39–14.25 Ma (SZIV and SZVG) to 17.49 Ma (IT-PUB), which exactly overlaps the age range shown by the BVF main eruption units. The two ignimbrite samples from Tar yielded the same zircon ages (14.9 ± 0.2 Ma and 15.0 ± 0.2 Ma), overlapping with those of the Demjén ignimbrite (14.88 Ma; Lukács et al., 2015). The SZVG and SZIV samples yielded eruption ages similar to the Harsány ignimbrite within uncertainty (14.4 Ma; Lukács et al., 2015). The calculated eruption ages of the Ipolytarnóc and Nemti localities overlap those the Eger and Mangó ignimbrites (17.5 and 17.1 Ma; respectively; Lukács et al., 2018a), whereas the bentonite sample from MSZ has younger eruption age akin to the Bogács unit (16.7 Ma; Lukács et al., 2018a).

Trace element composition of zircons from pyroclastic deposits in the Northern Hungary was analyzed along with the U-Pb isotopes (full data are presented in Supplementary Tables S1, S2). The compositional features overlap those of the BVF main units (Figure 2) suggesting that they belong to the same volcanism.

5 Discussion

5.1 Correlation of Eruption Units Based on Zircon Fingerprints

5.1.1 Zircon Trace Element Discrimination

Trace element composition of zircons from different eruption units shows a wide variety, although certain trace elements and trace element ratios have distinct values characteristic of the eruption units (Figure 2). The within-unit zircon chemical variation is partly due to the fact that LA-ICP-MS spot analysis takes a compound volume of several growth zones of the crystals and forms therefore an average composition of the spot. Although we focused on measuring the outer parts of the crystal grains, age-zoned crystals could yield distinct composition. Thus, outlier data occurs, when older domains or part of them are incorporated into the measured volume. In addition, tiny glass or mineral inclusions incorporated into the measured volume can also modify the trace element composition. However, they can be recognized and omitted from the data set. The distinct compositional characteristics of the main eruption units suggest that LA-ICP-MS spot data for zircons can be used effectively for correlation purposes. Within the trace elements and trace element ratios, we found that P, Y, Nb, Nd, Hf, U and the Th/U, Yb/Dy, Eu/Eu* ratios have the most discriminating role.

Multivariate mathematical techniques offer a powerful tool to examine datasets with large number of samples and variables and increase the interpretability with minimizing the information loss. Although they have already widely applied in tephrochronology (e.g., Lowe, 2011 and references therein), there are only a few applications in zircon studies (e.g., Bell and Harrison, 2013; Barth et al., 2017). We show here that such techniques could provide a promising perspective in the interpretation of zircon trace element data set.

In the first step, we conducted multivariate mathematical techniques to characterize the main eruption units defined by Lukács et al. (2018a) for the BVF and then, in the second step, we explored correlation of distal pyroclastic deposits with the main units. We applied Principal Component Analysis (PCA) to describe the variation of the data set and to recognize the outlier data. This was followed by the application of the Multiple Discriminant Analysis (MDA) to find predictive equations to distinguish successfully the main eruptive units by their zircon trace element abundances and trace element ratios. The U/Pb ages were not involved in this analysis; therefore, this result could be used as an independent check how temporally distinct eruptive events can be differentiated. The MDA can find the best solution to separate naturally occurring groups in a data set using discriminant functions, which are linear combinations of the variables. The first discriminant function maximizes the differences between groups, while the second, third and so on functions have subsequently less ability to do that. Discriminant scores of each sample are calculated based on the equations given for the discriminant functions and the score plots give a visual impression of how well the discriminant functions classify the data.

We run MDA using zircon compositions from samples of the five main eruption units (Eger, Mangó, Bogács, Demjén and Harsány) and used canonical score plots to check the separation of the compositional groups. We found that the yielded canonical scores effectively discriminate these units (and Supplementary Figure S3 in Supplementary Data Sheet S2). Based on the canonical coefficients on the discriminant functions (DF), the scores of the samples can be calculated as follows:

The discriminating ability of each variables on the specific discriminant function is expressed by the standard canonical coefficients. For instance, Nb, Eu/Eu*, Th/U and Yb/Dy ratios have the greatest separation ability for the studied zircons as indicated by the highest coefficient values on the discriminant function1, while Y, Nb and U have the greatest importance in the second discriminant function (Table 3). Based on the discriminant score1 and score2 values, the five ignimbrite groups defined by their eruption ages can be readily separated (Figure 3). Zircons from the 17.5 Ma Eger and the 14.9 Ma Demjén units show some overlap, although they can be distinguished based on the canonical score2 and score4 (Supplementary Figure S3 in Supplementary Data Sheet S2 along with the MDA result statistics).

FIGURE 3. Multiple discriminant analysis plot with the two most significant discriminant functions (DF1 and DF2) can effectively distinguish the 5 main eruptive units of the BVF (top left panel) based on combination of trace elements and trace element ratios from zircon. This plot is used to correlate zircons from further pyroclastic deposits in the BVF (top right panel) and distal occurrences in Northern Hungary to one of the main eruption units (two panels at bottom).

TABLE 3. Standard canonical coefficients showing the contribution of each variable to the separation ability of the five main groups (eruptive units).

A particular interest on this analysis is how zircons with broadly similar or close ages (e.g., Eger vs. Mangó units and Demjén vs. Harsány units) can be separated based on trace element characteristics. The canonical variate score1 and score2 values successfully fulfill this requirement suggesting that such plots can be effectively used to correlate unknown zircon samples with the main eruption units of Bükkalja (Figure 3). Zircons of the 17.5 Ma Eger and the 17.1 Ma Mangó ignimbrite units as well as the 14.4 Ma Harsány and the 14.9 Ma Demjén ignimbrite units are distinguished primarily based on the canonical variate score1 values, whereas the zircons of the 16.8 Ma Bogács units are also separated in this plot. Thus, this multivariate mathematical analysis supports the notion that each eruption unit has a specific geochemical fingerprint and this is reflected also in the zircon trace element compositions. Furthermore, 2-group discrimination analysis (2GDA) can be applied to separate more clearly the 17.5 Ma Eger and the 17.1 Ma Mangó ignimbrite units based on the following discriminant function (Figure 4):

Zircons of the two oldest ignimbrite units can be readily distinguished based on the DF, where the Nb concentration and the Th/U and Yb/Dy ratios have the largest contributions. Zircons from the Eger ignimbrite have DF < -2 values, whereas zircons from the Mangó ignimbrite unit have typically DF > -1 values. In summary, both the multi-group and the 2-group discriminant plots provide a fast correlation technique based on zircon trace element composition.

FIGURE 4. Two-group discriminant analysis plot to distinguish the two oldest ignimbrite units: the 17.5 Ma Eger ignimbrite and the 17.1 Ma Mangó ignimbrite can be readily distinguished based on combined trace element content and trace element ratios of zircon (left panel). This plot is used to find correlation of further ignimbrites in the BVF (middle panel) and distal pyroclastic rocks in Northern Hungary (right panel) with similar ages with one of the two oldest eruptive units.

Having the best discriminant functions for the five already identified eruption groups in the BVF, the discriminant score plots can be used to find the appropriate eruption units for deposits, which cannot be unambiguously classified because of lack of proper knowledge on the position in the stratigraphic column or the large uncertainties in geochronological data including the zircon U-Pb ages analyzed by LA-ICP-MS (i.e., ±200–300 ky). For instance, the well-exposed Szomolya ignimbrite (Lukács et al., 2005) is one of the oldest eruption products in the Bükkalja volcanic field, but it is not clear whether it belongs to the Eger or to the Mangó ignimbrite unit. Zircon trace element signature in the MDA and 2GDA plots suggests an Eger-affinity (Figures 3, 4). At Tibolddaróc, there is an almost continuous section with several volcanic eruption beds between the Bogács unit and the Harsány ignimbrite (Lukács et al., 2007, Lukács et al., 2015; Biró et al., 2020). The accretionary lapilli-bearing layer (sample Td-L in Lukács et al., 2007; Lukács et al., 2015) just above the Bogács unit has zircons, which show a clear affinity to the Bogács unit (Figure 3), i.e., this fall unit was part of the Bogács eruption phase. On the other hand, zircons from the Td-E lapilli tuff have trace element composition akin to the slightly younger Harsány ignimbrite (Figure 3), suggesting that this belongs to a precursor eruption before the large ignimbrite-forming volcanic event at 14.4 Ma. The 16.2 Ma Td-J zircons differ from the Bogács unit and has an own, distinct geochemical signature as does the oldest, 18.2 Ma Csv-2 lapilli tuff (Figure 3), i.e. they represent additional volcanic event apart from the five main ones. Although the number of analyzed zircon grains from these units is rather small, it is important to recognize that trace element fingerprint of zircons from the Csv-2 sample differs from that of the Eger ignimbrite unit, which is the oldest volcanic product on the surface within the BVF.

Correlation of ignimbrite occurrences in the surroundings of the BVF was also performed using the multiple and the 2-group discriminant analysis techniques described above. The results for representative cases are shown in Figure 4. The first-order conclusion is that all except for one of the distal ignimbrite occurrences can be correlated with one of the five main eruption units (Lukács et al., 2018a) identified in the BVF. The Ipolytarnóc ignimbrite has a particular importance because it covers a famous early Miocene fossil footprint sandstone and this site is awarded by a European Diploma as a protected area with outstanding geological feature. It was formerly suggested that this volcanic unit could belong to one of the oldest Miocene pyroclastic horizons in the Pannonian Basin, but the age has remained controversial (Márton et al., 2007a). Taking only the zircon composition data, the discriminant scores show a clear affinity to the 17.5 Ma Eger ignimbrite unit, the oldest volcanic unit on the surface in the BVF. In contrast, ignimbrites from the nearby area (e.g., Nemti, Mátraszele) can be correlated rather with the slightly younger, 17.1 Ma Mangó ignimbrite unit. Zircons from the SZVG locality in the same area indicate that this volcanic bed can belong to the youngest, 14.4 Ma Harsány ignimbrite unit, as do the slightly reworked ignimbrite occurrences in the northwest margin of the Bükk Mts. (e.g., SZIV locality). The thick ignimbrite at Tar, which was formerly regarded as the type-locality of the Middle dacite tuff horizon, has an unambiguous affinity to the 14.9 Ma Demjén ignimbrite corroborating the conclusion drawn based on trace element signature of the glass shards (Harangi et al., 2005; Lukács et al., 2018a). In summary, a wide variation of volcanic eruption products are exposed in the Northern Pannonian Basin, but zircon fingerprints help to find their relations.

5.1.2 Eruption Chronology

Combination of zircon U-Pb dates and trace element composition provides a firm tool to correlate scattered eruption products of large volcanic events. The new zircon geochronological data, in addition to the trace element characteristics of the distal silicic pyroclastic deposits in the northern Pannonian Basin, help to constrain better their ages and their relations to the main ignimbrite-forming eruption events.

The ignimbrites at Ipolytarnóc and Nemti were considered to belong uniformly to the oldest eruption unit in the Pannonian Basin, the so-called Gyulakeszi Formation (Hámor et al., 1980; Hámor 1985). K-Ar ages from these sites (Hámor 1985) fall in the wide age range of 19.0 to 17.5 Ma characterizes this formation although the oldest ages were preferred. However, Pálfy et al. (2007) provided ID-TIMS zircon age data for the Ipolytarnóc ignimbrite that was younger (17.42 ± 0.04 Ma) than the formerly postulated age of ∼19.0 Ma. They published also ⁴⁰Ar/³⁹Ar plagioclase dates of 17.02 ± 0.14 Ma for Ipolytarnóc and 16.99 ± 0.16 Ma for Nemti. Our new zircon ages support the younger age of the Ipolytarnóc and Nemti ignimbrites (17.5–17.2 Ma and 17.1 Ma, respectively). Furthermore, trace element composition of zircons suggests that they could derive from two distinct eruptions: the Ipolytarnóc ignimbrites have zircon trace element affinity to the 17.5 Ma Eger ignimbrite unit, whereas the unwelded ignimbrite at Nemti has zircon trace element features akin to the slightly younger, 17.1 Ma Mangó ignimbrite eruption unit. So far, no older eruption ages (based on zircon U-Pb chronology) have been obtained for silicic volcanoclastic rocks in this area, thus the general assumption being the oldest volcanoclastics remains true, but with a younger formation age. Although the zircon U-Pb ages and trace element compositions strongly suggest a connection between these ignimbrite occurrences, paleomagnetic data seem to yield a different interpretation. Márton and Pécskay (1998) and Márton et al. (2007b) used the paleomagnetic declination data to distinguish three major eruption events, following the classic three-fold subdivision of volcanoclastics in the northern Pannonian Basin (Hámor et al., 1980). This approach assumes very similar (identical) rotation values for the coeval volcanic units and vica versa; similar rotation would mean similar formation age. Structurally speaking, this model would induce homogenous rotation of completely rigid blocks or microplates (Márton and Fodor, 1995). Our results are at odds with this model in a few sites; e.g., the ∼30° counterclockwise rotation of Ipolytarnóc ignimbrite would imply the same rotation as the 16.8 Ma Bogács unit and therefore a coeval formation (Márton and Márton, 1996; Márton et al., 2007a), whereas the 17.1–17.5 Ma eruption units (Mangó and Eger) show larger rotation. Further complex studies are needed to dissolve this apparent controversy and the solution may imply heterogeneous block rotations or local remagnetization.

Zircon trace element characteristics of the Mátraszele bentonite do not give clear inference for its origin. The youngest zircon age population with the inferred eruption age of 16.67 Ma is close to the Bogács unit, yet this is not confirmed by the zircon trace element features. The relatively large compositional variation and different U-Pb dates could be consistent with the reworking character of the bentonite body.

The studied pyroclastic rocks in the northern Pannonian Basin confirm a relatively long gap in the silicic volcanism between 16.2 Ma and 14.9 Ma, which was already recognized in the Bükkalja volcanic field (Lukács et al., 2018a). One of the largest eruptions could have been related to the formation of the Demjén ignimbrite (14.880 ± 0.014 Ma; Lukács et al., 2015; Lukács et al., 2018a). The two ignimbrite occurrences at Tar can be clearly correlated with this eruption as noted already by Harangi et al. (2005) based on glass trace element composition and later by Lukács et al. (2018a) based on bulk rock geochemistry. This eruption age is older what Zelenka et al. (2004) proposed (13.0–13.9 Ma) and younger what was obtained previously for the TAR-Gh locality (15.9 ± 0.6 Ma) using K-Ar dating. Zircon trace element composition of both localities at Tar suggests connection to identical eruption event and a close correlation with the Demjén ignimbrite unit in the BVF. These volcanic products along with the thick pyroclastic sequence at Sirok as well as at Demjén (sites 5 and 8 on Figure 1, Lukács et al., 2018a) can define the lithostratigraphic unit of the Tar Dacite Tuff, but their age suggests a younger origin than previously determined by K/Ar technique (17–16 Ma, Hámor et al., 1980; Gyalog and Budai, 2004, Pécskay et al., 2006). Our results presented in this paper strengthen the conclusion of Lukács et al. (2018a) that the Demjén ignimbrite eruption event at ∼14.9 Ma yielded the Tar Dacite Tuff, which is the most extensive silicic pyroclastic horizon and a significant stratigraphic marker layer in the Pannonian Basin and probably also in Europe.

The Demjén ignimbrite eruption was followed by eruption of a more silicic magma yielding the Harsány ignimbrite at 14.36 Ma (Lukács et al., 2007; Lukács et al., 2015). Two pyroclastic deposits found several 10’s kilometres from one another (SZVG and SZIV) have zircon trace element composition readily showing a strong affinity to this eruption unit and this is confirmed also by the determined zircon eruption ages within uncertainty (14.4 ± 0.2 Ma and 14.3 ± 0.2 Ma). This suggests that the Harsány ignimbrite eruption could have also regional impact even in far distances as implied by the volcanic ash occurrences with similar ages in central and southern Italy (Wotzlaw et al., 2014; Lukács et al., 2018a; Rocholl et al., 2018). Trace element fingerprint of zircons could help to prove this correlation.

Petrogenetic Implications

5.1.3 Zircon/Melt Partition Coefficients

Trace element data of zircons can be used for tephra correlation purposes only when compositionally distinct magmas erupted. Trace element composition of zircons reflects the magma composition as well as the crystallization condition in the shallow magma reservoir such as the prevailing temperature and redox state. This is particularly true for crystal-poor evolved systems when zircon is a late phase and is in equilibrium with the melt (represented by glass phase in the erupted material). Knowledge on the mineral/melt distribution or partition coefficients could help to constrain better the chemical evolution of the silicic magma as well as characterizing the coexisting melt when no available glass phase exists in the rock. The partition coefficients are determined experimentally, although natural mineral and glass data from volcanic rocks can be used also to calculate these values although with inherently less precision than the former data. Nevertheless, the latter approach appears to be more promising as pointed out by Claiborne et al. (2018) and Burnham (2020). Partition coefficients between zircon and melt (glass) were published for dacitic to rhyolitic magmas from various tectonic settings (e.g., Nagasawa, 1970; Mahood and Hildreth, 1983; Fujimaki, 1986; Sano et al., 2002; Bachmann et al., 2005; Colombini et al., 2011; Padilla and Gualda, 2016; Claiborne et al., 2018) and these data show a considerable scatter (Figure 5). The Miocene silicic ignimbrites from the BVF have low crystal content and contain fresh glass shards and pumices and therefore they provide a good opportunity to determine zircon/melt partition coefficients for calc-alkaline rhyodacitic to rhyolitic systems. For such purposes, we used the trace element content of zircons measured at the outer margins of the crystals and combined these data with trace element composition of the glass phases in the same sample (Harangi et al., 2005) to calculate the partition coefficients for a set of trace elements.

FIGURE 5. Zircon-melt partition coefficients for the BVF rhyolites (left panel) compared to a set of published data (gray lines) for silicic volcanic systems: High-K calc-alkaline rhyolite from Highland Range HRL21 and Peach Spring Tuff (KPST01), Colorado River Extensional Corridor (CREC), western United States (Colombini et al., 2011), High-K dacite, Fish Canyon Tuff (Bachmann et al., 2005); dacite from Spirit Lake stage, Mount St. Helens (MSH; Claiborne et al., 2018); Torihama dacite pumice, Kyushu, Japan (Sano et al., 2002) and Króksfjördur calc-alkaline rhyolite from Iceland (Claiborne et al., 2018). Zircon/melt partition coefficients for the BVF rhyolites (gray field) are normalized to the average values represented with the ±standard deviation range (pink field right panel). They are mostly within the range of ±2-times of the average values (green field). A set of published partition coefficients are also presented for comparison.

Determination of zircon/melt partition coefficients was performed in a similar way as Claiborne et al. (2018). Zircon trace element data were carefully screened before using them for partition coefficient calculation. Only the outer rim data were selected and any outliers possibly reflect mineral or glass inclusions were omitted. We used the average along with the standard deviation of the remaining data set to characterize the zircon trace element composition assumed to be in equilibrium with the host melt represented by the glass data. In each eruption unit, only those samples where trace element content of the glasses was homogeneous were involved. Thus, samples from the Bogács unit were not involved, where different juvenile clasts exist and the zircon-melt relationships are not obvious. On the other hand, two compositional types of zircon are recognized in the Demjén unit and this seems to be consistent with the slight compositional difference in the glasses (particularly in the Eu-anomaly; Supplementary Figure S6 in Supplementary Data Sheet S2). Here, only the high Yb/Dy and high Eu/Eu* zircons were considered and used along with the glasses showing no or slight negative Eu-anomaly. These partition coefficients were similar to the ones obtained for the other units. In a few cases, partition coefficients were calculated also using bulk rock data (Lukács et al., 2018a) and these gave similar results within uncertainty (Supplementary Figure S4 in Supplementary Data Sheet S2). The result is presented in Table 4 and shown in Figure 5. The quality of the obtained zircon/melt partition coefficients was evaluated based on the lattice strain model (Blundy and Wood, 1994; Blundy and Wood, 2003a; Blundy and Wood, 2003b) using Onuma diagrams (Onuma et al., 1968).

TABLE 4. Calculated zircon/melt partition coefficients with the obtained strain-lattice model parameters for samples of three eruption units in the BVF. The partition coefficients have very similar values although zircons show distinct trace element characteristics. The average partition coefficients for the BVF silicic volcanic rocks are also given. Detailed data are presented in Supplementary Table S3.

The calculated zircon/melt partition coefficients (D_zircon/melt) of the Bükkalja rhyolites have very similar values to published data for dacitic to rhyolitic systems (Figure 5; Sano et al., 2002; Bachmann et al., 2005; Colombini et al., 2011; Padilla and Gualda, 2016; Claiborne et al., 2018). Remarkably, the same values were obtained within the ±standard deviation range, irrespective of the eruptive units even though zircon, glass and bulk rock compositions differ significantly in the BVF rocks. This implies that chemical composition of the Si-rich magmas (at SiO₂ >70 wt%) does not have major influence on the D_zircon/melt’s. In general, the partition coefficients show a steady increase from the light REE to the heavy REE (D_zircon/melt’s for HREE are >100); Th, U and Hf are strongly compatible (D_zircon/melt >10), whereas Nb is slightly incompatible (D_zircon/melt = 0.3–0.7). Although we calculated D_zircon/melt values for La and Pr, we note that these values have to be considered with caution, because both elements have very low concentration, close to the detection limit for LA-ICP-MS analysis in zircons. In fact, their partition coefficients show considerable variation in the published data sets and also in our results. Cerium is slightly incompatible (D_zircon/melt = 0.2–0.3) showing a positive anomaly (D_zircon/melt^Ce > D_zircon/melt^Ce*) in the partition coefficient pattern (Figure 5), whereas Eu has a slight negative anomaly similarly to the Mt St Helens data (Claiborne et al., 2018). The Th/U^zircon/melt data are fairly similar ranging from 0.11 to 0.16. There is also a generally good agreement with the partition coefficients of rhyolites, such as from the Highland Range and the Peach Spring Tuff (Colombini et al., 2011; Padilla and Gualda, 2016) that were formed during the early stages of Miocene extension in the southwestern United States, a similar tectonic setting characterized the Pannonian Basin during the silicic ignimbrite volcanism at the BVF. Furthermore, we got similar, although generally lower, values as those from the only calc-alkaline rhyolite in Iceland at the 11–12 Ma Króksfjördur (Jónasson et al., 1992; Claiborne et al., 2018; Figure 5).

Trace element partitioning between a solid phase and a coexisting melt is described by the crystal lattice-strain model (Blundy and Wood, 1994; Blundy and Wood, 2003a; Blundy and Wood, 2003b) providing that in equilibrium condition trace element distribution coefficients are strongly dependent on the ionic radius and ionic valence that can be described by a mathematical expression. This means that isovalent cations have a parabola-like function in an ln D_min/liq vs. Ionic radius plot, commonly referred as Onuma diagram (Onuma et al., 1968). The parabola-like curves have four freely varying parameters: the maximum D_min/liq (D_o(M)ⁿ⁺), the radius of that site (r_o(M)ⁿ⁺), the apparent Young ‘s Modulus (E_Mⁿ⁺) and the temperature, which make it difficult to find the best-fit result. This is also due to the fact that trivalent cations (REE and Y) fall on the right side of the parabola, while there is no strict constrain for the left side due to the lack of trivalent cations (D_zr/liq^Sc could help in this, but it is usually not determined, because of difficulties to measure Sc in zircon, particularly by LA-ICP-MS). In our calculation, we did not include La, Ce, Eu in the line-fitting and constrained the temperature between 700 and 750°C. The best-fit curves are presented in Figure 6 and the obtained parameters are shown in Table 4. In all cases the residuals, i.e., ln (D_min/liq^fit/D_min/liq^calc) as proposed by Colombini et al. (2011), remained low (close to zero).

FIGURE 6. Onuma diagrams showing the best-fit parabola-like curves on the obtained zircon/melt partition coefficients for the Harsány ignimbrite; Demjén ignimbrite; and Mangó ignimbrite. They have very similar parameters in the equation of the lattice strain model of Blundy and Wood (1994), the values are presented in Table 4. The best-fit curves are compared with those for published data (grey-lines; S02: Sano et al., 2002; B05: Bachmann et al., 2005; C11: Colombini et al., 2011; bottom right panel and general similarities can be recognized, particularly with that of Sano et al. (2002).

The best-fit curves with high R² (>0.92) for the Bükkalja samples (Figure 6) support the validity of the obtained zircon/melt partition coefficients. The position of the right part of the parabola is the same compared to other published data (Figure 6; e.g., Sano et al., 2002; Colombini et al., 2011), only the Fish Canyon Tuff line (Bachmann et al., 2005) deviates a bit. This suggests that relatively coherent zircon/melt partition coefficients characterize the calc-alkaline dacitic to rhyolitic systems and they are not significantly influenced by the silica content of the melt at SiO₂>70 wt%. In the Onuma plots, the r_o(M)^{3 +} values are similar (0.920–0.935) in all cases that are a little bit lower than what Blundy and Wood (2003a) and Burnham and Berry (2012) determined based on experimental D_zircon/melt values (0.968–1.1018). The E_M³⁺ is also in a fairly similar range (520–793 GPa) consistently with the obtained values by Claiborne et al. (2018), whereas the D_o(M)^{3 +} values range from 328 and 575. These values are in the range listed by Burnham and Berry (2012) for natural samples. Ce has a significant positive offset from the best-fit line, whereas Eu has a modest negative offset what are due to the coexistence of Ce³⁺ and Ce⁴⁺, and Eu³⁺ and Eu²⁺, respectively. The extent of the deviation is dependent on the oxidation state of the magmatic system.

5.1.4 Constraints on the Magma Evolution

Zircon trace element characteristics reflect closely the environment under the crystallization took place. Combining these data with the U-Pb dating results, important constraints on the evolution of the crustal magma reservoirs can be made (e.g., Claiborne et al., 2010; Schmitt et al., 2010; Reid et al., 2011; Klemetti et al., 2011; Wotzlaw et al., 2013; Chamberlain et al., 2014; Cooper and Kent, 2014; Klemetti and Clynne, 2014; Storm et al., 2014; Barboni et al., 2016; Reid and Vazquez, 2017; Szymanowski et al., 2017). Zircons from the BVF silicic volcanic rocks show distinct trace element compositional features (Figure 2) implying different crystallization conditions in their reservoirs before eruption. Noteworthy, these magma reservoirs existed and evolved partially coeval as shown by their U-Pb spot dates (Lukács et al., 2015; Lukács et al., 2018a; Harangi and Lukács, 2019).

Zircon Th/U ratio can be used as a proxy for the magma differentiation, whereas Hf and Ti reflect the thermal history of the magma during crystallization (Claiborne et al., 2006; Claiborne et al., 2010). The Th/U ratio in zircons depends closely on the chemical composition of the melt from which they crystallize (Belousova et al., 2002; Schmitt et al., 2003; Bindeman et al., 2006; Reid et al., 2011; Kirkland et al., 2015; Reid and Vazquez, 2017). Both elements are highly compatible to zircons and they usually show substantial variation as crystallization condition changes (Hoskin and Schaltegger, 2003). As melt is getting to be more differentiated, zircon Th/U ratio decreases, although this can be affected also by co-crystallization of other accessory phases, such as allanite. In the studied BVF rocks, the Th/U ratio is in the range between 2 and 5 both in the whole rocks and the glasses (Harangi et al., 2005; Lukács et al., 2018a). In the zircons, this ratio is <1.1, mostly between 0.3 and 0.8, which is typical of igneous zircons (Hoskin and Schaltegger, 2003), although the Highland Range zircons show higher values (Colombini et al., 2011). The corresponding D_zr/melt^(Th/U) is around 0.14 ± 0.02 consistent with zircons grew from high-silica rhyolitic melt (Bindeman et al., 2006; Reid and Vazquez, 2017; Supplementary Figure S5 in Supplementary Data Sheet S2). Only the hybrid, less evolved pumice and zircon pairs from the Bogács unit have higher D_zr/melt^(Th/U) values (0.2–0.3) suggesting crystallization from less evolved, dacitic melt (Schmitt et al., 2003; Bachmann et al., 2005; Bindeman et al., 2006). The progressively decreasing Th/U ratio in the zircons could reflect the co-precipitation of allanite, which is particularly visible in the Harsány and Mangó ignimbrites (Lukács et al., 2007). Appearance of allanite affects not only the Th/U ratio, but also the light REE to heavy REE ratio in the melt as recognized this in the Harsány ignimbrite (Lukács et al., 2007).

The Th/U ratios do not show any correlation with Hf, neither Ti suggesting that crystallization of main and accessory mineral phases controlled the Th/U ratio independently from temperature. Indeed, zircon saturation temperature calculated from the whole rock and glass composition following the procedure of Boehnke et al. (2013) shows a restricted range between 700 and 750°C (bulk rocks) and 640–740°C (glasses) for all the samples from different eruption units except for the Bogács unit, where we obtained larger temperature, mostly above 750°C (Figure 7). The low crystallization temperature of zircons is corroborated by the low Ti concentration of zircons (Ti < 7 ppm), while it is between 5 and 14 ppm in the samples from the Bogács unit. On the other hand, there is a negative correlation between Ti and Hf in zircons, supporting the notion that Hf variation is strongly dependent on crystallization temperature (Claiborne et al., 2010). Hf content is typically low in the Bogács unit (<9500 ppm), whereas it is always greater than 9000 ppm in the other eruption units (Figure 2). Zircons from the Mangó and Harsány ignimbrites have relatively wide variation in Hf (from 8300 to 12000 ppm), whereas the Eger and Demjén samples show usually much restricted range (mostly from 9000 to 11000 ppm).

FIGURE 7. Ti-in-zircon temperature Uncertainty of single temperature values is shown as ±1 sigma range variation in function of Eu-anomaly (expressed as Eu/Eu*) as an indicator of differentiation in the BVF major units. Temperatures are calculated based on the thermometer of Ferry and Watson (2007) using a_TiO2 = 0.6. Legend is the same as for Figure 2. The colored boxes represent the zircon saturation temperatures with glass and bulk rock composition data (Harangi et al., 2005; Lukács et al., 2018a) calculated using Boehnke et al. (2013) scheme for the crystal-poor rhyolites (Harsány, Demjén, Mangó, Eger) and for different clasts types of the Bogács unit (Czuppon et al., 2012).

The relatively low Hf, high Th/U and high Ti content of zircons in the Bogács unit is consistent with repeated replenishments by less evolved mafic to intermediate magmas in the crustal reservoir (Czuppon et al., 2012). The Ti-in-zircon thermometer (Ferry and Watson, 2007) yields large range of temperature values from 710 to 840°C for the Bogács zircons (at Ti activity of 0.6; Figure 7). The low temperature values (710–760°C) overlap the temperature range obtained for the zircons of other units (680–750°C) suggesting crystallization from an evolved melt in a silicic crystal mush environment. This is supported also by the low zircon saturation temperature values calculated from bulk rock and glass compositions (boxes in Figure 7). Thus, the Bogács magma reservoir was thermally heterogeneous, whereas in the other units zircon developed at relatively narrow temperature interval close to the solidus. In the Bogács unit, various juvenile clasts with less evolved scoria and A-fiamme glasses and rhyolitic B-fiamme and pumice glasses occur in the upper eruptive subunit implying complex mixing and mingling processes prior to the eruption as described in detail in the petrogenetic model of Czuppon et al. (2012). They assumed an extended magma reservoir where melt-dominated lenses with slightly different chemical compositions were formed. Mixing of melts and crystal mush material occurred periodically during recharge events. Although the analyzed zircon crystals were separated from single clast types, the crystallization environment of zircons is not straightforward. Applying the zircon/melt partition coefficients presented in this study, it is possible to calculate the melt composition, which was in equilibrium with the zircons (Figure 8). Some of the modeled melt composition resemble that of the most evolved, rhyolitic glasses of pumices, whereas the others show similarities to the glasses of the A-fiamme (Czuppon et al., 2012). These melts were generated in distinct pods, where the pumice glasses represent the most differentiated melt remained after crystal fractionation, whereas the A-fiamme melt was formed after mixing of dacitic and rhyolitic magmas. Subsequently, these distinct melt fractions were incorporated into the erupted magma just prior to the eruption.

FIGURE 8. Composition of various glass types, high-Si rhyolitic pumice glass and glass shards (left panel) and less evolved glass from A-fiamme (Czuppon et al., 2012 right panel), compared to the modeled melt composition (gray field: average value with standard deviation range) calculated from the zircon trace element data applying the average partition coefficients for the BVF silicic rocks presented in this study (Table 4).

Zircon trace element data could help also to refine the magmagenesis of the 14.9 Ma Demjén ignimbrite, one of the largest volume eruption products within the BVF. The bulk rock petrology and composition as well as glass trace element content seem to be homogeneous and distinct compared to the ignimbrites of the other units (Harangi et al., 2005; Lukács et al., 2018a). Occurrence of amphibole and the depletion in the Y and heavy REE suggest that amphibole crystallization had a major role in their magma evolution, a feature not present in the other units. However, zircon trace element concentrations show a relatively large variation and two groups can be distinguished, particularly based on the Eu/Eu* and the Yb/Dy ratios (Supplementary Figure S6 in Supplementary Data Sheet S2). This could be consistent with the presence of two melts with slightly distinct compositions in the magma reservoir and these melts reflect different degrees of amphibole and plagioclase fractionation. Melt trace element contents calculated from the zircon compositions using the average Kd value for the BVF rocks provided in this study give slight differences that can be observed also in the glass trace element compositions of the Demjén ignimbrite (Supplementary Figure S6 in Supplementary Data Sheet S2). Interestingly, zircons from the much older, 17.5 Ma Eger ignimbrite overlap exactly the less evolved zircon group of the Demjén ignimbrite in the trace element ratio plots (Figure 9).

FIGURE 9. Zircon compositions of the BVF eruptive units show distinct evolutionary paths. The Harsány, Mangó and Bogács units are characterized by relatively low Yb/Dy ratio and wide variation in Dy consistent with dominant plagioclase crystallization accompanied with minor allanite, zircon and amphibole fractionation (left panel), whereas the Demjén and Eger zircons have high Yb/Dy ratios and limited Dy variation what can be explained by dominant amphibole fractionation in addition to plagioclase crystallization (right panel). Occurrence of titanite could also increase strongly the Yb/Dy ratio and decrease Dy (e.g., Wotzlaw et al., 2013; Szymanowski et al., 2017), but titanite has not been observed in any BVF samples.

6 Conclusion

In this study, we demonstrate that zircon trace element content together with their U-Pb dates is a powerful tool to correlate eruption products and to constrain the magmatic differentiation of silicic magmas. During the early to mid-Miocene, a highly productive magmatic event, which produced more than 4000 km³ of erupted material (including far-traveled ignimbrites) occurred in the Pannonian Basin. These volcanic products serve as unique stratigraphic marker horizons in central and southern Europe. The newly determined zircon U-Pb eruption ages for the distal, previously ambiguously categorized pyroclastic deposits in the northern Pannonian Basin are between 17.5 and 14.3 Ma, which are comparable with published ages for the main eruptive events. The five major eruption units identified in the Bükkalja volcanic field have zircon with a characteristic and distinctive geochemical fingerprint. Multivariate discriminant analysis proved to be useful in distinguishing these main eruptive units based on selected trace elements and trace element ratios of zircon and to link the distal occurrences to the already defined eruption units. The new results give a better understanding on the temporal and spatial distribution of the volcanic products of the ignimbrite flare-up episode. This methodology with zircon fingerprinting offers a promising perspective to be applied for tephras and pyroclastic rocks where glasses as well as the main minerals are already altered or missing and only zircon preserves the original magma character. This is the case for instance for the far-distal ash deposits in various Paratethys sedimentary sequences in central and southern Europe (Lukács et al., 2018a; Rocholl et al., 2018). Revealing the zircon compositional similarities of ash beds with similar ages can confirm their common origin, which can lead to the refinement of the chronological framework for these early to mid-Miocene sediments.

Using our data set, and published glass trace element data, we determined zircon/melt partition coefficients for various volcanic units and found very similar values to published data for other volcanic systems. This suggests that zircon/melt partition coefficients in calc-alkaline dacitic to rhyolitic systems are not significantly influenced by the SiO₂ content of the melt at SiO₂ > 70 wt%. The obtained partition coefficients and zircon trace element data can be used to calculate the equilibrium melt composition. Since tephra correlation is mostly based on glass composition, this is particularly important for volcanic material in which glasses are thoroughly altered or missing. Nevertheless, the zircon proxy approach can be limited in cases when no systematic compositional difference is found between eruption products, although the latter problem similarly stands for glass chemistry-based tephrostratigraphic studies.

During the 3 Myr long Miocene magmatic phase in the Pannonian Basin, silicic magmas formed during a period of lithospheric extension. Magma reservoirs with distinct composition existed contemporaneously, but evolved at least partly separated. Most zircons from the explosive volcanic products crystallized from an evolved melt close to the haplogranitic solidus at relatively low temperature (<750°C) for protracted periods of time. Zircon trace element variation implies that amphibole fractionation played an important role in the evolution of the Eger and Demjén magma reservoirs. In contrast, it has a smaller or no contribution in the other units, where allanite crystallization can be recognized leading to contrasting evolved melt compositions. The calculated equilibrium melt data along with the already published glass composition from the same samples show heterogeneity in certain eruptive units. This confirms the view that silicic magma reservoirs are often heterogeneous, not only thermally, but also compositionally, since separated melt pods with different chemical signatures can develop in the magma reservoir as pointed out in Yellowstone and Taupo (e.g., Cooper et al., 2012; Ellis et al., 2014), among others.

Data Availability Statement

All datasets presented in this study are included in the article/Supplementary Material.

Author Contributions

RL designed the study as part of her NKFIH projects. RL performed all the geochemical analyses with the assistance and technical guidance of MG. RL processed and interpreted data and wrote the first draft of the manuscript. SH and RL performed the statistical analysis. LF took part in the regional geological interpretation. OB and SH participated in the petrogenetic discussions and helped finalizing the ideas. SH provided figures and wrote subsections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This study was financed by the National Research, Development and Innovation Office–NKFIH within the FK OTKA project (No. 131869) and partly within two postdoctoral projects for RL (PD112584 and PD 121048). RL was supported by the Bolyai János Research Fellowship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to József Pálfy and Orsolya Sztanó (ELTE), Péter Prakfalvi (MBFSZ) and Imre Szarvas (Bükk National Park) for their help in the field during sample collection. Discussions with Sándor Józsa and József Kovács (ELTE) are appreciated in the petrographic and the multivariate statistical interpretations, respectively. Support of the members of MTA-ELTE Volcanology Research Group is acknowledged. J-F Wotzlaw is appreciated for the ID-TIMS measurements of reference material AUSZ8-10. Constructive comments of the reviewers helped us to clarify our ideas and we thank Adelina Geyer for the careful editorial work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.615768/full#supplementary-material.

References