Chemical interdiffusion between Na-series tephritic and phonolitic melts with different H₂O content, temperature, and oxygen fugacity values

The diffusive exchange of major elements in Na-series tephrite–phonolite diffusion couples with compositions relevant to the Canary Islands magmatism was determined at 300 MPa and variable H₂O concentrations (0.3 wt % to 3.3 wt %), temperatures (1150 to 1300 °C), and fO₂ levels (NNO−1.5 to NNO+1.7). Composition-dependent effective binary diffusion coefficients were determined from concentration–distance profiles. Results show a wide range of diffusivities for different cations, consistently following the sequence Na ≫ Al ≫ K ≥ Mg = Fe = Ca > Si > Ti, with a mild diffusivity contrast (0.2–0.8 log units) between tephritic and phonolitic melts. Na is the fastest component, with diffusivities falling ∼1.0 log units above those of Si for any given condition. An anomalously fast Al diffusion is observed, with D_Al falling ∼0.4 log units above Si and ∼0.6 log units below Na, suggesting a prevalence of Al–alkali coupling across our range of run conditions. The relationships between log D and H₂O content in melt for all cations in an intermediate composition are strongly nonlinear and can be fitted using an exponential function with a convergence in diffusion coefficients for different temperatures with increasing H₂O content. Thus, Arrhenius analyses result in a decrease in activation energies from 222–293 kJ mol⁻¹ at 1.7 wt % H₂O to 48–112 kJ mol⁻¹ at 3.0 wt % H₂O. These results provide new data on chemical interdiffusion in highly alkaline Na-rich melts and suggest that H₂O content plays a key role in increasing the chemical efficiency of magma mixing at low temperatures. The obtained dataset is used to test chemical controls of magma mixing in the El Abrigo ignimbrite, Tenerife, where banded pumices involving basanitic–tephritic to phonolitic magmas are common in several compositionally bimodal ignimbrite units.

1 Introduction

The study of the major and trace element composition of the melt phase is an important tool for understanding the evolution of sub-volcanic reservoirs and their pre-eruptive dynamics and timescales. Magma mixing and mingling are common and exert an important influence on the petrological and geochemical evolution of melt and whole-rock compositions, depending on the structure of the magma plumbing system and mixing endmember proportions (Sliwinski et al., 2015; González-García et al., 2022). In crystal-poor systems, the magma hybridization process is mainly controlled by the interaction of melt phases, in which chemical exchange by major and trace cation diffusion plays a key role. Although diffusion is a slow process occurring over small spatial scales, it takes advantage of extensive and complex contact areas between magmas resulting from advection, making it a key process controlling the chemical diversity of hybridizing magmas (Perugini et al., 2010, 2015; De Campos et al., 2011; DeVitre et al., 2019). Chemical variations along diffusion profiles can also provide valuable information on magma interaction and ascent timescales, as recently demonstrated by the application of major element diffusivities to banded pumices and glassy tephra (González-García et al., 2023; Shamloo and Grunder, 2023), giving a complementary insight into results from diffusion modelling on zoned minerals. Moreover, cation diffusion plays an important role in several other magmatic processes such as crystal growth and dissolution (Zhang et al., 1989; Mangler et al., 2023).

Because of this substantial field of application, the study of diffusion in melts at conditions as close as possible to those of natural volcanic systems is of major importance. Indeed, a large dataset of diffusion coefficients and Arrhenius equations has been published in the last decades (Zhang et al., 2010; Zhang and Gan, 2022), covering major elements, trace elements, and volatiles in various conditions. However, some important gaps still exist in the database. Diffusion data in Na-rich alkaline melts typical of ocean island systems are scarce, and consistent information regarding the effect of H₂O on cation diffusivity is difficult to find (Baker et al., 2002). The presence of dissolved H₂O is known to drastically reduce melt viscosity and consequently increase diffusivities (Dingwell et al., 1996; Romano et al., 2003). Similarly, high Na content and Na $/$ K ratios, typical of alkaline ocean island magmas, lead to an additional reduction in melt viscosity compared to K-rich melts in response to the modification of their molecular structure by decreasing the degree of polymerization (Le Losq et al., 2021b). In turn, these chemical variations can influence eruptive styles (Giordano and Dingwell, 2003; Andújar and Scaillet, 2012). To address these issues, we experimentally characterized the diffusive exchange of eight major elements between hydrous tephrite and phonolite melts from the Canary Islands to constrain their diffusive behaviour as a function of H₂O content and temperature on the diffusion process. These experiments are relevant to conditions prevailing during magma mixing events recognized in the Las Cañadas edifice, Tenerife (i.e. 1150–1300 °C, 0.3–3.3 wt % H₂O, and 300 MPa). The experimental data obtained provide insights into magma mixing events on the island of Tenerife and can potentially be applied to other volcanic systems of analogous composition, such as the Laacher See (Tomlinson et al., 2020) or Mayotte (Berthod et al., 2021).

2 Materials and methods

2.1 Starting materials

Two natural eruption products, which are representative of the Canary Islands volcanism, were selected as experimental starting materials (Table 1). The mafic endmember (PF21) is a tephrite sampled in the Duraznero crater, which was built during the 1949 San Juan eruption on the island of La Palma (Klügel et al., 2000; Fuchs, 2014). Such a tephritic composition is common in the historical volcanism in the Cumbre Vieja volcanic ridge (Klügel et al., 2000, 2005; Pankhurst et al., 2022), and it also resembles tephritic lava flows in the Diego Hernández Formation (DHF) of the Las Cañadas edifice, Tenerife (Bryan et al., 2002), in both major and trace element content. These tephrites represent a widespread step in the differentiation of primitive basanitic magmas in the Canary archipelago and are considered a common endmember for magma mixing events in the DHF (Wolff, 1985; Bryan et al., 2002; Sliwinski et al., 2015; González-García et al., 2022).

Table 1 Composition of starting endmembers (both the original starting glass and those from the Fe-poor experiments DC-03 and DC-04) and the TP55 intermediate composition (tephriphonolite) to which most diffusivities are referred throughout this work (see text for details). n is the number of averaged EPMA analyses.

^* H₂O concentrations in DC-03 and DC-04 were determined by Fourier-transform infrared spectroscopy (FTIR).

The evolved endmember (ABP-F) is a phonolite originating from aphyric white pumices of the El Abrigo ignimbrite (González-García et al., 2022). Phonolites of similar composition are widespread in all ignimbritic units in the DHF (Bryan et al., 2002; Martí et al., 2020), as well as in the Teide–Pico Viejo system (Ablay et al., 1998; Dorado et al., 2021), and are the result of protracted magma differentiation and storage in shallow (<5 km) reservoirs of the Las Cañadas volcano and in the currently active Teide–Pico Viejo system (Andújar et al., 2008, 2010). Phonolites are also present in La Palma, are occasionally involved in magma mixing events, and are represented in historical volcanism on the island (Johansen et al., 2005; Klügel et al., 2022).

These starting materials were crushed and powdered in a ring mill and subsequently melted at 1600 °C for 4 h in a Nabertherm^® box furnace at ambient pressure. After quenching in water, the resulting glasses were ground in an agate ball mill and melted again under the same conditions to ensure chemical homogeneity. After three melting and crushing cycles, the resulting glass was crushed again to obtain the final starting material powders for the experimental capsules.

2.2 Experimental setup

Hydrous starting glasses were produced at high pressure and temperature using Au₈₀Pd₂₀ alloy capsules with an inner diameter of 5.0 mm and a length of 35 mm. The capsules were filled with glass powder and, when necessary, appropriate amounts of distilled water, in several steps. Glasses of both endmember compositions were produced as nominally dry (ND; i.e. no added water) and at nominal H₂O concentrations of 1.5 wt % and 3 wt %. After loading, capsules were welded shut, and glasses were synthesized at 300 MPa and 1200 °C in two different internally heated pressure vessels (IHPVs; Berndt et al., 2002). Temperatures were controlled with S-type (Pt-Pt₉₀Rh₁₀) thermocouples, and pressure was monitored with calibrated Burster pressure transducers. For all runs, the capsules were quenched isobarically after switching off the power supply of the furnace, resulting in an initial cooling rate of ca. 200 °C min⁻¹. After the synthesis, the glasses recovered from the capsules were cut to obtain cylinders with a length of ca. 3–4 mm, and each cylinder was finely polished on one side.

A first set of starting glasses (experiment set A) was synthesized at fO₂ conditions close to the Ni–NiO buffer (NNO, equivalent to 0.6 log units above the quartz–fayalite–magnietite, QFM, buffer) under H₂O-saturated conditions in an IHPV equipped with an H₂ Shaw membrane for fO₂ monitoring (Berndt et al., 2002) using an Ar–H₂ pressure medium, with a synthesis duration of 72 h. In the set A synthesis experiments, the H₂ partial pressure in the IHPV was adjusted to buffer fO₂ close to the NNO equilibrium. Based on the amount of H₂ introduced into the vessel, the prevailing oxygen fugacity is estimated to be between NNO−0.1 and NNO+0.8. However, these values only correspond to experimental fO₂ if the charges are H₂O saturated. In H₂O-undersaturated runs, the final fO₂ is shifted towards more reducing values by the term 2⋅log(aH₂O), where aH₂O is the H₂O activity in the charge (Botcharnikov et al., 2008). In order to apply this correction, aH₂O was approximated by XH₂O = H₂O_measured (mol %) $/$ H₂O_saturation (mol %). H₂O saturation values were obtained using the MagmaSat solubility model (Ghiorso and Gualda, 2015). Consequently, experimental oxygen fugacity in the synthesis runs varied between NNO−1.7 and NNO−0.3 (Table S1 in the Supplement).

A second set of glasses (experiment set B) was synthesized at the intrinsic fO₂ of the IHPV (Ar-only pressure medium), with fO₂ close to NNO+2.3 at H₂O-saturated conditions, as verified using CoPd redox sensors (Taylor et al., 1992; Marxer and Ulmer, 2019). As described above, the final experimental fO₂ ranged between NNO−0.1 and NNO+1.7 as melts were water-undersaturated. Thus, in both experimental series (set A and B), the most oxidizing conditions were reached in H₂O-rich experiments.

Diffusion experiments were run using the diffusion couple method (Baker, 1989, 1991; Nowak and Behrens, 1997), in which the polished sides of phonolitic and tephritic glasses with the same nominal water content are juxtaposed with each other and placed inside an Au₈₀–Pd₂₀ capsule with an inner diameter of 5.1 mm. The denser tephritic glass was placed in the bottom position to avoid gravitational instability during the experiment. The capsules were closed by arc welding and pre-compressed in a cold-seal pressure vessel (CSPV) to check for tightness.

The diffusion experiments were run in the same IHPVs that were used for glass synthesis, at temperatures ranging from 1150 to 1300 °C and a constant pressure of 300 MPa. Table 2 shows a summary of run conditions for all experimental charges. The vessels were pressurized with Ar–H₂ (set A) or Ar (set B), and the experiments were run for 2 or 4 h (plus one additional zero-time experiment). The heating ramps from ambient temperature to final dwell temperature were 30 °C min⁻¹ in the set A experiments and 50 °C min⁻¹ in the set B experiments. After each run, a rapid-quench device (Berndt et al., 2002) was activated, allowing the capsule to fall into the cold part of the vessel and cool quickly, with quenching rates up to 50 °C s⁻¹. A total of 13 experiments were run successfully and utilized for the present study (Table 2).

2.3 Analytical procedures

From the starting glasses, 10–20 mg glass chips from both ends of the capsule were used to measure H₂O content by pyrolysis and subsequent Karl Fischer titration (KFT; Behrens, 1995; Behrens et al., 1996). The iron oxidation state (Fe³⁺ $/$ ΣFe ratio) was measured by colorimetric wet-chemical analysis (Schuessler et al., 2008). Small glass chips (5–20 mg) from the two ends of the post-synthesis glass cylinders were powdered, dissolved in HF, and analysed by absorption spectrometry in the visual spectrum using a Shimadzu UV 1800 UV–VIS spectrometer. This method allows determination of the Fe³⁺/ΣFe ratios with high precision (typically within <3 %).

After the diffusion couple experiments, each capsule was cut longitudinally in two halves. A double-polished section was prepared from one half, with a thickness between 100 and 300 µm. The second half was mounted in a 1 in. epoxy resin mount and polished. H₂O concentrations were measured by Fourier-transform infrared spectroscopy (FTIR) in the double-polished sections, and an electron probe microanalyser (EPMA) was employed to determine major element concentrations in the epoxy-mounted experimental charges.

In the resulting diffusion couples, H₂O concentrations were obtained by FTIR using a Bruker IFS88 spectrometer coupled to an IR-Scope II microscope at the Institute of Earth System Sciences (IESW), Leibniz University Hannover. Absorption coefficients determined for a Na-rich phonolite glass from Teide volcano (Carroll and Blank, 1997) were applied to quantify H₂O concentrations, namely ${\mathit{\epsilon }}_{\mathrm{4500}}=\mathrm{1.25}\pm \mathrm{0.3}$ and ${\mathit{\epsilon }}_{\mathrm{5200}}=\mathrm{1.10}\pm \mathrm{0.10}$ L mol⁻¹ cm⁻¹. No absorption coefficients are published for tephritic glass, but we used the values published for compositionally close shoshonitic glass by Vetere et al. (2011), i.e. ${\mathit{\epsilon }}_{\mathrm{4500}}=\mathrm{0.80}\pm \mathrm{0.06}$ and ${\mathit{\epsilon }}_{\mathrm{5200}}=\mathrm{1.03}\pm \mathrm{0.03}$ L mol⁻¹ cm⁻¹. To estimate melt density, we used the equations provided by Carroll and Blank (1997) and Vetere et al. (2011). Using these data, water concentrations were estimated using the modified Beer–Lambert law (Stolper, 1982).

Major element oxide concentration profiles (SiO₂, TiO₂, Al₂O₃, FeO_t, MgO, MnO, CaO, Na₂O, K₂O, and P₂O₅) were obtained with a JEOL JXA-iHP200F Hyper Probe electron microprobe equipped with five wavelength-dispersive X-ray spectroscopy (WDS) spectrometers at the IESW, Leibniz University Hannover. Operating conditions were an accelerating voltage of 15 kV, a beam current of 4 nA, and a defocused beam diameter of 12 µm to minimize alkali migration. Accuracy is less than 3 % for concentrations >1 wt % except for CaO and K₂O, which have an accuracy of <5 %. For most experiments, two parallel transects were obtained with different lengths and resolutions to capture the widely varying diffusion rates of all major elements (especially Na). In transect 1, intervals between analytical spots were kept as small as possible (typically 15–20 µm). In contrast, the second profile spanned 3 to 4 mm across the interface with separation between analytical spots in the 30–40 µm range. This second profile also allowed us to check for possible disturbances caused by convection in the experiments. Precision and accuracy were determined by measuring VG-568 (rhyolite) and VG-2 (basalt) reference glasses (Jarosewich et al., 1980; Helz et al., 2014).

2.4 Correction for experimental times

Nominal experimental dwells of 2 and 4 h were defined, and actual dwell times were registered for each experiment. However, due to the activation of the diffusion process at a certain temperature before reaching the target experimental temperature, experimental dwell times needed to be corrected by incorporating the effect of heating ramps. Thus, effective experiment durations were calculated by assuming Arrhenian behaviour of diffusivity during heat-up ramps and integrating the diffusivity as a function of T with respect to t during that interval (Zhang and Behrens, 2000). Thus, the effective heating time (t_eff) can be calculated using the following expression (Koepke and Behrens, 2001):

$$\begin{array}{}\text{(1)}& t_{\mathrm{eff}}=\underset{\mathrm{0}}{\overset{h}{\int }}{\displaystyle \frac{D_{T\left(t\right)}}{D_{\mathrm{Texp}}}}\mathrm{d}t=\underset{\mathrm{0}}{\overset{h}{\int }}\exp \left({\displaystyle \frac{E_{\mathrm{a}}}{R}}\cdot \left({\displaystyle \frac{\mathrm{1}}{T\left(t\right)}}-{\displaystyle \frac{\mathrm{1}}{T_{\exp }}}\right)\right)\mathrm{d}t,\end{array}$$

where D_T(t) is diffusivity as a function of temperature, which is itself a function of time, T(t); D_Texp is the diffusivity at the experimental dwell temperature (T_exp); E_a is the activation energy and R is the ideal gas constant; and h is the duration of the heating ramp. As a first approach, an E_a value of 150 kJ mol⁻¹ was used, which is similar to the final obtained values from the Arrhenius analysis (see Sect. 4.4), resulting in calculated effective heating times on the order of 180 to 500 s. Final effective experiment durations are listed in Table 2. The high cooling rates resulting from rapid quenching allow us to assume instantaneous cooling with no effect on effective run duration.

2.5 Determination of composition-dependent diffusion coefficients

Diffusion in natural multicomponent systems is a complex phenomenon. In such cases, the diffusion of a component not only depends on its own gradient but is also dependent on the gradients of the other components in the system (Chakraborty et al., 1995a, b; Liang, 2010). This dependence can be rigorously described using a multicomponent diffusion matrix, requiring specific experimental strategies. Although multicomponent diffusion matrixes have been obtained for up to seven- and eight-component systems (Guo and Zhang, 2016, 2018), such an approach has never been carried out for systems with large compositional gradients, and thus it is outside of the scope of this work. Instead, we here use the effective binary diffusion (EBD) approach, a simplified procedure in which the diffusivities of each component are treated as binary diffusivities, with all the other components considered to be a single component. Diffusion coefficients (D) obtained via this approach are therefore apparent diffusivities that are not necessarily comparable to those obtained in compositionally different systems.

Due to the variable bulk glass composition along the analysed profiles, a concentration-dependent procedure for determining elemental diffusivities must be used. EBD coefficients were obtained from measured concentration–distance profiles by the modified Boltzmann–Matano methodology of Sauer and Freise (1962; hereafter SF), which does not require knowledge of the Matano interface. In particular, the analytical solution for one-dimensional molar volume independent diffusion from Sauer and Freise (1962) was used:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}D\left(x\right)& ={\displaystyle \frac{\mathrm{1}}{-\mathrm{2}t\left(\partial c/\partial x\right)}}\left[\left(\mathrm{1}-c\left(x\right)\right)\underset{x}{\overset{\mathrm{\infty }}{\int }}c\mathrm{d}x+c\left(x\right)\right.\\ & \left.\underset{-\mathrm{\infty }}{\overset{x}{\int }}\left(\mathrm{1}-c\left(x\right)\right)\mathrm{d}x\right],\end{array}\end{array}$$

where c(x) is the normalized concentration of the diffusing component, with c(x)=1 at $x=-\mathrm{\infty }$ and c(x)=0 at $x=+\mathrm{\infty }$. A modification of the Python programming language script (Oliphant, 2007) used by González-Garcia et al. (2018) was used to calculate composition-dependent diffusivities from the SF procedure.

The SF method requires that the variation in the molar volume of the melt along the diffusion profile is non-significant. To check this point, we used the Bouhifd et al. (2015) formulation to obtain the molar volume of tephritic and phonolitic melts under our experimental conditions, resulting in a >8 % variation. This value is within the error in the SF method itself and will have an insignificant impact on obtained diffusivities. The error in D was estimated by simple error propagation from uncertainties in involved data, including molar volume variation. Another limitation of the SF method is that D cannot be accurately obtained at the endmember compositions, and for this reason we restrict the calculation of diffusivities in our profiles from 20 % to 80 % of the endmember compositional range.

3 Results

3.1 Starting glasses and diffusion couples

All synthesis experiments resulted in homogeneous, bubble-free glasses. Phonolitic glasses were crystal-free; however, tephritic glasses showed variable numbers of dendritic quench crystals (<25 µm), likely clinopyroxene, produced during the experiment cool-down phase. In contrast, diffusion couple runs resulted in crystal- and bubble-free glasses, as observed in scanning electron microscope (SEM) backscattered electron (BSE) imaging. This was corroborated by a zero-time experiment performed at 1200 °C and 0.3 wt % H₂O, which showed that all quench crystals formed during the synthesis phase of the tephritic glasses effectively dissolved upon heating and before reaching the experimental temperature (Fig. S1 in the Supplement).

H₂O content measured by KFT and FTIR in starting materials and diffusion couples agrees well with intended concentrations, except for the nominally dry glasses of experiment set A. In these experiments, the water content measured in starting glasses was 0.8 wt %–1.1 wt %, resulting from a reduction in Fe₂O₃ in the starting material and H₂ influx through the capsule walls. In experiment set B, the same effect resulted in H₂O content of 0.3 wt %–0.5 wt %. In the diffusion couples, H₂O concentrations, measured by FTIR along a line centred on the interface of the experiments, proved to be constant or close to constant for each experiment. Variations along the line were not higher than 10 %, and hence the average and standard deviation of analysed spots were assumed to be the H₂O concentration value for each experiment (Table 2).

Measured molar Fe³⁺ $/$ ΣFe ratios in starting glasses are provided in the Supplement (Table S1). Both set A tephritic and phonolitic glasses show a similar trend, with Fe³⁺ $/$ ΣFe varying from 0.1 to 0.2 with increasing H₂O content (see Fig. S3). However, for set B glasses, the observed trends differ for phonolite and tephrite, with Fe³⁺ $/$ ΣFe varying from 0.25 to 0.5 for the tephrite and from 0.45 to 0.65 for the phonolite and again with an increasing trend with increasing H₂O owing to the more-oxidizing conditions in H₂O-rich melts. This is a result of the fO₂ buffering mechanism in the IHPV. Synthesis runs with low H₂O content performed under reducing conditions suffered from loss of iron to the Au–Pd alloy capsule (Fig. S4), resulting in a relative loss of ca. 50 % and 80 % of FeO in phonolitic and tephritic glasses, respectively (Table 1). This effect is only significant in the two ND experiments and diminishes progressively with increasing fO₂. Experiments run at intrinsic vessel fO₂ do not show significant Fe loss (Fig. S4).

Figure 1 Total alkali–silica diagram showing the compositions of initial tephritic (PF21) and phonolitic (ABP-F) endmembers, as well as Fe-poor endmembers from experiments DC-02 and DC-03 (PF21-R and ABP-F-R). The intermediate composition (TP55) for which diffusivities are derived is also shown. Glass compositions from El Abrigo pumices (González-Garcia et al., 2022) are plotted for comparison.

Figure 2 (a) Major element oxide concentration–distance profiles measured in experiment DCB-09, with concentrations normalized for endmember concentrations (C_min=0 and C_max=1). (b) MnO profile showing prominent uphill diffusion. (c) Full Na₂O profile showing the greater extension of the diffusive zone (note the different length scale). Where not shown, error bars are smaller than the symbol.

3.2 Diffusion data

Concentration–distance profiles of 9 of the 10 major oxides measured in the diffusion couple runs resulted in monotonic diffusion profiles with no evidence of uphill diffusion. Figure 2 shows an example of concentration–distance profiles measured with the EPMA, and the complete dataset is included in the Supplement File 2. As expected from compositional contrast, all profiles are asymmetric and extend deeper into the tephritic half of the diffusion couple except for Na, which shows the opposite behaviour. This dataset allows an effective binary diffusion (EBD) approach (Zhang et al., 2010). In contrast to other oxides, the MnO profile shows strong signs of uphill diffusion (Fig. 2b) resulting from multicomponent diffusion effects and thus was left out of the EBD analysis.

Figure 3 Concentration dependence of major element diffusivities in tephrite–basanite couples. (a) Diffusion coefficient vs. SiO₂ concentration along run DCB-09, resulting from the SF analysis. The thick vertical line indicates the TP55 composition used as a proxy for diffusion coefficients throughout the text. (b) Variation in the compositional dependence of diffusion coefficients with H₂O content at 1250 °C, using the ratio $D_{\mathrm{51}}/D_{\mathrm{59}}$ (ratio between diffusivities at 51 wt % and 59 wt % SiO₂) as an index.

The application of the SF equation to the diffusion profiles led to continuous variations in D with distance across the profile (Fig. 3) and hence with bulk glass composition. As already noted by observing the shape of the concentration–distance profiles, most components show an inverse correlation between SiO₂ content and D, with the only exception being Na, which shows the opposite behaviour. Due to the intrinsically fast diffusivity of Na, its diffusion coefficient could only be determined in the low-temperature and/or low-H₂O experiments (seven runs). In the remaining experiments, Na diffusion reached the ends of the experimental charges, and therefore the experiment lost its characteristic of a semi-infinite medium, leading to strongly underestimated Na diffusivity. In the experiments where D was determined for Na, it showed the fastest diffusion among major elements, being consistently ∼1 order of magnitude faster than Si. Surprisingly, Al has the second fastest diffusivity, falling approx. 0.5–0.6 log units below Na and ∼0.4 log units above Si for all experiments. This is an unexpected result, since the network-forming character of Al³⁺ cations should result in diffusivities close to those of Si and Ti, as observed in most data in the literature (see Sect. 4.5). The Mg–Fe–Ca–K group displays similar diffusivities but always falls between Al and Si. Finally, Ti is the slowest diffusing component, showing values similar to or lower than Si. In comparison to the diffusion couple experiments carried out by González-Garcia et al. (2017) using shoshonite–rhyolite couples, results from tephrite–phonolite couples are highlighted by (1) a lesser prevalence of diffusive coupling to Si, resulting in a wider dispersion of diffusivities, and (2) a lack of uphill diffusion phenomena, with the only exception being Mn.

4 Discussion

4.1 Compositional dependence

One of the main factors influencing cation diffusion in silicate melts is the melt chemical composition. Diffusion of divalent, trivalent, and tetravalent cations is usually faster in silica-poor depolymerized melts than in silica-rich polymerized melts, while the reverse is usually true for monovalent cations (e.g. Zhang et al., 2010). As already visible in the asymmetry of composition–distance profiles, there is a clear compositional dependence in our tephrite–phonolite couples, and the SF analysis clearly shows this effect (Fig. 3a). To quantify the degree of compositional dependence of major cation diffusion, we can use the ratio of the diffusion coefficients for 51 wt % SiO₂ and 59 wt % SiO₂ ($D_{\mathrm{51}}/D_{\mathrm{59}}$). For the 1250 °C experiment series, we find that compositional dependence is higher in the nominally dry experiments, where $D_{\mathrm{51}}/D_{\mathrm{59}}$ values are between 1.5 and 5 (equal to a 0.2 to 0.7 log unit difference between the tephritic and phonolitic compositions; Fig. 3b). For H₂O content more than 1 wt %, the ratio $D_{\mathrm{59}}/D_{\mathrm{51}}$ decreases to 1.1–1.5 (<0.2 log unit difference). Again, the exception is Na, which follows the expected behaviour for monovalent cations and diffuses up to 0.2 log units faster in phonolite than in tephrite.

Due to the rather mild compositional dependence for most elements, we have chosen an intermediate composition equivalent to 55 wt % SiO₂ (tephriphonolite) as an index, for which the diffusive behaviour of major components will be described from now on. This melt composition will be hereafter referred to as TP55 (Table 1), and the full dataset of diffusion coefficients is given in Table 2.

Table 2 Experiment parameters and diffusivities obtained for the TP55 composition.

All diffusivities correspond to a composition of 55 wt % SiO₂ (19.7 mol % Si), referred to as TP55 in the text. D is in m² s⁻¹. ^a H₂O content is the average of several FTIR measurements in each experiment. ^b XH₂O is calculated as H₂O(experiment) $/$ H₂O(saturation) in mol %. H₂O saturation was calculated using the MagmaSat model (Ghiorso and Gualda, 2015).

Figure 4 Mg-normalized Fe diffusivities ($D_{\mathrm{Fe}}/D_{\mathrm{Mg}}$) as a function of (a) experimental H₂O content and (b) the estimated Fe³⁺ $/$ ΣFe ratio of each experiment. The discontinuous line is a linear fit of the data.

4.2 Fe diffusion

The varying fO₂ experimental conditions resulted in variable Fe oxidation states (Fe³⁺ $/$ ΣFe ratio) in experimental sets A and B, as already observed in Fig. S2. The prevailing structural role of Fe²⁺ as a network modifier and Fe³⁺ as a weak network former (with some Fe³⁺ possibly acting as a network modifier; Le Losq et al., 2021a; Moretti and Ottonello, 2022) may result in different diffusivities for both species, although the uncertainty in the behaviour of Fe³⁺ does not allow us to make a reliable prediction. In any case, due to its higher bond strength, Fe³⁺ diffusivities are expected to be lower than those of Fe²⁺ (Zhang et al., 2010). Hence, the measured D of total Fe is expected to be between that of Fe²⁺ and Fe³⁺. Our measured diffusion dataset can be used to put some first-order constraints on the differential diffusivity of Fe²⁺ and Fe³⁺in the tephrite–phonolite system. Since Fe³⁺ $/$ ΣFe was not determined in the diffusion couples, we used its relationship with estimated fO₂ in the synthesis runs (Fig. S3) to obtain a tentative estimate of Fe³⁺ $/$ ΣFe in the diffusion couples (Fig. 4). To avoid the effects of different experimental conditions, Fe diffusivities from each run were normalized to those of Mg from the same experiment. The results suggest that the effect of oxygen fugacity in Fe diffusion is resolved in our experimental dataset (Fig. 4a). $D_{\mathrm{Fe}}/D_{\mathrm{Mg}}$ values for set B are below 0.9, while in the more-reducing experiments (set A), $D_{\mathrm{Fe}}/D_{\mathrm{Mg}}$ increases to 0.9–1.3. From these data and the linear fit in Fig. 4b, the difference in diffusivity of Fe³⁺ and Fe²⁺ can be estimated. Although the associated errors are large, our results suggest that Fe²⁺ diffusivity in the tephrite–phonolite system is on the order of 3 to 5 times faster than that of Fe³⁺ for a given temperature and given H₂O content (Fig. 4b). Therefore, our dataset seems to indicate a measurable difference in diffusivity between both Fe species in alkaline silicate melts.

4.3 H₂O dependence

Dissolved H₂O in melt has already been demonstrated to strongly enhance diffusivities in silicate melts, for both major elements (e.g. Baker, 1991; Baker and Bossányi, 1994; González-García et al., 2017, 2019) and trace elements (Watson, 1981; Baker and Bossányi, 1994; Behrens and Zhang, 2001; Zhang et al., 2010; González-Garcia et al., 2018; Spallanzani et al., 2022). Our experimental dataset confirms this observation, and for the TP55 composition major element diffusion is enhanced by 1.2–1.5 log units (corresponding to a factor of 20–30) when H₂O content increases from 0.3 wt % to 3.3 wt % at 1250 °C. At 1200 °C, the enhancement factor amounts to 0.8–1.0 log units for an increase in H₂O content from 1.1 to 3.0 wt % H₂O (Fig. 5). For the dataset obtained at 1150 °C, with two successful runs at 1.6 wt % and 3 wt % H₂O, the enhancement is on the order of 0.5 log units for all elements. Here, runs at different fO₂ levels were treated together since oxygen fugacity is not expected to have a major effect on diffusivity, except for Fe (see above).

Figure 5 Water-dependence of major element effective binary diffusivities for a 55 wt % SiO₂ compositional term.

Table 3 Parameters a, b, and c obtained from fitting the H₂O-dependent diffusivities with an exponential law in the form log $D=a-b\cdot e^{c\cdot w}$, where w is the H₂O content in weight percent. RMSE: root-mean-square error.

Moreover, the increase in D as a function of H₂O in weight percent is clearly nonlinear, and instead it can be modelled by an exponential function in the form $\log D=a-b\cdot e^{c\cdot w}$, where a, b, and c are fitting parameters, and w is the H₂O concentration given in weight percent. Parameter a gives the asymptotic log D value, i.e. the diffusivity at the saturation level. This fitting equation was applied to our diffusion datasets at 1200 and 1250 °C using an instrumentally weighted Levenberg–Marquardt least-squares algorithm in the software QtiPlot v. 1.0.0. Fitting parameters for each component are provided in Table 3. The 1200 °C dataset, spanning 1.1 wt % to 3.0 wt % H₂O, can be fitted with a root-mean-square error (RMSE) between 0.02 and 0.13 log D units (with D in m² s⁻¹). The 1250 °C dataset, spanning 0.3 wt % to 3.2 wt % H₂O, yielded an RMSE between 0.04 and 0.09 log D units. Considering H₂O content in mole percent instead of weight percent did not change the fitting shape.

Such strongly nonlinear H₂O–diffusion relationships are common in the literature (e.g. Baker, 1991; Watson, 1994; Baker et al., 2002; Zhang et al., 2010), but fitting equations are rarely provided. For small variations in H₂O content (0.3 wt %–2 wt %), González-García et al. (2017, 2018) found linear relationships for both major and trace elements in shoshonite–rhyolite couples, but the inclusion of diffusivities under dry conditions seemed to make a square-root function more appropriate (González-García et al., 2019). However, different behaviours arise when the effects of water species OH and molecular H₂O (H₂O_m) are considered. FTIR analysis shows that the proportion of OH and H₂O_m varies with total H₂O in melt, with OH being the predominant species at H₂O $<\sim \mathrm{2}$ wt % (∼8 mol %). This trend (Fig. S5) is equivalent to those established for other melts (e.g. Lesne et al., 2011). We find that log D follows a similar exponential function with H₂O_m (in mol %) to that found for total H₂O, and a clear relationship with OH cannot be established. However, the best fit comes with the molar ratio [OH] $/$ [H₂O_m], where linear correlations can be observed for all major oxides for H₂O >1 wt % (Fig. S6).

It is remarkable to note the convergence of elemental diffusivities towards high H₂O content in both datasets, which also includes the data point at 1150 °C and 3 wt % H₂O. All experimentally measured diffusivities at such high H₂O concentrations are within 0.4 log units, compared to a variation of 0.6–0.8 log units at 1.6 wt % H₂O. A consequence of this behaviour is that the diffusivity of the slowest components (in our case, Ti and Si) is more sensitive to the addition of H₂O than the fastest ones are, confirming previous observations (e.g. Watson, 1994). The exponential dependence of log D on H₂O content is likely controlled by melt viscosity, which follows a similar non-linear evolution with H₂O content. Our observation suggests that H₂O plays a major role in increasing the efficiency of chemical mixing in interacting melts, virtually equalizing diffusivity variations across a temperature range. Therefore, water-rich magmatic environments, such as those present in the Canary Islands, could be more susceptible to developing efficient mixing between melts than water-poor melts, not only diminishing the viscosity but also reducing the viscosity difference between endmembers. On the contrary, H₂O-poor or dry magmas with higher viscosity and slower diffusion kinetics would favour physical mingling over chemical hybridization.

4.4 Arrhenius relations

The temperature dependence of diffusion is well established and is determined by the Arrhenius equation:

$$\begin{array}{}\text{(3)}& D=D_{\mathrm{0}}\times e^{-E_{\mathrm{a}}=}/RT,\end{array}$$

where D₀ is the pre-exponential factor (m² s⁻¹), E_a is the activation energy of the diffusion process (J mol⁻¹?), R is the gas constant (J mol⁻¹ K⁻¹), and T is temperature (K). By plotting 1000 $/$ T vs. ln D and fitting measured diffusivities with a linear regression, it is therefore possible to calculate D₀ and E_a and characterize the T dependence of major element diffusion. The dataset was fitted using the same procedure and software detailed above for obtaining H₂O vs. D equations.

From our experimental dataset, Arrhenius equations relating temperature and diffusivity were obtained under specific conditions. The major parameter influencing diffusion in our dataset is H₂O content, but its measured value varies up to 0.1 wt %–0.3 wt % for experiments with the same nominal H₂O content. Therefore, a correction of diffusion coefficients to a common H₂O concentration must be made before plotting the Arrhenius relations. We used the H₂O-dependent equations obtained in the previous section (Table 3) to obtain diffusivities at 1200 and 1250 °C for fixed H₂O content of 1.66 wt % and 2.99 wt %. These values were determined by the H₂O content of the experiments at 1150 °C, where an H₂O-dependent equation was not obtained.

Figure 6 Arrhenius relations of TP55 for six major elements at 1.66 wt % and 2.99 wt % H₂O. Due to the anomalous behaviour, a linear fit was not obtained for Al at 2.99 wt % H₂O. The determined ln D values of the 0.3 wt % series are plotted for reference.

Figure 6 shows Arrhenius fits for Si, Ti, Al, Mg, Ca, and K, and the parameters obtained are summarized in Table 4. Fe is not included due to inconsistencies arising from the different Fe³⁺ $/$ ΣFe ratios in experiments, and Na data are insufficient to obtain Arrhenius parameters of good quality. Results show a coherent picture, indicating a decrease in E_a with increasing H₂O content from 220–290 kJ mol⁻¹ at 1.66 wt % H₂O to 47–112 kJ mol⁻¹ at 2.99 wt % H₂O. Activation energies and ln D₀ for both H₂O concentrations plot to a single compensation law linear fit (Hart, 1981; Fig. S7). This is in line with previously published data, where H₂O content is known to reduce the activation energy of diffusion (Watson, 1981, 1994; Behrens and Zhang, 2001; Spallanzani et al., 2022). Only Al shows anomalous behaviour, where the fitting errors are larger for the 1.66 wt % H₂O fit, and the 2.99 wt % H₂O dataset does not allow us to obtain a good Arrhenian fit.

Table 4 Arrhenius parameters calculated for 1.66 wt % and 2.99 wt % H₂O.

Figure 7 (a) Al and (b) Si effective binary (EB) diffusion data obtained in this work compared to bibliographic data. The dashed coloured lines indicate Eyring diffusivity for TP55 at H₂O content of 0.3 wt %, 1.7 wt %, and 3.0 wt %, with viscosities calculated using the Giordano et al. (2008) model. Bibliographic data are from (1) Baker and Bossányi (1994), (2) Chen and Zhang (2008, 2009), (3) Lundstrom (2003), (4) Zhang et al. (1989), (5) González-García et al. (2019), and (6) Koyaguchi (1989).

4.5 The case of aluminium

Throughout our complete dataset, Al is the second-fastest-diffusing element after Na, and Al diffusivities are 0.4 log units faster than Si. This behaviour can be considered anomalous, given that the expected network-forming character of Al³⁺ in the silicate network should result in diffusivities comparable to those of Si⁴⁺ and Ti³⁺. A comparison of our Al diffusivities with published data (Fig. 7) confirms this anomaly. Al diffusion in TP55 melt composition with 2.7 wt % to 3.3 wt % H₂O plots between 1 and 2 orders of magnitude above diffusivities measured by Baker and Bossányi (1994) in wet dacite with the same range of H₂O content and is comparable to those observed in dry basalts and basanites at 1500 °C (Lundstrom, 2003; Chen and Zhang, 2008, 2009), i.e. 200 to 350 °C hotter than our experimental conditions. This behaviour is not observed in Si, for which diffusivities do not differ strongly from data in the literature for similar water content. The comparison to Eyring diffusivity (i.e. viscosity-related diffusion usually applied to network-forming cations; Glasstone et al., 1941) of TP55 at different levels of H₂O content yields similar differences, with Al EBD plotting significantly above Eyring D values for the corresponding melt H₂O content (Fig. 7). Although Si data also plot above Eyring D, their deviations are lower and consistent with previous studies (Fanara et al., 2017). The large departure of Al from Eyring diffusivity also contrasts with the comparatively low deviations observed by Fanara et al. (2017) for trivalent cations. The only dataset supporting high Al mobility in silicate melts comes from dynamic-mixing experiments performed with K-series alkaline melts from Campi Flegrei. In such experiments, Perugini et al. (2013, 2015) showed that the relaxation of concentration variance (R, which is highly dependent on diffusivity) of Al is larger than that of most major components, falling between that of Na and K. On the other hand, this behaviour is not observed in chaotic-mixing experiments with subalkaline melts (Morgavi et al., 2013), where Al mobility is similar to those of Si and Ti. Overall, these observations point to a fast Al diffusion mechanism operating in highly alkaline melts, in agreement with our diffusion experiments.

In the silicate melt framework, Al usually plays the role of a network-forming cation (Mysen et al., 1981), and as such its diffusivity is expected to be comparable to that of other network-forming cations in tetrahedral (i.e. 4-fold) coordination (Si, Ti) and to Eyring diffusivity. Departures from this behaviour have been suggested for high-P and peraluminous magmas, suggesting that Al may not be in tetrahedral coordination if the local charge balance is not attained with M⁺ (Na, K) and M²⁺ (Ca) cations (Mysen et al., 1981). In molecular dynamics simulations, alkali cations (especially Na⁺) tend to cluster around Al³⁺ to compensate for the charge balance (Le Losq et al., 2017), which could lead to increased Al diffusivity, where highly mobile Na cations result in an enhancement effect in Al mobility in the melt. This effect is greatly facilitated in melts with a large Na $/$ K ratio, while in K-rich melts the large size of K⁺ promotes polymerization, which in turn increases viscosities and decreases diffusivities.

A complementary explanation for the observed Al behaviour could potentially reside in multicomponent diffusion effects. Diffusion in complex multicomponent systems is subject to diffusive coupling between oxides, usually resulting in uphill diffusion of one or more components (e.g. Guo and Zhang, 2016, 2018). This is the case for Al in basalt–rhyolite couples or in other systems with similar endmembers (e.g. Koyaguchi, 1989; González-García et al., 2017), but Al uphill diffusion is notably absent in our tephrite–phonolite couples. This behaviour is, however, expected, since the occurrence of uphill diffusion is highly dependent on the arrangement of endmembers in the compositional space, as observed in simple three-component systems (Chakraborty et al., 1995b). Multicomponent diffusion experiments in the system K₂O–Al₂O₃–SiO₂ (Chakraborty et al., 1995a, b) suggest again that Al diffusion is coupled to alkalis. The degree of coupling depends on composition (peralkaline vs. peraluminous) and the nature of the gradient of all remaining elements. In our system, an indication of Al–K diffusive coupling comes from their behaviour in the Arrhenius plots (Fig. 6). Al and K show here anomalous behaviour, both deviating from linearity. Diffusivities of Al and K show an increase from 1200 to 1150 °C, especially evident in the 3.0 wt % H₂O data. This behaviour could potentially be related to a change in the Al–K diffusive coupling at lower temperatures due to the variation in the intrinsic diffusivity of both elements.

Figure 8 Comparison between experimental diffusivities for TP55 composition and normalized concentration variance (${\mathit{\sigma }}_n^{\mathrm{2}}$) for three mixed clasts from the El Abrigo ignimbrite. Plotted diffusivities are from experiment DCB-11 (1200 °C, 2.93 wt % H₂O). The inset shows a microphotograph of a mingling and mixing area in El Abrigo banded pumice. The scale bar is 200 µm long.

4.6 Implications for magma mixing in Tenerife

The experimental results obtained in this work can be used to obtain some insights into the petrogenesis of bimodal ignimbrite sheets in the Diego Hernández Formation in Tenerife. In several of these units, the presence of banded pumices suggests the occurrence of magma mixing and mingling shortly before the eruption (González-García et al., 2022). During magma mixing, the combined effects of diffusion and advection progressively reduce the compositional variability in the system and tend to produce a hybrid composition. Perugini et al. (2015) proposed that the concentration variance normalized to that of the starting endmembers (${\mathit{\sigma }}_n^{\mathrm{2}}$) of elemental concentrations in the mixing melts is reduced with time, following an exponential decay function with a decay parameter R that depends, among other parameters, on the diffusion rate of the element considered. In consequence, in systems where magma mixing is the primary process controlling compositional variability, ${\mathit{\sigma }}_n^{\mathrm{2}}$ can be used as a measure of the chemical hybridization of the system.

Here we use the compositions of groundmass glass in three crystal-bearing (<30 % crystals) banded pumice clasts from the El Abrigo ignimbrite, Tenerife (González-García et al., 2022), to obtain the values of ${\mathit{\sigma }}_n^{\mathrm{2}}$ for Si, Ti, Al, Fe, Mg, Na, and K. Concentration variance was normalized using the compositions of the two endmember melts proposed by González-García et al. (2022), i.e. a phonotephrite represented by the most mafic glass in banded pumices and a high-Zr phonolite from aphyric white pumices. The obtained ${\mathit{\sigma }}_n^{\mathrm{2}}$ values were later compared to diffusivities obtained in run DCB-11 (1200 °C, 2.93 wt % H₂O). The choice of this experiment for comparison is based on the availability of Na diffusivities and on its high H₂O content, close to that expected for the phonolitic melt of the El Abrigo ignimbrite, but analogous results are obtained using any other experiment. However, we should note that experimental temperatures are still 150–200 °C higher than expected for the natural system.

The results (Fig. 8) show a visible correlation between diffusion coefficients and ${\mathit{\sigma }}_n^{\mathrm{2}}$, with slow-diffusing elements showing a lesser homogenization degree than fast-diffusing elements (particularly Na). However, Al and K outliers are present in this trend. The Pearson correlation coefficients for the three products range between 0.21 and 0.45 (in absolute value), but when the Al and K outliers are removed, the correlation improves to 0.84–0.88. These data confirm that the main control of melt chemistry in the banded pumices is the diffusive exchange between endmembers. However, the presence of Al and K outliers is an anomaly in this context and might provide additional insight into the issue of Al diffusivity. Figure 8 shows that the ${\mathit{\sigma }}_n^{\mathrm{2}}$ value of Al in the mixed-eruption products does not correlate with the fast Al diffusivity measured in our experiments but instead suggests a slow diffusion mechanism in the natural system. In one of the banded pumices, K follows similar behaviour, with one of the ${\mathit{\sigma }}_n^{\mathrm{2}}$ values higher than expected from its measured diffusivity. The experimental results from this work and from Perugini et al. (2013), where high Al mobility was observed, were carried out at temperatures between 1150 and 1300 °C. In contrast, mixing temperatures in the El Abrigo system were significantly lower, likely in the range of 900–1050 °C (González-García et al., 2022) and with the presence of significant proportions of crystals. Experimental studies in simplified Ca- and Na-aluminosilicate glasses have shown that the coordination number of Al³⁺ is temperature dependent, with lower temperatures favouring 4-fold coordinated Al (^IVAl) over 5-fold coordinated Al (^VAl; Stebbins et al., 2008; Le Losq et al., 2014). In such aluminosilicate systems, it has been suggested that the formation of ^VAl results in a reduction in melt viscosity (e.g. Kim et al., 2022) and likely Al-alkali clustering, while the opposite is true for ^IVAl-dominated melts. We speculate that a variation in the coordination number of Al towards a predominance of ^IVAl at lower temperatures, like that estimated for the natural system, could potentially result in additional viscosity increase and/or decreased Al-alkali clustering, favouring a comparatively slower Al diffusion compared to at higher temperatures.

Alternatively, the reasons for this behaviour could reside in crystal dissolution and assimilation processes in the natural products. The banded pumices in the El Abrigo deposit are mostly rich in crystals, with two coexisting mineral assemblages in disequilibrium. The partial melting and assimilation of enough quantities of alkali feldspar, as documented in several ignimbrites in the DHF, could potentially result in local Al–K inhomogeneities, which increase their concentration variance relative to a mixing-only scenario. However, our current dataset does not allow us to distinguish between these two possible scenarios, and more experimental data at lower temperatures would be necessary to confirm the occurrence of varying diffusive mechanisms for Al.

5 Conclusions

Diffusion coefficients of major elements have been obtained from tephrite–phonolite couples with compositions relevant to Canary Islands magmatism. Diffusivities show a wide variation range, with Na and Al being significantly faster (ca. 0.5 and 1 log unit, respectively) than Si. All log D values for an intermediate composition (55 wt % SiO₂) have a strongly non-linear relationship with dissolved H₂O, which is best captured by an exponential expression. In addition, an apparent convergence of diffusivities is evident towards the water-rich end (3 wt %) of our dataset, which suggests that the addition of H₂O could increase magma mixing efficiency by equalizing diffusivities of major elements. Arrhenius equations obtained at 1.66 wt % and 2.99 wt % are also convergent with increasing temperature, and E_a values decrease with H₂O content. The anomalous behaviour of Al, which is second in diffusivity only after Na, may be explained by a change in behaviour from network-forming to network-modifying cations by complexation with Na⁺ and K⁺ cations. Such behaviour may be a common occurrence in highly alkaline melts. This result suggests the occurrence of a fast diffusion mechanism for Al in our experimental dataset, different from that observed in the literature for sub-alkaline and mildly alkaline melts.

Our results may have important implications for the study of kinetic processes in highly alkaline magmas and particularly for melt hybridization during magma mixing events. The comparison of our experimental results with the chemical distribution in mixed pyroclasts from the El Abrigo ignimbrite confirms that diffusion is the main mechanism controlling the chemical variability in the melt phase. However, the natural Al distribution does not support the fast Al diffusion mechanism observed in the experiment. We suggest that at the low temperatures inferred for the mixing event (900–1050 ° C), Al does not show the fast diffusivities observed in experiments, possibly due to a change in Al coordination and hence a variation in its coupling to Na⁺ and K⁺. Alternatively, this could be an effect of the melting and assimilation of a K-feldspar rich mush, which is documented from the El Abrigo eruption.