Introduction

Silicic volcanism not only poses serious threat to life, socio-economics and aviation, but also contributes to long-term changes in climate and global temperature (Self and Blake 2008; Robock et al. 2009; Sigl et al. 2015). Constraining magma evolution and eruption triggers of silicic volcanoes is one of the key objectives in petrology. Among others, one of the crucial aspects in volcanic studies is the state and duration of storage of magma and its crystals. However, despite the gravity, the understanding of this aspect remains fairly deficient (National Academies of Sciences 2017). Various studies have explored and developed different chronometric methods in the last two decades to address this issue, but all possess certain adavantages as well as disadvantages over each other. U-series absolute dating of magmatic crystals robustly constrains timescales of magma storage, transport and evolution (e.g. Zellmer et al. 2000; Hawkesworth et al. 2004; Schmitt et al. 2010; Rubin et al. 2016; Cisneros de León and Schmitt 2019; Cisneros de León et al. 2019 and many others; see reviews and examples in Reid 2008; Schmitt 2011; Cooper 2015). However, this method does not allow to distinguish between times of residence of eruptible magma and the age of crystals that may spend a major part of their lifetime in uneruptible crystal mushes. U-series dating also lacks temporal resolution when it comes to short timescales of eruption triggers just prior to eruption. On the other hand, diffusion chronometry (e.g. Morgan et al. 2004; Chamberlain et al. 2014; Kahl et al. 2015 and many others) provides high temporal resolution, but suffers from higher degree of uncertainties (Spear 2014; Rout et al. 2020).

Diffusion modeling measures the duration of diffusive homogenization of compositional contrasts formed within the crystal due to changes in growth (or magmatic) conditions (see Costa et al. 2008 and Dohmen et al. 2017 for many examples). However, since the rate of diffusion or diffusivity is highly susceptible to temperature (and sometimes moderately by other parameters also e.g. pressure and H2O-content) (Costa et al. 2008; Spear 2014; Dohmen et al. 2017) constraining temperature and the thermal history of magma and its crystals becomes exceedingly important, if storage and evolution timescales are to be estimated.

When it comes to thermal conditions of magma storage, contrasting models have been proposed over the years making it a highly debated topic in volcanology (National Academies of Sciences 2017). One model, known as the “warm storage” model (Barboni et al. 2016), states that magma resides in a reservoir in a hot (> 750 °C) liquid-rich eruptible state for tens to hundreds of thousands of years (Kaiser et al. 2017). The temperature is maintained by periodic recharges of fresh and hot magma (Gelman et al. 2013; Huber et al. 2012). The other model, known as the “cold storage” model (Cooper and Kent 2014; Rubin et al. 2017), argues that instead, magma is stored as a cold (at near or subsolidus temperature) rigid non-eruptible crystal-rich mush for up to hundreds of thousands of years (Reid 2008; Cooper 2019). This kind of a system needs a process that would melt the mush to produce eruptible magma and trigger an eruption (Bachmann and Bergantz 2008; Girard and Stix 2009; Wolff et al. 2015). There is evidence for both models and different chrononmetric studies have fallen into either of the two catagories. For example, for some systems, the diffusion times obtained for temperatures of more than 750 °C are limited to only a few hundreds to a few thousands of years, which are very short (by 100 times) compared to the absolute U-series ages (e.g. Cooper and Kent 2014; Rubin et al. 2017). This suggests that the actual temperature of storage was much lower making it a case of cold storage. Similarly, for some systems, diffusion times obtained for temperatures as high as 780 °C are hundreds of thousands of years in magnitude (e.g. Morgan and Blake 2006). This clearly suggests a very long storage at high temperature making it a case of warm storage.

Here, we study the Laacher See volcano in the context of applicability of concepts of “cold” and/or “warm” storage. We study the pumice and the cognate crystal-rich cumulate ejecta from this phonolitic system. The Laacher See volcano erupted 12.9 ka (Baales et al. 2002) and is one of the largest Quaternary eruptions in central Europe. Although the current level of volcanic and magmatic activity in the crust below the Eifel volcanic field is assumed to be very low, increased gas flux (Goepel et al. 2015) and deep low-frequency earthquakes at 10–40 km depth are detected underneath Laacher See, suggesting ongoing lower-crustal level magmatic activity (Hensch et al. 2019). It is therefore essential to know the triggers as well as the type and timing of potential precursor activity in case Laacher See becomes active again in furture.

The Laacher See reservoir has been well studied in terms of chemical stratigraphy, evolution of the resident phonolite and its crystalline equivalent i.e. the cognate crystal-rich cumulate ejecta (e.g. Wörner and Schmincke 1984a, b; Wörner and Wright 1984; Tait et al. 1989; Harms and Schmincke 2000; Schmitt et al. 2010 and others). For some of these rocks, absolute crystal ages have been determined by U-series dating (Bourdon et al. 1994; Schmitt et al. 2010) which show that crystallization preceded eruption by up to 20 ky. Here we apply diffusion chronometry to diffusive compositional gradients across sanidine crystals from the phonolitic pumice and the cumulates to extract information on crystal residence times and the storage conditions of the main phonolite body (“warm” storage) and the “colder” crystalline cumulate mush (near-solidus storage) prior to eruption.

Geology and petrography

Geological setting

The Laacher See volcano, located in western Germany, erupted 12.9 ka (Baales et al. 2002) and is the youngest and largest volcano in the East Eifel Volcanic Field (EEVF). It is part of the central European Upper Tertiary to Quaternary intra-plate alkaline volcanism, and is located in the western part of the uplifted Paleozoic Rhenish Massif (Fig. 1) (Schmincke 2007, 2008).

Fig. 1
figure 1

Schematic map (modified after Rout and Wörner 2018 and van den Bogaard and Schmincke 1984) of East Eifel volcanic field with Laacher See (~ 12.9 ka) in the east and Rieden (~ 460–360 ka) and Wehr (~ 215 ka) centers in the west and north-west. The green areas represent Laacher See Tephra thickness (in m)

The Upper Tertiary to Quaternary magmatism in this region primarily produced mafic magmas with composition ranging from basanite and alkali basalt to quartz tholeiite as small scattered eruption centers or as large shield volcanoes (e.g. Vogelsberg, 15–16 Ma, 600 km3) (Wilson and Downes 1991; Wedepohl and Baumann 1999; Bogaard and Wörner 2003). There are > 100 eruption centers in EEVF over the past ~ 440 ky (Duda and Schmincke 1978; Wörner and Schmincke 1984a; van den Bogaard et al. 1989a, b; van den Bogaard and Schmincke 1990; Schmincke 2007). Poorly preserved oldest ash deposits erupted in the NW of the EEVF date back to 450–650 ka (van den Bogaard and Schmincke 1990).

The eruption history of EEVF consists of mostly varieties of mantle-derived alkaline magmas (predominantly basanite, nephelinite, leucitite). Mafic eruptions generally start with a phreatomagmatic phase forming maar deposits followed by growth of scoria cones and emission of blocky lava flows (Schmincke et al. 1973; Schmincke 2007). Large volumes of evolved phonolitic magmas formed intermittently that are identified to be the differentiation products of the mafic magmas (Duda and Schmincke 1978, Viereck 1984 and Wörner and Schmincke 1984a, b). Evolved magmas formed the Rieden volcanic complex that produced more than 10 km3 of phonolite and leucititic tephra 380–430 ka (Viereck 1984; van den Bogaard et al. 1989a; van den Bogaard and Schmincke 1990). After Rieden, the Wehr volcano erupted phonolite to trachyte and evolved phonolitic tephra ~ 215 and ~ 150 ka, respectively (Wörner et al. 1988; van den Bogaard and Schmincke 1990). The third and the largest evolved volcanic center is the phonolitic Laacher See volcano. The evolution of this phonolitic reservoir began ~ 30 ka (Schmitt et al. 2010) and ultimately fed the gigantic eruption 12.9 ka (Baales et al. 2002).

The Laacher See eruption

The Laacher See produced ~ 6.3 km3 DRE of phonolitic magma within a period of days to weeks (van den Bogaard 1995) whereas the total eruption period may have lasted for several months (Schmincke 2007). The eruption is documented as a major Quaternary tephra marker across central and northern Europe (van den Bogaard and Schmincke 1984; van den Bogaard et al. 1989b; Baales et al. 2002). Millimeter-scale ash layers have been observed in southern Sweden and northern Italy, more than 1000 km from the Laacher See vent (Blong et al. 2018; Hensch et al. 2019). Several petrological and experimental studies have estimated the top of the phonolitic reservoir to be at 4–6 km depth (Wörner 1982; Harms and Schmincke 2000; Berndt et al. 2001; Harms et al. 2004) based on experimental and barometric pressure constraints. The Laacher See reservoir is a classic example of a zoned magma reservoir which is evident from its well-stratified and compositionally zoned tephra deposit (Laacher See Tephra or LST, Fig. 2) (Wörner et al. 1983; Wörner and Schmincke 1984a, b). The Lower Laacher See Tephra (LLST) is characterized by highly differentiated composition and enrichment of incompatible trace elements and represents the top of the reservoir. Middle Laacher See Tephra (MLST) is the moderately differentiated and most volumous mid-part of the reservoir while the Upper Laacher See Tephra (ULST) is the least evolved base of the reservoir. The top most part of LST (above ULST) that was deposited towards the end of the eruption derives from the base of the magma reservoir and contains heterogeneous clasts of a hybrid composition. They contain disequilibrium mineral assembleges with (1) mafic minerals such as forsteritic olivine and magnesian clinopyroxene, even small mantle peridotite xenoliths next to (2) phonolitic minerals like sanidine, sodic plagioclase, haüyne, amphibole, fassaitic clinopyroxene and titanite, which are typical of MLST and ULST. These hybrids are products of mixing and mingling between the resident phonolite and an intruding mafic basanite magma (Fig. 3; Wörner and Wright 1984). It is unclear wheather this mixing was pre-eruptive (perhaps part of the trigger) or rather syn-eruptive (a by-product of the eruption). A SiO2 vs alkalis plot is presented in Fig. 3 highlighting the resident phonolite, intruding basanite and hybrid compositions.

Fig. 2
figure 2

Compositional variation across the Laacher See Tephra (modified after Wörner and Schmincke 1984a). Samples from the whole rock data set that are studied here are noted specifically and marked next to their stratigraphic position

Fig. 3
figure 3

The cumulates and the volcanic samples (including phonolitic pumice, basanite and hybrid samples) plotted on a alkalis vs SiO2 plot (after Cox et al. 1979). The three points for the cumulates denote the representative compositions for the felsic, intermediate and mafic cumulates. The figure also shows the Laacher See phonolite differentiation, basanite-tephrite differentiation and basanite-phonolite mixing trends

Apart from phonolite, the ULST has a variety and abundance of contact-metamorphosed and fenitized country rocks as crustal xenoliths, e.g. partially melted Devonian slates, greenschists and rare lithics from the amphibolite-facies metamorphic basement (Wörner et al. 1982; van den Bogaard and Schmincke 1984). They were incorporated from the thermal aureole around the reservoir during the gravitational destabilization in the final stages of the eruption. The observed thermal metamorphism indicates that the liquid-rich phonolite had continuously supplied heat to the surrounding country rock over its 20 ky long evolution. The results of this heat loss are also evident from the cognate plutonic and cumulate ejecta (also present only in ULST), which represent the cooling exterior of the reservoir. The cumulates are crystal-rich (> 70%), and petrographically and texturally very diverse (see sample description for details) indicating a very complex crystallizion history of a cumulate carapace that surrounded the Laacher See reservoir (Rout and Wörner 2018).

Apart from the cumulates, the ULST also contains, although rarely, highly crystalline cognate carbonatitic syenites. These carbonatites are distinctive in terms of their intrusive nature and presence of interstitial highly differentiated phonolitic glass similar to LLST (Liebsch 1996). Like the cumulates, these rocks are also interpreted to be fragments of the surrounding crystalline carapace (Tait et al. 1989). Crystals in these rocks share highly faceted idiomorphic boundaries with the interstitial vesicular melt. The interstitial glass, like LLST, is also poor in Ba and Sr, and highly enriched in incompatible trace elements. These observations argue against re-melting of the carbonatitic syenites to generate the interstitial melt and, instead, suggest a close genetic relation with the main phonolitic magma from the reservoir interior. U-series dating of cumulate minerals and zircons contained in the carbonatitic syenites (Schmitt et al. 2010) indicate, that this evolved crystal-rich carapace of the Laacher See phonolite magma existed for about 20 ky prior to the 12.9 ka eruption. Rout and Wörner (2018) specifically focussed on these carbonatitic syenites and further constrained their residence and thermal history. Their study concludes that the syenites were stored as a crystalline mushy carapace for up to 20 ky at 630–670 °C until the final destabilisation of the system that occurred only 15–50 days before the eruption.

However, the thermal history and eruption trigger of the main volume of erupted phonolite still remained unknown. Therefore, in the present study, we analyze (1) sanidine phenocrysts from phonolitic pumice representing the molten core of the magma system, (2) sanidines in the cumulates and (3) sanidines in the final hybrid eruption products to fully document the thermal history of the entire Laacher See magma system prior to its eruption.

Samples and petrography

Pumice

We have investigated sanidine phenocrysts from seven pumice samples that were chosen as representatives from different parts of the compositionally zoned deposit (Fig. 2). Figure 2 shows the chemical gradient across LST including the position and whole rock composition of the samples studied here.

Detailed mineralogical descriptions of these samples have been given by Wörner and Schmincke (1984a). The sample (1002) from LLST are white, micro-vesicular and crystal-poor pumice from the very base of the deposit (top of the magma chamber). Rare phenocrysts (≤ 1%) are mainly sanidine but nepheline, plagioclase, haüyne, clinopyroxene, amphibole, titanomagnetite, titanite, apatite, rare zircons and cancrinite also occur in traces. Samples from the middle and top of the MLST (1034, 1050: middle and base of the central part of the magma chamber) are white to light grey pumices. Phenocryst-content varies from rare (1–5%) to relatively more abundant (5–15%) up-section (i.e. down the magma body). They are mainly feldspars (sanidine and plagioclase), large euhedral haüyne, titanite, apatite, clinopyroxene, amphibole and titanomagnetite. Clots, 1–3 cm in size, of a more mafic phonolite with crenular margins, petrographically identical to the more mafic pumice found in ULST (i.e. more mafic phonolite with 57–58 wt% SiO2) occur in various states of disaggregation. Samples from the ULST (1060, 1088) are mafic phonolites from the base of the magma chamber that have a phenocryst content of 30–50%, mainly feldspars (plagioclase and sanidine), amphibole, clinopyroxene, titanite, apatite and titanomagnetite. There is no evidence for significant mixing with a mafic magma (only few small phlogopite crystals and rare Fo-rich olivines are found in the matrix of the most mafic phonolite). Samples 1088-H1, 1088-H3 and 1099, which come from the very top of the deposit, are hybrid between the phonolite and the basanite (Wörner and Wright 1984), and contain, in addition to mineral phases similar to those of phonolite sample 1088, significant amounts of crystals from the basanitic endmember: olivine, pyroxene megacrysts and phlogopite. However, feldspars are the most abundant phase (~ 30%, up to 50%). Large, resorbed nepheline crystals are also present. Because their feldspar mineralogy is the same as in mafic phonolite sample 1088, the feldspars in the hybrid lavas are considered to be derived from (and equivalent to) those observed in crystal-rich ULST magma.

Cumulates

The Laacher See cumulates are crystal-rich fragments of highly variable crystal contents that range from mafic (clinopyroxene- and/or amphibole-rich) to evolved ("sanidinites”) and show a large variability in petrographic textures (Tait 1988). The cumulate clasts are present only within ULST. Tait et al. (1989) have shown, on the basis of the petrography and chemical composition of minerals and interstitial glasses, that these cumulates are comagmatic with the LST phonolite, representing the materials crystallizing at the chamber margins prior to eruption.

Cumulates studied here typically consist of interlocking crystal frameworks with a glassy matrix. Individual crystals are 0.5–5 mm in size and sometimes show parallel alignment textures typical of layered cumulate rocks. Mineral assemblages range from mafic to felsic, and interstitial phonolitic glass occupies between ~ 1 and 30 vol% of the cumulates. Figure 3 shows the three representative cumulate compositions (felsic, intermediate and mafic) in a TAS plot along with other Laacher See products. Mafic cumulate samples consist of predominantly amphibole (25–55 vol%), clinopyroxene (25–50 vol%), magnetite (5–15 vol%) and apatite (< 10 vol%) with small amounts of feldspar (< 5 vol%). An average of 7 vol. % matrix is present with a variable clinopyroxene/amphibole ratio. Intermediate samples contain between 50 and 80 vol.% feldspar, the remainder being composed of clinopyroxene, amphibole, apatite, phlogopite and titanite. The felsic cumulates are commonly of syenitic composition with sanidine occurring as the only type of feldspar. Haüyne (10–15 vol%) and minor amounts of clinopyroxene, amphibole, titanite and magnetite are also generally present. Some of the other felsic samples contain sanidine plus a triclinic feldspar which is usually anorthoclase but sometimes lies in the adjacent oligoclase field of the plagioclases (Tait 1988). The matrix composition in the felsic and intermediate cumulates span the same range as MLST and LLST pumice matrix (i.e. mafic to evolved phonolite, Tait et al. 1989). Full details of the petrography of cumulates are described in Tait (1988) and Tait et al. (1989).

Analytical methods

Thin sections from seven individual pumice samples were investigated, whereas mineral separates were used for samples 1002 and 1034 because of their low crystal contents. In each of the samples, 15–50 sanidine crystals were investigated by optical microscopy and 4–5 zoned crystals (30 in total) were chosen for detailed microprobe analyses. From the cumulates, we studied 15 samples for petrography and microscopic analysis, and chose six representative samples (two each from felsic, intermediate and mafic types) for microprobe analysis. In each sample, 3–4 of the identified zoned sanidine crystals (which are very rare, < 1%; see description in the next section) were analyzed by electron microprobe.

Quantitative analyses were done at GZG, Universität Göttingen with a JEOL JXA-8900R microprobe. The measurements were performed with an accelerating voltage of 15 kV, beam current of 15 nA and beam diameter of 10 μm. For all the elements other than Ba and Sr, counting times were 15 s on the peak and 5 s on the background. For Ba and Sr, up to 120 s on the peak and 60 s on the background were used, but concentrations in the sanidines from the most evolved pumice and cumulates were still below the detection limit (79 and 220 ppm for Ba and Sr, respectively). Calibration standards used during the measurements were albite for Na, sanidine for K; anorthite for Si, Ca and Al; hematite for Fe; SrTiO3 for Sr and celsian for Ba. The relative standard deviation for all major oxides was < 2% and the calculated absolute error for the minor oxides was 0.003–0.03 wt%. Accumulated back scattered electron (BSE) images were obtained by superimposing ten BSE image accumulations acquired using map analysis in COMPO mode with 20 kV accelerating voltage and 20 nA beam current and with a slow scanning beam with acquisition time of 120 s per accumulation (e.g. Figs. 4, 5, 6, 7 and 8). Compositional profiles across zone boundaries and from crystal core to rim were obtained perpendicular to the zone boundaries while referring to the BSE images. To estimate the lateral spatial resolution of the BSE images, a Monte Carlo simulation was done. The simulation suggests a lateral resolution of < 500 nm with 95% of the BS electrons coming from a depth of < 450 nm with lateral diameter of ≤ 400 nm (online resource 3). The obtained quantitative data and BSE images for the sanidine phenocrysts in the phonolite and in the cumulate samples are provided in online resource 1 and 2. The specifics of the obtained results are discussed below.

Fig. 4
figure 4

Accumulated backscattered electron images of selected crystals from the three major stratigraphic units of Laacher See Tephra [a & b: LLST (sample 1002); c and d: MLST (sample 1034 and 1050); e and f: ULST (sample 1088-H and 1060–51) illustrating the intra-grain compositional zonation and textural features. a and c are from Ginibre (2000). b and d are from Ginibre et al. (2004) with permission of Oxford University Press (License Number 4858160241332)

Fig. 5
figure 5

Examples of the compositional variation in the sanidines from the pumice samples. a selected crystal from ULST sample 1088-H. The red line denotes the measurement points over which the EMP quantitative line analysis is done and compositional profiles of b Or-content, c Ba (wt%) and d Sr (wt%) are obtained. The yellow rectangle in a denotes the area over which e the grey-scale profile is taken. The type of compositional variation in a is representative of all the sanidines from MLST and ULST. f Selected crystal from LLST sample 1002 (from Ginibre 2000) with highlighted compositional difference in terms of Or-content. The Ba- and Sr-content in these samples is below detection limit of EMP ((79 and 220 ppm for Ba and Sr, respectively)

Fig. 6
figure 6

Accumulated backscattered electron images of selected crystals (1088-H sanidine 2 and 1088-H1 sanidine 2) from the hybrid samples (upper ULST). The thin (< 10 μm) Ba-rich overgrowths (light grey in color) are evident in the crystals

Fig. 7
figure 7

Accumulated backscattered electron images of selected crystals from felsic and mafic cumulates. a an example (281 san5, felsic) from the unzoned crystals which constitute majority (> 99%) of the population in all types of cumulates. b Rare example of zoned sanidine crystal (280 san2z) from felsic sample 280 showing variation in terms of Or-content (just like LLST sanidines, Fig. 5f). Ba- and Sr-contents are generally below detection limit. c Rare example of zoned sanidine crystal (294 san1z) from mafic sample 294 showing variation in terms of Ba-content and homogeneity in terms of Or-content (just like MLST and ULST sanidines, Fig. 5a)

Fig. 8
figure 8

a Accumulated BSE image of a selected sanidine grain [sanidine 2 in sample 1088-1H (a hybrid sample from upper ULST)]. The solid red rectangle encloses the area over which the grey-scale swath profile (c) for inner boundary is taken. b The zoomed-in view of a, the solid red rectangle covers the area of the swath profile (d) for grey-scale values used for diffusion modeling across the outer-most boundary. The red curves in c and d are the model curves obtained through a least square fitting (see discussion below)

Compositional zonation and interpretation

Results from compositional analysis

Sanidines from pumice

The sanidine phenocrysts in the sample from LLST, i.e. the most evolved phonolite, are characterized by 65–68 wt% SiO2, < 0.3 wt% CaO, ~ 0.15–0.47 XOr and < 100 ppm Ba and Sr. More than 60% of the crystals are patchy-zoned in orthoclase-content (Or15–Or47) (Figs. 4a, b and 5f). The zone boundaries, if identifiable, are always over printed by exsolution (between the orthoclase- and albite-rich phases). These observations make these sanidine crystals compositionally and texturally identical to the sanidine grains in the carbonatitic syenites analyzed by Schmitt et al. (2010) and Rout and Wörner (2018). Therefore, we interpret these grains to be xenocrysts coming from the overlying crystalline syenitic carapace which were incorporated into the underlying phonolite during the final destabilization of the system (probably weeks before eruption; Rout and Wörner 2018). A similar interpretation was also proposed by Bourdon et al. (1994) to explain the disturbed U-Th isochrons. The resorbed rims of these crystals suggest disequilibrium with the host magma and thus, support this interpretation.

Sanidine phenocrysts from MLST are characterized by alternating growth and resorption bands (Fig. 4c, d). They have 63–66 wt% SiO2, 0.4–1 wt% CaO and ~ 0.63–0.69 XOr, with distinctly higher Ba- and Sr-content (between 0.4 and 1.6 wt%, and 0.1 and 0.4 wt%, respectively; Fig. 5) compared to LLST sanidines. Unlike the sanidine phenocrysts from LLST, they are nearly homogenous throughout crystal in terms of major element composition. Thus, the observed zonation is only due to variation in Ba and Sr (Fig. 5). Each sanidine grain has 2–4 significant or well-identified growth bands that are characterized by a change of at least 0.5 wt% in Ba content. Apart from these growth bands, the crystals occasionally (< 5% of the population) also show oscillatory zonation with a general increase in Ba from the core followed by a decrease towards the rim (Fig. 4d). The sanidine phenocrysts from the bottom of the reservoir or ULST (Fig. 4e, f) are compositionally very similar to the ones in MLST with Ba concentration (0.5–2.4 wt%) overlapping (online resource 2). Like the sanidines in MLST, the overall major element composition show minor variations compared to larger ranges in Ba and Sr contents (e.g. Fig. 5). An important observation is that 60–70% of the sanidine phenocrysts from the upper ULST samples, which are the hybrid between the phonolite and the basanite, have an additional very thin (2–10 μm wide) Ba-rich (at least 1.8–2.5 wt%) overgrowth (Fig. 6) with 60–70% change in Ba compared to the preceding zone. This kind of overgrowth is never observed in the sanidine phenocrysts from MLST and non-hybrid ULST samples.

Sanidines from the cumulates

The sanidines in the cumulates are compositionally equivalent to their stratigraphic pumice counter-parts. The composition of the sanidines in the felsic cumulates are similar to those in the LLST and syenites, characterized by ~ 0.2–0.46 XOr and extremely low (< 200 ppm) Ba and Sr concentration (Fig. 7, online resource 2). Less than 1% of the crystals are zoned. The zonation is mainly patchy with minute exsolution texture (see examples in online resource 1). Sanidines in the intermediate and mafic cumulates contain ~ 0.6–0.7 XOr, 0.1–0.5 wt% Sr and 0.1–1.8 wt% Ba, similar to MLST and ULST (Fig. 7; more details in online resource 2). Few of the sanidines in all three types of samples, share ragged outlines with surrounding glass and have holes filled with glass suggesting possible resorption and melting.

In the following sections, we discuss the causes of zoning patterns and their correlation with respect to differences between different samples, and link these to the crystal residence times obtained from Ba-diffusion modelling.

Zonation patterns

The zonation patterns in the sanidine phenocrysts from MLST and ULST samples can be broadly divided into two categories: multiple growth bands (up to 200 μm in width; Fig. 4c, e, f) and occasional oscillatory zoning (up to only 20 μm width, Fig. 4d). Growth bands are defined by different Ba-contents with or without evident resorption surface in between. They often show a week saw-tooth pattern where the Ba-content starts to decrease gradually within a band after a relatively sharp increase at the preceding resorption interface. This gradual decrease continues until the next resorption event, which is again followed by a sharp increase in Ba-content. To interpret the texture and cause of different kinds of zonation, several questions need to be answered:

What causes the growth bands in the crystals?

The growth bands may have been caused by one or combination of the following: (1) local kinetic control (2) external factors e.g. melt composition, temperature, pressure and water-content (Shore and Fowler 1996; Ginibre et al. 2007).

Local kinetic control involves the development of a diffusion boundary layer that forms over time at the crystal–melt interface. However, Ginibre et al. (2004) observed correlations between all minor and trace elements, especially Ca with Sr and Ba, i.e. both compatible and incompatible elements in sanidine (KBa and KSr > 5, and KCa < 1; Stix and Gorton 1990; Ginibre et al. 2004). A diffusion boundary layer at the crystal-melt interface cannot explain such correlations. This indicates that the factors controlling these zonations are external e.g. temperature, melt composition, water-content and/or pressure. However, given the small range of pressure inside the Laacher See reservoir (a gradient of 45 MPa between LLST and ULST; Harms et al. 2004), the effect of pressure on the partitioning of Ba would be negligible (Guo and Green 1989). That leaves only temperature, melt composition and water-content as the plausible options. Similarly, the effects of temperature and water-content will be small on the partitioning of trace elements (Ba, Sr) between sanidine and melt, although such effects cannot be completely excluded. For crystals growing inside a magma reservoir, changes in melt composition (and subordinate temperature and water-content) can occur in several ways including: (1) by intrusion of a magma of different composition and temperature, and (2) by movement of the crystals into different parts of a chemically and thermally zoned reservoir.

Compositional effects on growing crystals by intrusion of a compositionally different magma followed by mixing, however, would change the bulk major element composition of the melt affecting the composition of parts of the crystals that grow post mixing. In this case, we would see changes in the major element composition of the sanidines, which is not the case. Thus, any possible role of recharge would have to be limited to changes in the temperature regime of magma storage. Such a scenario is supported by various reservoir models (e.g. Wolff et al. 2015), where a fresh hot mafic recharge underplates the reservoir while significantly affecting the system thermally with negligible or very limited chemical mixing. The thermal effects (i.e. heating) from such a process would lead to resorption of the crystals enriching the surrounding melt in compatible elements like Ba and Sr, which, in return, will enrich parts of the crystals that grow from this melt when temperature starts to drop afterwards. In fact, the strong resorption events recorded in the sanidine crystals from MLST and ULST (e.g. Figs. 4c, e, f and 3b), that reflect drastic temperature variations of > 100 °C (based on crystallization temperatures calculated from ternary composition of sanidine by Ginibre et al. 2004), are most likely to be effected by hotter recharges. A recent study (Tomlinson et al. 2020) also proposes periodic mafic recharges in Laacher See based on the major and trace element data of LST glasses suggesting periodic addition of basanite magma that significantly contributes to the heat budget of the system. A similar pre-eruptive basanite recharge is also extensively documented in the hybrid samples from upper ULST (Wörner and Wright 1984). Thus, hot mafic recharge is certainly a plausible reason behind the observed zonation. However, although hot mafic recharges can explain the growth of Ba-rich zones preceded by resorption, they cannot explain the growth of distinct Ba-poor zones with sharp boundaries (e.g. Fig. 4c, f).

Since major element compositions within each crystal remain almost constant and only trace element composition vary significantly and over several zones (online resource 2; Ginibre et al. 2004), we argue that variable composition of a phonolitic melt, in paticular different in Ba and Sr contents, plays a significant role in these zonations. This is supported by significantly different Ba and Sr concentration in different parts of the zoned Laacher See phonolite melt (Fig. 2; Wörner and Schmincke 1984a).

The movement of the crystals to different parts of reservoir can happen through convection and in the case of Laacher See reservoir, convection inside the liquid-rich core can definitely be a possibility provided the crystallinity is low (below the rheological lock-up i.e. 40–60 vol% crystals, Bachmann and Bergantz 2008; Cooper and Kent 2014, Rubin et al. 2017). The vol% of crystals within MLST samples is between 1 and 15%. With convection, settling, and recapture, crystals could have been repeatedly transferred down into hotter and drier (high Ba, high Sr) phonolite melt and back to colder and water-rich parts of the zoned magma reservoir (low Ba, low Sr). This interpertation is based on the observed compositional zonation in which Ba- and Sr-rich, hotter and drier magmas reside at lower levels in the magma chamber, overlain by Ba- and Sr- poor, volatile-enriched colder phonolite magma (Fig. 2; Wörner and Schmincke 1984a; Berndt et al. 2001; Ginibre et al. 2004). When transported to a different part of the zoned magma column with different melt composition, the crystals may first undergo dissolution and, upon cooling, will start to grow again from a melt rich in Ba and Sr. When these crystals are transported back into to a colder and more evolved part of the zoned magma reservoir, they continue to grow from a Ba-poor melt, but without resorption. The compositional zoning shown in Fig. 4f is an example of such a case. Such growth bands are indeed observed in most sanidines from MSLT and also from ULST where the crystallinity of the samples is between 30 and 50%. At 40–50% crystallinity (i.e. the threshold of rheological lock-up), the melt would start to behave as a rigid highly viscous crystal mush making convection and crystal transfer not a feasible process (Bachmann and Bergantz 2008; Cooper and Kent 2014; Rubin et al. 2017). However, petrographical textures clearly show that accumulation and cumulate compaction must have followed crystallization, in which case the growth of crystals would actually occur, and zonations would have formed, initially at lower crystallinities also for ULST sanidines. However, the preserved overall chemical zonation within the overall deposit (including cumulates), suggests against any large reservoir-scale convection. Any possible convection and crystal transfer between melts of different composition would have to be small-scaled and localized i.e. limited to neighboring magma batches.

Thus, there is evidence for periodic recharges and small scale convection that are both needed to explain the complex intra-crystal and intra-reservoir zonation in the phenocrysts of the Laacher See magma system. We, therefore, propose a combination of the two processes i.e. repeated recharges of hotter basanite magmas that provide heat and set off small/local convective plumes to facilitate crystal transfer within the zoned magma reservoir. The recharges would cause temperature variations that result in resorption of sanidine crystals, followed by Ba-rich growth whereas alternatively or additionally, crystal transfer between melts of distinct trace element composition would result in distinct Ba-poor and Ba-rich zones with sharp (and non- or minorly-resorbed) boundaries between growth bands.

The gradual decrease in Ba within each Ba-growth band needs a different explanation. Since Ba is compatible in sanidine, the host phonolite melt will be continuously depleted in Ba during sanidine crystallization. This depletion in the host melt will reflect in the crystal growth leading to the observed intra-zone gradual decrease in Ba. Whether this depletion only affects the melt boundary layer near the crystal due to sluggish Ba diffusion (e.g. at 750–850 °C, diffusivity of Ba in phonolitic melt is estimated to be of the order of ~ 10–13–10–15 m2/sec, when extrapolated from Zhang et al. 2010) or affects the entire surrounding melt body is difficult to establish. The gradual depletion in Ba in the growing sanidine would only be interrupted when the melt composition around the crystal changes e.g. through dissolution or movement of the crystal into a different melt (see above) causing high-Ba or low-Ba zones (Figs. 4c, e, f, 5 and 6).

What causes the fine scale oscillatory zonation?

This kind of zonation is only represented either in the mantle or rim of < 5% of the crystals and is associated with very small Ba-changes (~ 0.2–0.5 wt%) at higher spatial frequencies (< 20 μm). These relatively fine-scale partially euhedral oscillatory zones trace the crystal shapes and probably reflect incremental diffusion-controlled growth (Haase et al. 1980, Pearce and Kolisnik 1990, Singer et al. 1995). This type of growth implies steady-state cooling, low undercooling conditions, and negligible melt convection within the diffusion boundary layer. The growth may also have been associated with slight dissolution as the zone boundaries, being typically euhedral, often also have slightly rounded corners and sporadically low amplitude undulating resorption surfaces. Such zonation has also been observed in the plagioclases from Laacher See (Ginibre et al. 2002a). Shore and Fowler (1996) and Ginibre et al. (2002a, b; 2007) have summarized some of the existing models that are used to explain such zonation. This process might not be fundamentally different from the process mentioned above explaining the gradual Ba-decrease within each growth zone, but the spatial and temporal scales over which both the processes operated, might have been very different. However, since their formation does not reflect, and does not allow for estimating diffusion (or crystal residence) times at large, we do not model or analyze this type of zonation further.

What is the reason behind the thin high-Ba overgrowth in the hybrid samples?

The thin (2–10 μm) overgrowths (Fig. 6) in the sanidine crystals in the samples from the very top of the tephra (ULST), which are also the hybrid between the phonolite and the intruding basanite, contain the maximum Ba-content (at least 1.8–2.5 wt%) compared to other parts of the crystals. This overgrowth is often preceded by strong resorption and a sharp Ba jump (Fig. 8b, d, and and online resource 1). Other crystals from these hybrid samples also contain strongly resorbed boundaries without any overgrowth. These observations suggest that this overgrowth must have preceded by a very strong, short and sharp heating event just before the eruption, which we interpret to be the arrival of the basanite that mixed and mingled with the resident phonolite.

It is, however, unclear whether the overgrowths preceded or succeeded the actual mixing process. If preceded, the heat (for resorption) must have come only from the underplating of the basanite (before mixing) and the overgrowths would grow just from the surrounding Ba-enriched melt (just like other Ba-rich zones in the crystal interiors; see above). However, after heating and resorption, it seems unlikely that the temperature would drop significantly (for regrowth to begin) between underplating and mixing. Several observations, e.g. non-resorbed exteriors in some of the overgrown crystals (Fig. 6a; online resource 1) and higher Ba-content in the overgrown parts compared to rest of the crystals suggest a post-mixing origin of the overgrowth. So, could these sanidines crystallize a high-Ba rim after basanite-phonolite mixing? To evaluate this question, (1) the Ba-content and major element composition and (2) the solidus-relations for sanidine of the hybrid host's matrix need to be considered. Wörner et al. (1983) reported the Ba-content in the matrix of the hybrid sample 1099 to be at ~ 1300 ppm which, based on Ba-partition coefficients (Stix and Gorton 1990), is in equilibrium with the overgrowth. However, such high Ba contents in the hybrid matrix cannot be simply explained by phonolite-basanite mixing considering the Ba contents of the pure endmembers. Ba in the matrix of ULST samples 1060 and 1088 (just above the basanite-phonolite mixing zone in the magma chamber) is 595–599 ppm (Wörner et al. 1983). The basanites in the East Eifel Volcanic Field have typical Ba contents of 942–1026 ppm e.g. Eppelsberg volcano (unpublished data) and Nickenicher Sattel (Wörner et al 1983). Therefore, the Ba-content of a hybrid matrix would be between the ~ 600 and ~ 1000 ppm. Thus, the basanite-phonolite mix alone cannot explain the ~ 1300 ppm Ba in the matrix and the resulting Ba jump in the outermost rims. Therefore, an additional Ba-enriched component such as of melted sanidines is required to have been incorporated into the melt (as suggested by Wolff et al. 2020 for certain silicic systems) either during compositional mixing or even before mixing during an earlier stage (underplating) of the recharge. Either way, the hybrid melt most likely contained components from both the basanite and the melted sanidines which prompted the Ba jump in the outermost rims that grew from this hybrid melt. Wörner and Wright (1984) also reported similar overgrowth of phonolitic clinopyroxene and amphiboles from the hybrid matrix that formed after the mixing event. Other hybrid samples show only minor modal sanidines which are often of completely rounded and strongly resorbed shape. In such cases, sanidines were dissolved in the hybrid but never started to grow during later cooling. Therefore, depending on how far such Ba-enriched hybrid melts cool before eruption, they may or may not crystallize sanidine on the previously melted crystals.

There are two major conclusions from these considerations that are highly relevant for the diffusion modeling and the significance of the diffusion timescales: (1) Sanidine resorption and Ba-rich overgrowth in the hybrid magma are direct consequences of mixing, and constraining the timing of this overgrowth should constrain the maximum time between the onset of mixing-induced regrowth and eruption, and (2) the temperature of ~ 1050–1170 °C estimated for the hybrid magma (Wörner and Wright 1984; Sundermeyer et al. 2020) is not relevant for the Ba-diffusion across the core-rim overgrowth because the sanidines can only grow from the hybrid magma after cooling to its liquidus temperature (890 ± 10 °C, see below).

How can the zonation patterns in the cumulates be interpreted?

In cumulates, only very few sanidines (< 1%) are zoned. The partly patchy core-rim zonation of orthoclase content (i.e. Na- and K-content) in sanidines from a felsic cumulate (Fig. 7b and online resource 1) is most probably the result of carbonate–silicate unmixing as interpreted by Rout and Wörner (2018) for similar zonation patterns in sanidines from carbonate-bearing syenites. However, unlike the syenites, the zonation with Or-rich and Ab-rich phases are generally not overprinted by exsolution (Fig. 7b). This is probably due to faster cooling to temperatures below the solvus (712–725 °C; Brown and Parsons 1989) for the felsic cumulates. The zoned crystals in the intermediate and mafic cumulates (Fig. 7c) are identical to the ones in the pumice samples and hence must have formed by the similar mechanism as discussed above. However, the majority (> 99%) of the cumulate crystals are completely unzoned (Fig. 7a). These crystals must have either been homogenized completely by diffusion or grew unzoned from the same phonolite melt without Ba-depletion and then accumulated almost immediately before undergoing further growth or any thermal changes. In any case, the majority of the cumulate crystals must have had a different growth and/or storage history than the phenocrysts observed in the pumice samples.

Ba-diffusion modeling

Each sanidine crystal from the pumice samples has 2–4 distinct growth zones with at least 40% change in Ba content. Modeling of these profiles for diffusive homogenization should constrain the pre-eruptive residence times of these crystals. However the eligibility of these Ba-gradients for diffusion modeling needs to checked first in terms of necessary boundary conditions.

Boundary conditions for diffusion modeling

Starting profiles

All major zone boundaries relate to dissolution interfaces as shown by their curved shapes. This petrographic observation suggests that the compositional gradients (up to 50 μm in length) were resorption surfaces formed by dissolution followed by regrowth (Figs. 6 and 8). The compositional boundary formed by regrowth should be a sharp step-like profile. But, given the high compatibility of Ba in sanidine, the possibility arises for in-diffusion of Ba from Ba-rich surrounding melt into the crystal during the dissolution itself (Chamberlain et al. 2014). We modeled the extent of Ba in-diffusion that would happen until the Ba-rich surrounding melt layer persists. As it turns out, due to the high diffusivity of Ba in melt compared to sanidine (by eight orders of magnitude; Zhang et al. 2010), the Ba-rich melt layer will diffuse away before the in-diffusion proceeds significantly into the crystal to affect the initial profile (see online resource 3 for details). This implies that the initial compositional boundary formed as soon as the regrowth starts was a rather sharp and close to an ideal step profile. This assumption is also supported by the quantitative analysis of two boundaries (Fig. 5 and online resorce 2), that shows ~ 10 times wider gradients of Sr compared to Ba (which has a ~ 100 times slower diffusivity; Cherniak 2010). This satisfies one of the basic criteria of applying diffusion modeling, and makes these compositional gradients across resorption boundaries ideal for diffusion modeling.

Starting concentration

The next basic requirement to apply diffusion modeling lies in constraining the initial concentrations, which can be done using the compositional plateaus on both sides of a diffusion profile. Similar to the classification by Chamberlain et al. (2014), we divide our diffusion profiles (Fig. 9) into two categories: type 1: Profiles that fit perfectly to a one-dimensional diffusion profile between two semi-infinite bodies and thus can be efficiently modeled for the same, and type 2: Profiles that do not follow an ideal diffusion profile because of, e.g., lack of compositional plateaus on either side. Unlike type 1, the initial concentrations in type 2 profiles are not well-constrained. However, due to the small number of suitable, well-zoned crystals in LST, especially MLST, we cannot completely exclude these “imperfect” diffusion profiles from our analysis. Therefore, we also modeled these profiles for a one-dimensional diffusion between semi-infinite bodies using a multi-parameter curve fitting algorithm that works for unknown initial concentrations. However, as Chamberlain et al. (2014) suggested, diffusion modelling of type 2 profiles would most probably lead to an overestimation of diffusion time. Thus, the timescale estimates from type 2 profiles need to be considered as the maximum possible timescales.

Fig. 9
figure 9

Two types of diffusion profiles modelled in this study: a type 1 (with visible compositional plateaus on both sides) and b type 2 (without evident compositional plateaus on either side)

Most compositional gradients in our samples are less than 40–50 μm wide (online resource 1 and 2) and could be better analyzed using grey scale profiles from BSE images to achieve better spatial resolution (< 500 nm, see analytical methods for details). In Laacher See sanidines, the intensity of the BSEs only depends on Ba-content in all the samples analyzed here for diffusion chronometry because the major element contents are fairly constant. Therefore, the grey scale values in the BSE images correlate directly and only to Ba-contents. This has been tested by comparing X-ray intensities from Ba element maps with grey scales from BSE images over the same crystal transects (Fig. 10). Such correlation for sanidines was also previously confirmed by Ginibre et al. (2004) and Iovine et al. (2017). Consequently, grey-scale values can be safely used as proxy for Ba-content for modeling of Ba-diffusion.

Fig. 10
figure 10

a Correlation between Ba-content (cps) and grey-scale values from a selected crystal (1099 sanidine 4). The values are obtained using swath profiles (yellow rectangles) in the BSE image (b) and Ba-intensity X-ray map (c). b and c are from Ginibre et al. (2004) with permission of Oxford University Press (License Number 4858160241332)

Constraints on critical parameters

Temperature

Harms et al. (2004), through a series of experiments, estimated the temperature of crystallization of phenocrysts in the phonolite magma to be 750–760 °C at 115–145 MPa pressure for the evolved phonolite melt composition i.e. top of the magma chamber or LLST. Their experiments also constrain the sanidine-in temperature for LLST at ~ 750–775 °C for pressure between 75 and 125 MPa. Berndt et al. (2001), who used the composition equivalent to the samples used in this study, confirmed similar sanidine crystallization temperatures for MLST and ULST phonolites at 760–780 °C and 840–860 °C, respectively. Wörner (1982) also determined the sanidine-plagioclase crystallization temperature within ULST to be < 890 °C. Considering these constraints, we used 750 °C for LLST, 760 °C for MLST and 840 °C for ULST samples. We have made sure to use the lower limits for the temperatures in order to obtain maximum diffusion times. An uncertainty of ± 10 °C is assumed on the temperatures based on the experimental uncertainty and comparison between different data sets.

Temperature for hybrids: The temperature for the hybrid (upper ULST) needs different constraints as opposed to the rest of ULST. For the pre-mixing temperature of the hybrid’s phonolitic parent, we use the maximum temperature documented for ULST i.e. 870 ± 10 °C (see above; Wörner 1982; Wörner and Schmincke 1984a, b; Berndt et al. 2001). This temperature is used for the diffusion boundaries in the cores and mantles of the hybrid sanidines, which are interpreted to have formed pre-mixing in crystal-rich ULST-type mafic phonolite. However, this temperature may not apply to the outermost rims as they are interpreted above to have formed post mixing with the basanite. The temperature of the basanite is constrained at ~ 1050–1172 °C (Wörner and Wright 1984; Sundermeyer et al. 2020). Assuming that the hybrids contained ~ 50% basanite by mass (Wörner and Writght 1984) we get 960–1021 °C (i. e. the average of 870 and 1050–1172 °C) as the post-mixing temperature. However, sanidine would not grow at such high temperature and overgrowth must have started after the basanite/hybrid lost some of its heat to the overlying reservoir which is observed in resorption of crystals in LLST and MLST phonolite. Therefore, for the post-mixing outermost rims, we use, as a minimum, the highest liquidus temperature documented for upper ULST sanidines i.e. 890 °C (Wörner 1982; Berndt et al. 2001) with an uncertainty of ± 10 °C. Using this minimum temperature of course implies that the estimated diffusion times will be maximum.

Diffusivity

Diffusion coefficients are calculated using the Arrhenius equation between the diffusion coefficient, diffusivity parameters (activation energy and pre-exponential factor) and temperature as follows.

$$D={D}_{0}{e}^{\frac{{-E}_{A}}{RT}}$$
(1)

where D is diffusion coefficient (m2/sec), EA (kJ/mol) is the activation energy, R [8.3145 J/(mol.K)] is the universal gas constant, D0 (m2/sec) is the pre-exponential factor and corresponds to the value of D (m2/sec) at infinite temperature and T is the temperature in K. We use diffusivity data for Ba given by Cherniak (2002). Cherniak (2002) experimentally determined the diffusivity data for Ba in natural sanidine under dry 1 atm conditions. No compositional dependence of Ba diffusion on sanidine major element composition was reported. Also, no significant diffusion anisotropy was documented as diffusion coefficients for the direction normal to (001) was comparable to the direction normal to (010). Thus, diffusivity data only for the direction normal to (001) was provided with activation energy of 455 ± 20 kJ/mol and pre-exponential factor of 2.9 × 10−1 m2/sec, which are used in this study. These parameters are valid for a temperature range from 775 °C to 1124 °C, showing an Arrhenius relationship (Eq. 1) between diffusion coefficient and temperature. The temperatures used in this study lie between 750 and 870 °C. We assume the Arrhenius relationship to continue down until 750 °C and thus extrapolate to obtain the diffusion coefficient at 750 and 760 °C. The dependence of Ba diffusion on pressure and fO2 has not been detected (Cherniak 2002, 2010) and thus, pressure and fO2 are not constrained or used in our diffusion modeling.

Diffusion modeling

Each zone boundary is then modeled using a least square curve fitting (using MATLAB and ORIGIN) of the diffusion profiles. We use the analytical solution for diffusion between semi-infinite bodies (Crank 1975) i.e.

$$C\left(x,t\right)=\frac{({C}_{high}+{C}_{low})}{2}-\frac{\left({C}_{high}-{C}_{low}\right)}{2}\times erf\left(\frac{(x-{x}_{0})}{\sqrt{4Dt}}\right)$$
(2)

where C(x,t) is the Ba concentration at position x after time t (sec). Clow and Chigh are the initial Ba concentrations (observed as the plateaus on both sides of the gradient), D is the diffusion coefficient (m2/sec) and x0 (m) is the position of the diffusive interface. The curve fitting of a profile gives the corresponding diffusion time (t). Figure 8 shows an example of such diffusion modeling applied to the grey-scale profiles at a crystal core and rim. To constrain the minimum resolvable timescale, we modelled a profile obtained across a crack (thus, a known sharp boundry) for Ba-diffusion (see online resource 3), which gives apparent timescales of ~ 90 days at 890 °C, 7 months at 870 °C, 2 years at 840 °C and 34 years at 780 °C. Results from the diffusion modeling of all the samples are presented in Table 1. The errors on the timescales (1σ) comprise of the standard error (1σ) on curve fitting (which also includes the uncertainty in position i.e. the spatial resolution of BSE images) and the uncertainty in temperature (refer to online resource 3 for details on error calculation).

Table 1 Estimates of time scales obtained for diffusion across zone boundaries in the sanidines from MLST and ULST

Results and implications for magma residence and mixing times

Phonolite pumice

Calculated Ba-diffusion times do not vary significantly between sanidines from MLST and ULST phonolite (Table 1). Total residence times as calculated from the inner-most boundaries near the crystal cores range from 4.3–23.5 ky for MLST and 1.5–20.3 ky for ULST. The range of up to tens of thousands of years of crystal storage and growth at 760–840 °C is identical to the time estimated from U-series dating on the syenitic crystal-rich carapace around the phonolite magma body (Schmitt et al. 2010) and the crystal age of cumulates (Bourdon et al. 1994). Applying diffusion modeling to the outer zones and the rims of the crystals gives much shorter diffusion times of 1.4–9.3 ky (Table 1). This suggests that the onset of sanidine crystallization dates back to ~ 24 ky prior to eruption (Table 1) with crystallization continuing until < 1.4 ky before eruption, i.e. essentially up to the time of eruption. The distribution of diffusion times (Fig. 11) also suggests that repeated heating and dissolution events that formed the observed Ba jumps and zone boundaries occurred throughout the lifespan of the sanidine crystals. This recurrence of heating events (or recharges) is observed in both MLST and ULST, although ULST seems to document more of the younger events compared to MLST. However, due to the scarcity of crystals in MLST, these events may not be recorded or documented in our analysis.

Fig. 11
figure 11

a Frequency histogram showing the diffusion times obtained from Ba-diffusion modeling of 39 zone boundaries in 27 crystals (from both MLST and ULST). The distribution of the diffusion times indicates the occurance of heating events throughout the crystal residence. The histogram excludes the durations from the outer-most boundaries. Frequency histograms of cyrstals from only MLST (b) and only ULST (c) are also included for comparison

While most of the younger crystals (residence time < 10 ky) suggest time difference of < 3–6 ky between the core and the rim, the older crystals suggest a difference of ~ 12–24 ky (Table 1). This could easily be the result of intermittent strong dissolution of older cystals that dissolved their relatively younger exteriors preserving only the older cores. After dissolution, the rims which then are much younger compared to the cores would form leading to a large age gap between the cores and rims. Such strong intermittent dissolution is also supported by the observation that only six crystals out of 27 have residence times longer than 15 ky suggesting that very few of the older crystals survived all the heating and dissolution events while most of them could not. Therefore, most (up to 75%) of the erupted crystals are relatively young (< 15 ky before eruption). With these results including the prolonged crystallization history, and considering the low solidus temperature of Laacher See phonolite (Harms et al., 2004; Berndt et al. 2001), we conclude that MLST and ULST, representing the main body of the Laacher See reservoir, were kept at relatively high temperatures and sigificantly above the phonolite solidus for very long durations of at least 20 ky. Thus, the main part, namely the liquid core of the Laacher See phonolite magma reservoir, which was enclosed by a crystal-rich carapace, represents an example for “warm storage” (Barboni et al. 2016).

The thin outermost Ba-rich rims of the sanidines in the hybrid samples tell a completely different story because Ba-diffusion modeling gives a very short timescale of only ~ 1.5–3.4 years (Fig. 8 and Table 1). This implies that the last and critical recharge and heating event happened just one to three years (maximum) prior to eruption. The lowermost crystal-rich layer at the base of the phonolite magma reservoir (i.e. upper ULST) was evidently affected by this event (both chemically and thermally). Although we never observed such late Ba-rich overgrowth in samples from MLST, the occurrence of strongly resorbed crystals most probably is the result of the thermal effects of this recharge. Thus, the effects of this pre-eruptive heating event were probably not limited to the lowermost part of ULST. Rather, heating at the base of the magma column, overturn and initiation of degassing may also have been the result of this event, which eventually triggered the eruption.

Cumulates

No reliable temperature constraints exist for the cumulates. We modeled the rare diffusion profiles assuming similar temperature as for the corresponding pumice samples of similar major element composition of the contained crystals. For a temperature of 760–840 °C, the calculated durations were at only ~ 80–215 years (online resource 2). This represents the maximum diffusion time that these crystals might have spent at these high temperatures, but they are considerably shorter compared to diffusion times obtained from similar phenocrysts in MSLT and ULST pumice. This suggests that the cumulates might be very young. However, this is unlikely since Bourdon et al. (1994) dated a few of these cumulates (using U-series radiometric dating) and found them to be as old as 19–28 ka, implying pre-eruption storage times of at least 6–15 ky. Alternatively, the cumulate crystals may have rapidly cooled to temperatures below the diffusion threshold where diffusion would be extremely slow, which is a more plausible explanation for these very short diffusion times in the cumulates. We estimated the apparent diffusion threshold for Ba (below which diffusion is negligible) to be below ~ 700 °C (online resource 4) at which the cumulates must have spent most of their storage time. Therefore, the cumulates document storage at a comparatively lower temperature than the phonolite.

Crystal storage conditions prior to eruption

Figure 12 summarizes the thermal and temporal constraints obtained from three distinct kinds of eruption products from the Laacher See volcano: phenocrysts in phonolite magma, cumulates and carbonatitic-syenites (Rout and Wörner 2018). For the phonolite magma, crystal residence and storage times up to ~ 24 ky at temperatures of 760–840 °C suggests long-term storage and growth well above the phonolitic minimum. The crystallinity at these temperatures will be significantly below 50%, indicating that this part of the magma reservoir (MLST and LLST) was stored for thousands of years as a hot liquid-dominated body. It also implies that the magma reservoir remained in an eruptible state for thousands of years while it was replenished (and heated) from below by mafic magmas. Such frequent recharges apparently did not trigger an eruption during this long time of storage. However, a last recharge, that resulted in mingling and mixing, did trigger the eruption eventually only shortly (maximum 1.5–3 yrs) after the intrusion of mantle-derived basanite magma into the base of the resident phonolite magma reservoir. This timescale would drop to weeks to months (minimum) if a higher temperature (e.g. > 900 °C) is considered.

Fig. 12
figure 12

Contrasting zonation types in sanidine and results of diffusion modeling in syenite, phonolite pumice and cumulates from the Laacher See Tephra

Crystals in the cumulates are generally old (formed 6–15 ky before eruption; Bourdon et al. 1994), but give very short diffusion time scales. This is interpreted to imply that the cumulates were cooled below the diffusion threshold early in their history and were stored as a rigid crystalline mush before being heated and reactivated by basanite recharge. This interpretation is similar to that for the carbonatitic syenites, also crystalline ejecta, which are also proposed to have been rapidly cooled to, and stored at, temperatures as low as 630–670 °C (Rout and Wörner 2018).

We suggest (as per the model by Rout and Wörner 2018) that after crossing the threshold for rheological lock-up, the cumulate crystals underwent conductive cooling, where the temperature rapidly dropped to temperatures below diffusion threshold within a period of tens to a few hundreds of years. The rocks continue to be stored at these low temperatures as the crystalline carapaces around the liquid-rich magma body until eruption brings them to the surface. However, the presence of interstitial phonolitic melt in these rocks suggests that these rocks were never cooled below the solidus and must have stored at as low as near-solidus temperatures. Such a largely crystalline mush surrounding a liquid-rich reservoir at the core has been proposed by several reservoir models (e.g. Hildreth and Wilson 2007; Annen 2009; Gelman et al. 2013; Caricchi et al. 2014; Annen et al. 2015).

The storage conditions at the bottom of ULST are slightly different. With crystal content reaching ~ 30–50%, the mafic phonolite at the base of the phonolite magma reservoir was at or close to the threshold for rheological lock-up. However, frequent recharges, kept this part of the system also thermally active and maintained at temperatures of 840–870 °C. In either case, the time between the final recharge and the eruption was short, between a maximum of a few years and a minimum of a few weeks.

We summarize the history for the Laacher See reservoir in Fig. 13 in terms of the crystal residence and temperature through time. The three different curves in Fig. 13 denote the three major parts of the reservoir: the ULST (red), MLST and LLST (green) and the cumulates (purple). The model begins at the onset of crystallization and ends at eruption. Crystal-poor ULST and MLST (with < 30–50% and < 15% crystal-content respectively) phonolite magma is stored at temperatures of 840–870 °C and 760–780 °C, respectively. Storage was interrupted by episodes of heating as a result of mafic recharges that caused the intermittent dissolution and growth of the sanidine crystals. Such recharges also maintained the system at relatively high temperatures. However, recharge is not documented in compositional mixing with the phonolite and thus is limited to thermal effects only. Only at the very base of the magma reservoir, mafic and crystal-rich ULST phonolite experienced a distinct heating event just ~ 1.5–3 years before eruption which also resulted in compositional mingling and mixing between basanite and phonolite. Cumulate rocks, both carbonatitic syenites and mafic to felsic cumulates from phonolite melts, were stored as a rigid crystalline mush at < 700 °C. After crystallization, the temperature of these rocks dropped to temperatures below the diffusion threshold of Ba in sanidine within tens of years. Whether or not these cumulates have experienced previous heating events is unclear, petrographical evidence, e.g. in the form of cumulate melting does not exist. However, Rout and Wörner (2018) concluded that at least the syenites experienced one short-lived heating event that led to the formation of the exsolution textures in sanidines from those samples.

Fig. 13
figure 13

Schematic diagram (modified after Rout and Wörner 2018) summarizing the thermal history [temperature (T) versus time (t)] of the Laacher See reservoir

Mafic recharge, no mush rejuvenation and CO2 as eruption trigger

Mafic recharge, mixing and mingling between phonolite and basanite are documented in hybrid magmas produced at the very end of the eruption and this may indicate that the eruption of the resident “warm” phonolite magma may have been triggered by this (last) recharge event (Sparks et al. 1977). Partial melting of crystals in the hybrids are also inferred to have been caused by this recharge (see above). Similar recharge-induced heating might have also partially melted and remobilized older cumulates (aka "mush rejuvenation"; Bachmann and Bergantz 2008; Girard and Stix 2009; Cooper and Kent 2014; Wolff et al. 2015) to contribute to the eruptible magma. Melting sanidine-rich mushes would produce phonolitic melts that are in fact observed as interstitial glass in many cumulates. However, in case of Laacher See syenites and cumulates, the interstitial glass is in equilibrium with faceted crystals (i.e. melting is excluded), and the composition of the interstitial glass is highly fractionated, i.e. low in Ba and Sr (Tait et al. 1989; Liebsch 1996). This is the opposite of what would be expected from melting of cumulate feldspars, and argues strongly against any late-stage re-melting of cumulates as a mechanism to produce significant volume of largely crystal-free, eruptible phonolite magma. However, the basanite recharge heated the lower portion of the phonolite, which may have triggered overturn, degassing and eruption. Recently, Caricchi et al. (2018) suggested a mechanism by which a mafic recharge would cause CO2 flushing from the basanite into the phonolite that could cause massive degassing of H2O and subsequent destabilisation of the magma reservoir. A similar CO2-driven explosive mechanism is proposed by Weidendorfer et al. (2019) for carbonatite-bearing systems. There is ample evidence (carbonatites, pure CO2 fluid inclusions in syenite minerals) of abundant CO2 in the Laacher See magma system and therefore, CO2-flushing could well be the final trigger for the eruption, irrespective of any additional physical and thermal effects of the recharge event.

Summary and conclusions

Different erupted products of the Laacher See magma reservoir are studied to extract their individual and the complete thermal histories of the magma system. The mostly crystalline carapace of the Laacher See magma system spent much of its ~ 24 ky of storage time (maximum) at relatively low temperatures at near-solidus conditions while the largely molten core remained hot and well above the phonolite solidus for similar residence times. Although, in terms of storage conditions, the crystalline carapace shows similarity to “cold” storage (by definition; Cooper and Kent 2014), the lack of evidence for mush melting and the evident “warm” storage conditions of phonolite (which dominates the eruption volume) make the largely molten core of Laacher See reservoir a clear case of “warm” storage, i.e. magma storage well above its solidus. This reinforces the idea that magmas and their contained crystals can be stored at eruptible conditions for thousands of years in silicic magma reservoirs (Morgan and Blake 2006; Barboni et al. 2016), even in alkaline systems (e.g. Scheibner et al. 2008; Arienzo et al. 2011). While the crystal-rich carapace remained relatively cool, the temperature of the molten core was maintained well-above the solidus by regular heating and recharge events. These heating events caused dissolution of the long-stored crystals, but apparently could not initiate an eruption.

The last basanite recharge event at Laacher See was either sufficiently forceful or encountered an eventually volatile over-saturated phonolite magma and, therefore, resulted in eruption. This final eruption trigger was either by forceful injection, mixing, mingling, and overturn or by CO2-flushing and massive H2O degassing (Caricchi et al. 2018). This last recharge happened—at most—only a very few years before the eruption. Such short timescales between recharge and eruption are maximum estimates and fully in accordance with diffusion modelling of zoned olivine crystals contained in late Laacher See hybrid lavas by Sundermeyer et al. (2020) who conclude that mixing and eruption occurred "within a few to a maximum of 14 months after the intrusion of a new mafic basanite magma".