# The role of magma ascent rates and viscosity in explosive and dome eruptions (Kelud volcano, Indonesia)

Sri Budhi Utami, Fidel Costa, Hanik Humaida

► To cite this version:

Sri Budhi Utami, Fidel Costa, Hanik Humaida. The role of magma ascent rates and viscosity in explosive and dome eruptions (Kelud volcano, Indonesia). Bulletin of Volcanology, 2024, 86, 10.1007/s00445-023-01698-1. insu-04462219

HAL Id: insu-04462219

https://insu.hal.science/insu-04462219v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The role of magma ascent rates and viscosity in the explosive and dome eruptions (Kelud volcano, Indonesia)

Sri Budhi Utami (✉ sri.budhi.utami@gmail.com)

University of Liege: Universite de Liege https://orcid.org/0000-0003-4338-0157

Fidel Costa

Institut de Physique du Globe de Paris

Hanik Humaida

Balai Penelitian dan Perkembangan Teknologi Kebencanaan Geologi

## Research Article

**Keywords:** Kelud, Kelut, magma ascent rates, timescales, viscosity, apatite

**Posted Date:** August 16th, 2023

**DOI:** https://doi.org/10.21203/rs.3.rs-3133512/v1

**License:** cc © This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Version of Record:** A version of this preprint was published at Bulletin of Volcanology on January 12th, 2024. See the published version at https://doi.org/10.1007/s00445-023-01698-1.
The role of magma ascent rates and viscosity in the explosive and dome eruptions (Kelud volcano, Indonesia)

Sri Budhi Utami<fnref n="1" />,<fnref n="2" />*, Fidel Costa<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, Hanik Humaida<fnref n="5" />

<footnote n="1">¹ Asian School of the Environment, Nanyang Technological University, 50 Nanyang Avenue, N2-01B-27, 639798, Singapore. ORCID: 0000-0003-4338-0157.</footnote>

<footnote n="2">²Present address: Département de Géologie, Bâtiment B20, Quartier Agora, Allée du six Août 12, 4000 Liège 1, Belgium</footnote>

<footnote n="3">³Earth Observatory of Singapore, 50 Nanyang Avenue, N2-01A-15, Singapore, 639798.</footnote>

<footnote n="4">⁴Université de Paris Cité, Institut de Physique du Globe de Paris, CNRS, F-75005, Paris, France. ORCID: 0000-0002-1409-5325</footnote>

<footnote n="5">⁵Balai Penyelidikan dan Perkembangan Teknologi Kebencanaan Geologi, Jalan Cendana, 15, D.I. Yogyakarta, Indonesia.</footnote>

<footnote>* Correspondence: sri.budhi.utami@gmail.com</footnote>
1 Abstract

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

Magma ascent rate can control the hazard potential of an eruption, but it is difficult to directly determine. Here we investigate the variations in timescales of magma ascent and rates of magma ascent for the three most recent explosive and effusive eruptions of Kelud volcano in Indonesia (1990, 2007, and 2014) using the zoning of volatile elements (OH, Cl, F) in apatite.

We found that crystals from the 2007 dome show chemical gradients and increasing concentrations (reverse zoning) in chlorine and/or fluorine towards the crystals' rims whereas those of the 1990 and 2014 explosive eruptions are unzoned. Diffusion modelling of the volatile elements in zoned apatite of the 2007 dome rocks give magma ascent times of up to 3 months, although 71% of them are ≤ 60 days. In contrast, the maximum magma ascent timescales inferred from apatite of the 1990 and 2014 explosive eruptions are about 7-8 hours.

Using the pre-eruptive magma storage depths obtained from petrological and phase equilibria studies, we calculate ascent rates of about > 0.4 × 10⁻³ m s⁻¹ for the 2007 dome, and > 2.6 × 10⁻¹ m s⁻¹ for the 1990 and 2014 eruptions. We also calculated the magma viscosities for each eruption, which when combined with the magma ascent rates and magma mass discharge rates correspond well with the expected eruptive styles. Our study illustrates the robustness of modelling apatite zoning in volatile elements to constrain timescales and magma ascent dynamics, and highlights the important role of magma ascent on eruptive styles.

**Keywords:** Kelud · Kelut · magma ascent rates · timescales · viscosity · apatite.
## Introduction

21

The eruptive styles of active volcanoes are difficult to anticipate, but they are key for proper hazards assessment and crises response of impending eruptions. At many active volcanoes, eruptions can transition from effusive to explosive and vice versa between eruptions for example at Merapi (Costa et al. 2013; Preece et al. 2014), or within one eruptive period such as at Rabaul (Bernard et al. 2022), Redoubt (Bull and Buurman 2013), and more recently at Soufriere St. Vincent (Joseph et al., 2022). Such transitions in eruption style depend on several parameters, including melt volatile contents, melt compositions, magma rheology, and/or magma ascent rates, and these can rapidly change during transport in the conduit and ultimately determine whether magma “flows or blows” (Dingwell 1996; Rutherford 2008; Cassidy et al. 2018). Among the parameters that affect magma explosivity, magma ascent rates appear to exert critical control (Cassidy et al. 2018), but they are generally difficult to determine. Using geophysical data, the timing, and the location of seismicity that corresponds to the run-up to eruption have been used to obtain constraints on magma ascent rates (e.g. Passarelli and Brodsky 2012; Rasmussen et al. 2018; Roman and Cashman 2018), but most estimates are from experimental and petrological studies. These include quantitative analysis of crystal size distribution and microlites density (Couch et al. 2003; Rutherford 2003; Toramaru et al. 2008), and melt decompression experiments (Martel and Schmidt 2003). In addition, chemical zoning of volatile elements in minerals or melt embayments can be modelled using diffusion laws, from which timescales and magma ascent rates can be obtained (e.g. Myers et al. 2019; Li et al. 2020; Georgeais et al. 2021). Combination of micro-textures and geochemistry of the crystal cargo can yield timescales that are comparable to some volcano monitoring techniques (e.g. seismicity) (e.g. Saunders et al. 2012; Kilgour et al. 2014; Albert et al. 2015). These timescales of magma ascent can then be converted to ascent rates assuming magma storage depths are also known. Such depths below the volcano can be obtained from phase equilibria experiments (e.g. Andújar et al. 2013; Cadoux et al. 2014; Utami et al. 2022), or from mineral-melt geobarometry (Blundy and Cashman 2008; Putirka 2008). Alternatively, geophysical studies employing tomography can be used to calculate storage depth of many volcanoes (Puspito and Shimazaki 1995; Widiyantoro and van der Hilst 1996; Widiyantoro et al. 2018), although the uncertainty of magma storage depths of some of these studies yield larger ranges of potential magma ascent rates.

52
Among the various petrological tools available to constrain volatile budgets and magma ascent rates, apatite is an increasingly popular mineral (e.g. Stock et al. 2016; Humphreys et al. 2021; Bernard et al. 2022), as it can incorporate several volatile elements in its crystal structure [Ca₅(PO₄)₃(Cl, F, OH)]. Apatite has been used to investigate a range of processes, including the pre-eruptive melt volatile concentrations (Scott et al. 2015); estimate volatile budgets and ascent rates (Li et al. 2021; Bernard et al. 2022); and reveal the relationship between volatile dynamics and an eruption's explosivity (Stock et al. 2016); the record of metasomatism and fluid movement in the crust (Harlov 2015); track arc magma evolution (Nathwani et al. 2020); record carbonate assimilation into magma (Wang et al. 2014), or to investigate estimate lunar water contents (Boyce et al. 2014).

Here we use apatite to constrain magma ascent rates from three historical eruptions at Kelud (Kelut) volcano in East Java, Indonesia. Kelud has produced some of the most hazardous historical eruptions in Southeast Asia (Goode et al. 2018; Maeno et al. 2019; Utami et al. 2022). Hazards from explosive eruptions of Kelud are capable of wreaking havoc across the island of Java. For example, the 2014 eruption produced voluminous ash and ejecta that severely impacted several cities. Hazards included building collapse from tephra accumulation, destructive lahars, extensive building damage, and ash inhalation (Blake et al. 2015; Williams et al. 2020). On the other hand, dome emplacements such as in 2007 were less immediately hazardous but may have enabled accumulation of excess volatiles underneath the dome, and thus increase the explosivity potential of the next eruption (Utami et al. 2022). In this study we characterized the volatile chemistry of apatite crystals from the 1990, 2007, and 2014 explosive and effusive eruptions of Kelud. We also modelled the zoning of Cl, F, and OH in apatite crystals to obtain timescales of magma ascent and degassing. From this parameter we calculated magma ascent rates that, combined with magma viscosities, are used to infer the possible magmatic processes associated with each eruptive style.

79

## Kelud volcano and its historical eruptions

Kelud (Kelut) is an active volcano in East Java, Indonesia classed as a basaltic andesite semi-plugged strato-cone (Whelley et al. 2015) that has produced >30 explosive eruptions and effusive dome extrusions since 1000 CE, with a mean recurrence interval of ~23 years (Global Volcanism Program 2013). The volcano is permanently monitored with a dedicated
observatory (Global Volcanism Program 2013) by the Centre of Volcanology and Geological Hazards Mitigation of Indonesia (CVGHM, or PVMBG in Bahasa Indonesia). The explosive eruptions in 1990 and 2014 were short-lived, with unrest signals escalating within hours from the onset of unrest to eruption, making them challenging to forecast (Nakamichi et al. 2017). In contrast, the 2007 dome extruded effusively after a month of unrest, which included elevated seismicity and carbon dioxide emissions into the lake (Hidayati et al. 2009; Caudron et al. 2012). Unrest signals for all three eruptions consisted of heightened seismicity, deformation, and changes of the lake temperature above background levels (Lesage and Surono 1995; Vandemeulebrouck et al. 2000; Hidayati et al. 2019) (see Supplementary Material). Both the 1990 and 2014 explosive events were preceded with an almost exponential increase in seismicity around 12 hours before each eruption (Lesage and Surono 1995; Hidayati et al. 2019), but this trend was not observed before the 2007 dome. Several short-lived unrest episodes or crises occurred after the 1990 eruption and until the 2014 eruptions (Global Volcanism Program 1991) including 1996 and 2001, but they did not result in an eruption (Global Volcanism Program 1996; 2001).

Previous studies on timescales of magmatic processes at Kelud used the Mg-Fe interdiffusion in reversely zoned pyroxenes before the 1990 eruption (Utami et al. 2021), which they compared with timings of seismic and hydroacoustic unrest and run-up times (Lesage and Surono 1995; Vandemeulebrouck et al. 2000). Timescales of Mg-Fe interdiffusion in pyroxenes span weeks to months, and were attributed to degassing and fluid accumulation. These timescales also overlapped with the onset of seismic and hydro-acoustic unrest that began around three months before the eruption. The exponential run-up in seismicity that began ~12 hours before the onset of the 1990 eruption gave magma ascent rates of ~0.2 m s⁻¹ based on magma storage depths of 6-7 km (Lesage and Surono 1995; Utami et al. 2021, 2022). Other estimates of magma ascent rates for the more recent 2014 eruption are around 4-5 m s⁻¹ (Cassidy et al. 2019), which were calculated using shallower storage depths of 2-3 km.

## Methodology

We studied apatite crystals from one main-stage and one late-stage pumices of the 1990 deposits, one fragment from the 2007 lava dome, and two pumices from the sub-Plinian phase of the 2014 eruption. These are the same samples as in Utami et al. (2021) and Humaida
(2013), where more details about sample locations and descriptions can be found.

Chemical Analyses of Apatite

We obtained backscattered electron (BSE) images of carbon-coated thin sections for each sample, and identified apatite based on their Electron Dispersive Spectra (EDS) using a JEOL-JSM 7800F scanning electron microscope (SEM) at the Earth Observatory of Singapore (EOS). We then mapped the apatite grains for Mg, Fe, Cl, F, and S in Wavelength Dispersive Spectra (WDS) mode on the JEOL JXA-8530F field emission electron microscope (EPMA) at the Facility for Analysis Characterisation Testing and Simulation (FACTS-NTU), Nanyang Technological University, Singapore (NTU). Guided by these elemental maps, we selected target crystals for EPMA analyses.

Target apatite crystals were analyzed in spot mode broadly following the methodology outlined in Li et al. (2021), using a 15kV beam and 10 nA current, but with 4 µm beam size. Profiles of zoned apatite crystals from the 1990, 2007, and 2014 eruptions were analyzed using a 10 nA current and 1-4 µm beam size. The structural formula of each apatite was calculated using approach #1 of Ketcham (2015), which was adapted based on an estimate of the [OH] content. We used a threshold of < 5% deviation from the ideal cation ratio of 5/3 to discard non-stoichiometric apatite analyses and/or mixtures. The mole fraction of volatile elements in the M2 column anion sites are given as X<sub>Cl</sub>, X<sub>F</sub>, and X<sub>OH</sub>, respectively. The OH cation was calculated by subtracting the atom per unit formula of Cl and F from the site total and normalizing the value according to:

$$X_{OH} = 1/2 \times [2 - (Cl + F)]. \qquad (1)$$

Here we refer to the calculated OH as OH<sub>c</sub>. Relative errors from the EPMA analyses are 1-3% for X<sub>Cl</sub>, 1-4% for X<sub>F</sub>, and the propagated errors are 1-33% for the calculated X<sub>OH</sub>. Some of the compositional profiles in the 1990 apatite crystals had major element oxide totals that were < 98 wt. %, perhaps due to thicker carbon coating (> 15 nm) present in some thin sections. We determined the orientation of the crystallographic axes of apatite (α and γ for a and c-axes, respectively) using the Electron Back-Scatter Diffraction (EBSD) mode in the SEM at EOS. These crystallographic orientations were used to calculate the angle between each EPMA profile and the crystallographic axes needed for diffusion coefficient correction (Costa et al. 2008), with an in-house EBSD conversion code employing MTEX Texture Analyses Toolbox.
151

152

Modelling timescales and ascent rates from apatite

153

154

155

156

157

158

159

160

We modelled the concentrations of Cl, F, and OHc in apatite assuming volume diffusion through the crystal by using the ApTimer application (Li et al. 2020). Due to the large anisotropy of diffusion for many elements in apatite, we used the acute angles between the profiles and the crystallographic a- and c-axes (α and γ, respectively) to correct for Cl, F, and OH diffusion coefficients from Li et al. (2020). The anisotropy is such that the diffusion parallel to the c-axis is considerably faster than parallel to the a-axes. Thus a correction factor was applied to the diffusion coefficients (Dᵢ) of three profiles that were parallel to the a-axes:

$$D_{Cl \parallel a} = 166 \times D_{Cl \parallel c}, \text{ where } D_{Cl \parallel a} = D_{Cl \parallel b} \qquad (2)$$

$$D_{F \parallel a} = 19 \times D_{F \parallel c}, \text{ where } D_{F \parallel a} = D_{F \parallel b} \qquad (3)$$

$$D_{OH \parallel a} = 46 \times D_{OH \parallel c}, \text{ where } D_{OH \parallel a} = D_{OH \parallel b}. \qquad (4)$$

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

These factors were determined by first calculating the diffusion coefficients at 950°C for Cl, F, and OH parallel to the c-axis, and then dividing their values by diffusion coefficients parallel to the a-axis at the same temperature. This was necessary because Li et al. (2020) only determined the diffusion coefficients of these volatile elements parallel to a-axis at 950°C. Here we assume that the activation energy for these elements does not change with crystallographic orientation, which is the case for some minerals and elements (e.g. Fe-Mg in olivine; Dohmen and Chakraborty 2002). We mostly used the Cl profile to obtain best fit timescales as it was the most prominently zoned and precisely measured element, but in a few cases, we used the F profile instead. The initial boundary conditions were chosen by visual inspection of the 1-D profile shape and used to determine the most suitable core and rim compositions. A constant rim composition was used, and diffusion towards the core generated the observed concentration gradients (Li et al. 2020). This is a simplification of the diffusion process, as the composition of the apatite crystal rim likely changed in a complex manner during magma ascent due to loss of volatiles and microlite crystallization. Using compositions that change with time at the boundary would give comparatively shorter times (Costa et al. 2008). We used temperatures of 974-980 ± 39°C for the three eruptions, as constrained through geothermobarometry and phase equilibria results reported in Utami et al. (2022). All inputs into the ApTimer application are listed in Table 1.

For unzoned apatite from the 1990 and 2014 eruptions, we obtained the maximum time
of diffusion using the Cl concentration on the condition that if there was a given diffusion length (see below), we would have been able to detect it in our elemental Cl maps and EPMA profiles. We thus used a simplified diffusion equation for the penetration depth from the crystal rim towards the interior to calculate time (Costa et al. 2008):

$$x = \sqrt{4D_{cl}t}. \qquad (5)$$

Where x is the diffusion distance, which we chose to be 4 µm. D<sub>Cl</sub> is the diffusion coefficient for Cl from Li et al. (2020); no correction factor necessary as several grains we investigated had profiles that were parallel to the c-axis. Magma ascent rates were calculated using the timescales and magma storage depths of 7 km, or storage pressures of 175 MPa derived from phase equilibria experiments and geothermobarometry (Utami et al. 2021, 2022), along with a mean crustal density of 2640 kg m<sup>-2</sup> (Jeffery et al. 2013).

191

## 192 Magma viscosity calculations

193

To quantitatively compare the magma physical properties between the different eruptions, we first calculated the melt viscosity using the formulation of Giordano et al. (2008). For this, we used the calculated matrix glass composition calculated from Utami et al. (2022) for a range of pre-eruptive temperatures (average T = 974-980 ± 39°C) and the same pre-eruptive water content in the melt as our starting point (4 wt. %), which was chosen based on previous studies on the same samples (Utami et al. 2021, 2022). To incorporate the effect of crystals, we calculated the effective viscosity for all three eruptions using the Einstein-Roscoe equation (e.g. Einstein 1906; Roscoe 1952; Gonnermann and Manga 2013):

$$\eta_{eff} = \eta_{melt} \times (1 - \frac{\phi}{\omega})^{-2.5}. \qquad (6)$$

Where η<sub>eff</sub> is the effective viscosity, η<sub>melt</sub> is the melt viscosity, φ is the volume fraction of crystals, and ω is the threshold volume fraction of crystals required for to transition to a rigid system, for which we used the value of 0.6. The phenocryst content of the 1990 and 2014 eruption deposits range from 28 to 47 vol. %, which allowed us to use the Einstein-Roscoe equation. The phenocryst content of the 2007 dome are up to 70 vol. %, which precludes accurate estimation of the effective viscosity for the dome rocks with Eq. 6. We thus compared several viscosity models capable of calculating effective viscosities for volume fraction up to 0.8 (Caricchi et al. 2007; Costa et al. 2009) and followed a conceptually accurate model (Faroughi and Huber 2015; Moitra and Gonnermann 2015; Moitra et al. 2018). These models

210
also require the ascent rate and conduit diameter, which are a priori unknown. However, the mass discharge rate (MDR) for the three eruptions were calculated by Maeno et al. (2019), and these can be related to the exit velocity ($v_e$) in m s⁻¹, magma density at discharge ($\rho_{discharge}$) in kg m⁻³, and conduit cross-sectional area ($A_{conduit}$) in m² through the following equation:

$$ MDR = \rho_{discharge} \times v_e \times A_{conduit}. \qquad (7) $$

The magma exit velocity is in principle different and higher than the magma ascent velocity we have calculated from the apatite diffusion models, and the magma density at pre-eruptive conditions is also different and higher than that at exit, when magma can be full of gas bubbles due to volatile exsolution. However, for the 2007 dome eruption the magma ascent velocity and the density we calculate at depth are very similar to the values at exit, and we have used them to obtain a range of conduit diameters and thus calculated the viscosity of the 2007 magma. For the 1990 and 2014 eruptions we used the Einstein-Roscoe equation to calculate their respective viscosities. Moreover, as cross-check we used the MDR values and measured density (1400 kgm⁻³) for the 1990 and 2014 eruptions reported by Maeno et al (2019), as well as Fig. 4.10 from Gonnerman and Manga (2013) to estimate the ascent velocity, from which we also obtained a range of conduit diameters. These values, together with the ascent rates from the apatite, were used to calculate the magma viscosity of the 1990 and 2014 magmas using various viscosity models that would better account for the crystal content (Caricchi et al. 2007; Costa et al. 2009; Faroughi and Huber 2015; Moitra and Gonnermann 2015; Moitra et al. 2018). For the Caricchi et al. (2007) and the Costa et al. (2009) viscosity models, the strain rate for the 2007 dome was needed. Thus the strain rate was calculated by dividing the magma ascent rate with the conduit radius (Moitra and Gonnermann 2015; Moitra et al. 2018).

## Results

Apatite crystals from the 1990 main-stage pumices are found in the groundmass (Fig. 1a), and as inclusions in pyroxenes and partially reacted amphibole phenocrysts (Fig. 1b, 1c). These apatite crystals are typically subhedral with rounded edges to euhedral, with sizes of around 100-500 µm. We find apatite in similar textural positions in the late-stage pumice, except for the apatite found in the quartz-bearing enclaves, which could have previously been groundmass apatite (Fig. 1d, 1f). Some apatite crystals are found in close textural association with magnetite and ilmenite or are found attached to phenocrysts. The 2014 eruption apatite has similar textural relations as that of the 1990 main-stage pumices (Fig. 1j, 1l), but it is also found inside
242

melt inclusions. Groundmass apatite crystals in the 2007 dome are significantly larger (< 1000 µm) than either of the 1990 and 2014 eruption deposits (Fig. 1g, 1h). Apatite in the 2007 dome rocks is also subhedral to euhedral and found as inclusions in melt embayments in glomerocrysts (Fig. 1g) or wedged between phenocrysts (Fig. 1h).

243

244

245

246

247

Apatite compositions from the 1990 and 2014 deposits mostly overlap, but contain some chemically distinct populations. For example, single 2007 dome apatite crystals mostly overlap with those from the other eruptions in F (1.4-1.7 wt. %) and Cl (0.6-1.0 wt. %) with only two grains recording lower Cl (0.6-0.8 wt. %). The 2007 dome apatite crystals in contact with the melt, for example those found in the groundmass or wedged in between phenocrysts, show reverse apatite zoning with rims are generally slightly higher in Cl than in cores (Fig 2b). This is confirmed in X-ray elemental maps, where we see that Cl increases from core to rim (Fig 2c, e). Apatite not in contact with the melt is unzoned (Fig. 2d).

248

249

250

251

252

253

254

255

256

The 1990 main-stage apatite volatile contents largely overlap irrespective of their textural positions (groundmass or inclusions), with 0.9-1.6 wt. % F, and 0.6-0.9 wt. % Cl. Apatite in the late-stage pumice is, however, enriched in F (2.0-2.4 wt.%), and in Cl (1.1-1.3 wt. %) compared to the main-stage apatite (Fig 3a). This is consistent with previous evidence of the late-stage interstitial melt being more evolved than in the main-stage pumice (Utami et al. 2021). The 2014 apatite volatile compositions overlap with those of the 1990 main-stage apatite in F (1.0-1.7 wt. %) and Cl (0.6-1.7 wt. %). The exception are amphibole-hosted apatite crystals that have significantly higher Cl concentrations (1.7 wt. %) compared to the rest (0.6-1.0 wt. %).

257

258

259

260

261

262

263

264

265

266

None of the apatite crystals from the 1990 nor the 2014 pumices are visibly zoned in volatile elements with the range of analytical techniques we have used (Fig. 3b-e). The 2007 dome apatite crystals are different from the 1990 and 2014 eruption in that they are reversely zoned in Cl, F, and OH<sub>c</sub> over lengths of about 9-25 µm. Rim compositions of zoned apatite are 0.9-3.2 wt. % F and 0.5-1.6 wt. % Cl, whereas core compositions have 1.4-1.7 wt. % F and 0.6-0.9 wt. % Cl. We identified two main types of zoning in the 2007 apatite crystals. One common feature in the groundmass apatite is reverse zoning in Cl, which increases towards the rims (reverse zoning) (Fig. 2). F and OH<sub>c</sub> are mirror images in concentration of those of Cl from core to rim (Fig 4a-d). Other apatite crystals show a reverse zoning with an increase in F towards the rims (Fig. 4e-h), alongside a slight increase in Cl, whereas OH<sub>c</sub> decreases. This

267

268

269

270

271

272

273

274

275
276

277

278

279

280

281

# Discussion

282

## Interpretation of the volatile zoning patterns in apatite

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

zoning is found in crystals that are not exposed to the melt or only partly exposed (e.g. wedged between phenocrysts, melt tubes/re-entrants, or unsealed melt pockets in glomerocrysts). Only the apatite crystals from the 1990 late-stage pumice, and the amphibole-hosted apatite from 2014 display significantly different compositions.

We propose that the zoning of Cl, F, and OH in apatite during magma ascent can be related to the combined effect of magma degassing (e.g. water loss from the melt) during ascent and the associated microlite crystallization (e.g. Blundy and Cashman 2000; Couch et al. 2003; Martel, 2012). However, due to the coupled partitioning and uphill diffusion of OH, Cl, and F into the same crystallographic site of apatite (Li et al. 2020), their concentration gradients are not straightforward to interpret. Water loss from the volatile exsolution typically leads to crystallization of microlites of minerals including plagioclase, pyroxenes, and Fe-Ti oxides (e.g. Blundy and Cashman 2000; Toramaru et al. 2008; Bernard and Bouvet de Maisonneuve 2020). Given that these minerals do not incorporate significant amounts of Cl and F, we propose that the concentrations of these volatiles would increase in the melt, and thus they would also increase in apatite. The amount of increase may depend on several interrelated factors, including the value of the partitioning coefficient of volatile elements between apatite and the surrounding melt, and their dependence on intensive melt parameters such as temperature (Li and Costa, 2019). The relationship could also change if there is pressure- and temperature-dependent partitioning of these elements between the melt and apatite in the presence of an exsolved fluid phase or brine (Li et al., 2021). Despite this range of possible scenarios, we hypothesize that the volatile zoning in Kelud apatite is likely driven by coupled magmatic-volatile processes, including degassing, microlite crystallization, and in the case of the 2007 zoned dome apatite, diffusion from slow ascent. Crystals with F reverse zoning found as inclusion in other minerals, or only partly open to the melt may behave similarly to volatile trends in melt inclusion embayments (e.g. Myers et al. 2019; Myers et al. 2021; Saalfeld et al. 2022). In this case, the apatite volatile zoning may not be fully controlled by degassing-induced crystallization. Rather, the zoning could reflect a localized temperature decrease, as F partitioning between apatite and melt is more sensitive towards temperature dependence than
308 both Cl and OH (Li and Costa, 2019).

309

310 Other studies of zoned apatite have shed light on the role of volatile on pre-eruptive
311 processes. For example, Li et al. (2021) compared the volatile budget of zoned apatite from the
312 2006 versus unzoned apatite from 2010 rocks of Merapi volcano, and proposed that input of
313 deep volatile into shallow magma column before the 2010 eruption led to fast magma ascent
314 rates, and thus higher degree of explosivity. At other magmatic systems including Campi
315 Flegrei and Laacher See, volatile zoning in apatite are interpreted to reflect the progressive
316 volatile saturation in the system as a result of fractional crystallization, and re-organization of
317 the magmatic system (Stock et al. 2018; Humphreys et al. 2021).

318

319 Apatite timescales and their relation to volcanic unrest

320

321 2007 Dome eruption

322

323 A total of 17 diffusion timescales were obtained from 12 different crystals of the 2007 dome.
324 Best fit timescales of up to 175 days were found, although 71% of them are ≤ 60 days (Table
325 1). Most of the zoned apatite crystals from the 2007 eruption were in contact with the melt.
326 Moreover, three profiles have anomalously short timescales. For example, one zoned apatite
327 crystal was included in an amphibole pseudomorph gave a timescale of 8 days, which may
328 record amphibole breakdown. Another profile was short (4 µm) and yielded a model time of 1
329 day. Although these timescales may be related to late-stage processes, we propose instead that
330 because they are apatite inclusions in minerals, their exposure to the melt is limited compared
331 to those fully exposed in the melt. These observations supports our assumption of an initial
332 homogeneous composition, and the hypothesis that the zoning we observe in apatite is related
333 to ascent and associated changes in the melt.

334

335 Apatite compositional profiles closer to the crystallographic *a*-axis systematically record
336 timescales that are at least 2.5 times longer (63-175 days) than those closer to the *c*-axis (1-58
337 days). This difference is most noticeable where multiple profiles were measured on the same
338 crystal. Such discrepancy could stem from a paucity in the knowledge and limited data on
339 diffusion coefficients of volatile elements along the *a*-axis. We further observe no systematic
340 difference between the times obtained from crystals that are mainly reversely zoned in Cl vs.
341

those that are mainly zoned in F. Due to the potential kinetic delays between magma degassing, microlite crystallization and reaction of the apatite's F-Cl-OH concentrations to these processes, the calculated times are likely to be the minima especially for crystals that are not directly in contact with the external melt.

342

343

344

345

346

In relating the apatite timescales to pre-eruptive unrest signals observed before the 2007, we interpret the cluster of timescales < 30 days as the final stage of magma ascent that began in early October and continued until the eruption date (see Supplementary Material). These dates overlap with the onset of different types of unrest in October 2007, such as seismic swarms, increases in seismic energy released, and inflation/deflation (Hidayati et al. 2019). The relationships between the timing from the apatite zoning and the monitoring data thus indicate magma movement and ascent towards the surface prior to the dome emplacement. Moreover, timescales > 60 days belong to profiles that were measured along the *a*-axis. Given the limited experimental constraints on the diffusion coefficients along this axis, we are cautious about directly correlating these timescales with any unrest signal and ascribing them to magma ascent.

347

348

349

350

351

352

353

354

355

356

357

358

The 1990 and 2014 eruptions

359

360

The apatite crystals from the 1990 and 2014 we observed (six in total) are unzoned in Cl, F, and OH at the scale of the EPMA maps and profiles. This lack of zoning likely reflects fast magma ascent and thus limited transport time to be able to develop measurable chemical zoning. Notwithstanding, we can only calculate the time that would be required to generate an observable zoning length with the analytical techniques we have used. We found that using a length of 4 µm implies a maximum time of diffusion of 7-8 hours, which means that the magmas would have been transported from their reservoir to the surface in less than this time. These maximum timescales are shorter than the 1990 and 2014 total unrest (1-3 months - Lesage and Surono 1995; Hidayati et al. 2019; Utami et al. 2021) but are of the same order of magnitude as the run-up times or exponential increase in seismic unrest directly before the eruptions (12-48 hours; Lesage and Surono 1995; Hidayati et al. 2019). The relationship between the maximum timescales and the duration of unrest suggests that for the 1990 and 2014 eruptions, degassing and magma ascent occurred within 12-48 hours before eruption. Such short timescales of ascent likely account for the fast escalation of seismicity before

361

362

363

364

365

366

367

368

369

370

371

372

373
374 eruption.

375

376 **Magma ascent rates for explosive and dome eruptions**

377

378 We can estimate the magma ascent rates of Kelud by combining the timescales from apatite
379 and the magma storage pressure estimated from phase equilibria experiments, mineral
380 geobarometers, and volatile saturation pressures (Utami et al. 2021, 2022). The pressure was
381 converted into depths below the edifice using an averaged density of the crust of 2640 kg m⁻³,
382 which was previously proposed for the upper crust under Kelud (Jeffery et al. 2013). Utami et
383 al. (2021; 2022) reported a range of storage pressures from 160 to > 400 MPa depending on
384 the minerals, which suggests a melt and mineral distribution across a range of depths
385 underneath Kelud, with most crystals around 100-200 MPa. Here we chose an intermediate
386 pressure of 175 MPa (150-200 MPa), corresponding to a depth of ~ 6-7 km below the crater,
387 that best matches the phase assemblage and phenocryst content of the three magmas according
388 to phase equilibria experiments (Utami et al. 2022).
389

390

391 We find that the calculated ascent rates in this manner vary between 0.4-78 × 10⁻³ m s⁻¹
392 but most are between 2.6-39 × 10⁻³ m s⁻¹ for the 2007 dome (Table 1). These ascent rates from
393 the apatite crystals are comparable or faster than the extrusive rates from other dome building
394 eruptions such as Mt St Helens and Soufriere Hills (1-8 × 10⁻³ m s⁻¹) (Devine et al. 1998). Our
395 calculated ascent rates are only estimates of their mean values, as magma ascent rate changes
396 depending on processes occurring as the magma moves towards the surface, and the magma
397 storage depth is for the 2007 is less well defined than for the 1990 and 2014 eruptions (Utami
398 et al. 2022).
399

400

401 For the 1990 and 2014 magmas we only obtained minimum magma ascent rates of about
402 0.26 m s⁻¹ from the maximum transfer times as discussed in the previous section. These are
403 more than an order of magnitude higher than those of the 2007 dome magma. Such a difference
404 in ascent rates would account for the contrasting eruption styles, where the magmas feeding
405 explosive eruptions would be able to retain their dissolved volatile content, but the slow ascent
406 would allow the dome to release the volatiles through permeable degassing (Gonnermann and
407 Manga, 2013). These ascent rates are significantly lower than the rates of 4-5 m s⁻¹ previously
408 proposed for the 2014 eruption (Cassidy et al. 2019), but we note that the difference in ascent
409
407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

rates also depends on the chosen depths as a starting point. Still, our estimates are still comparable with the range of ascent rates for explosive eruptions from volcanoes including Mt St Helens and Soufrière Hills (Rutherford, 2008). We find that ascent rates for explosive eruptions of Kelud overlaps with rates from other historical eruptions from basaltic andesite volcanoes, such as 1974 eruption of Volcan de Fuego, and the 1986 eruption of Izu-Oshima (Cassidy et al. 2018).

## Relating ascent rates to eruptive styles

Many volcanoes show a range of eruptive styles, and these are likely controlled by magmatic processes resulting in rapid changes in pre-eruptive and syn-eruptive parameters (e.g. Cassidy et al. 2018; Annen and Burgisser 2021; Annen et al. 2023). Ultimately, magma fragmentation is involved in the explosivity of a given eruption (Dingwell et al. 1993), a process that is highly dependent on the physical properties of the magma and fluids, including the effective viscosity. We hereby attempt to relate the ascent rates to the effective viscosities to investigate the dominant processes controlling explosive and effusive eruptions at Kelud volcano, as proposed by Gonnermann and Manga (2013).

The magma viscosities of 1990 and 2014 explosive events are quite different from the 2007 dome, mainly because of the much higher phenocryst content of the latter (Table 3; Fig. 6). The magma viscosity models of Faroughi and Huber (2015), Costa et al. (2009), and Caricchi et al. (2007) yield comparable effective viscosity ranges for the 1990 ($10^{3.8}-10^{7.5}$ Pa s) and 2014 ($10^{3.6}-10^{7.2}$ Pa s) eruptions for a minimum magma ascent rate of $> 1 \times 10^{-1}$ m s$^{-1}$, respectively. Magma viscosities calculated using the viscosity model of Moitra et al. (2015; 2018) for the 1990 ($10^{5.1}-10^{8.6}$ Pa s) and 2014 ($10^{4.8}-10^{8.3}$ Pa s) give slightly higher values compared to the previous models, but still comparable with those calculated using the Einstein-Roscoe equation for the same eruptions (1990 => $10^{4.2}-10^{7.8}$ Pa s; 2014 => $10^{4.0}-10^{7.5}$ Pa s). In contrast, the 2007 dome span a wide range of effective viscosities for the given magma ascent rate ($10^{-4}$ to $10^{-2}$ m s$^{-1}$). The magma viscosity models of Faroughi and Huber (2015), Costa et al. (2009), and Caricchi et al. (2007) yield effective viscosities for the 2007 dome of $10^{7.8}-10^{12.5}$ Pa s. The viscosity for the 2007 dome ($10^{6.4}-10^{10.0}$ Pa s) using the model of Moitra et al. (2015; 2018) give slightly lower effective viscosity ranges relative to other models.
440

We used the relationship proposed in Gonnermann and Manga (2013) to connect the magma viscosity and magma ascent rates to eruptive styles (Fig. 6). We find that the initial processes dominating the 1990 and 2014 explosive eruptions are inertial fragmentation and bubble nucleation, straddling the Hawaiian and sub-Plinian fields. Given that both the 1990 and 2014 eruptions were sub-Plinian (e.g. Maeno et al. 2019), they eventually evolved towards sub-Plinian field where shear localization, viscous heating, and brittle fragmentation dominate (Gonnermann and Manga 2013). These processes contrast with the dominant process for the 2007 dome, which is primarily permeable outgassing. The data shows some overlap with the sub-Plinian field, but the combination of slow ascent and outgassing would make the 2007 magma so degassed that processes such as brittle fragmentation were unlikely to happen.

441

442

443

444

445

446

447

448

449

450

451

The ascent rates we calculated for the 1990 and 2014 events imply mass discharge rates >10<sup>5.5</sup> kg s<sup>-1</sup>, which overlaps with discharge rates (10<sup>7</sup>-10<sup>8</sup> kg s<sup>-1</sup>) determined for the same eruptions by Maeno et al. (2019). Combining the magma viscosity with the discharge rate for the 2007 dome (5.8 × 10<sup>3</sup> - 1.9 × 10<sup>4</sup> kg s<sup>-1</sup>) (Maeno et al. 2019), we find that the corresponding velocity becomes slightly higher at around 10<sup>-1.5</sup> m s<sup>-1</sup>. This velocity overlaps with the upper range of the ascent rates determined from apatite, but as the MDR is estimated at existing conditions, the corresponding exit velocities of the 1990 and 2014 eruptions could be significantly higher than the one we estimated. Moreover, these relationships do not typically consider the presence of an exsolved volatile phase at pre-eruptive conditions as we have proposed for the 1990 and 2014 eruptions (Utami et al. 2022), and which should have an influence on the viscosity, magma ascent rates, and eruptive style relations (Cassidy et al. 2018). Notwithstanding these limitations, we find it remarkable that magma ascent rates derived from apatite are consistent with other parameters such as discharge rates, which potentially enables us to characterize different eruption styles.

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

Conclusions

467

468

469

470

471

472

The duration and rates of magma ascent are considered to be key parameters that control the eruption style. Here we show that using the volatile zoning (F-Cl-OH) in apatite from the last three eruptions of Kelud volcano (1990, 2007, 2014), it is possible to gain insights into explosive and effusive eruptions at Kelud volcano, Indonesia using apatite. Unzoned apatite from the explosive eruptions of 1990 and 2014 record evidence for minimum ascent rates of >
2.6 × 10⁻¹ m s⁻¹. In contrast, reversely zoned apatite from the 2007 dome rocks record evidence for considerably slower ascent rates (> 0.4 × 10⁻³ m s⁻¹). The estimated ascent rates for the 1990 explosive eruptions overlap with magma ascent rates constrained through the onset of exponential increase in seismicity directly prior to the eruption. Our ascent rates are consistent with ascent rates from other volcanoes with effusive and explosive eruptions (e.g. Mount St Helens, Soufriere Hills, Soufriere St. Vincent), and thus we can establish a lower bound on the ascent rates of explosive eruptions for Kelud and other analogous volcanoes. The reversely zoned apatite from the 2007 dome allows us to estimate timescales of buoyancy driven magma ascent towards dome emplacement, which overlaps with the seismic and geodetic unrest. Viscosity and magma ascent relationships show that whereas the 1990 and 2014 eruptions are typical sub-Plinian eruptions with dominating processes including brittle fragmentation, bubble nucleation, and shear localization, dome emplacement is likely dominated by permeable outgassing. Our study demonstrates the utility of apatite to shed light on variations in eruption style for volatile dominant processes at active volcanoes such as Kelud to record timescales of ascent and rates that are consistent with eruptive style at active volcanoes.

## Author Contributions

SBU: Conceptualization, investigation, project administration, methodology, formal analyses, data curation, visualization, writing (original draft; review and editing). FC: Supervision, conceptualization, investigation, resources, methodology, writing – review & editing, supervision, project administration, funding acquisition, final approval of manuscript. HH: Writing – review & editing, final approval of manuscript.

## Acknowledgements

We thank S.J. Lim for assistance with sample preparation; J. Oalmann for assistance with the EBSD analyses of apatite and providing the code to calculate crystallographic angles from EBSD data; P. Moitra and M. Pistone for discussions on magma viscosities and assistance with viscosity model calculations; W. Li and C. Widiwijayanti for discussions on apatite, ApTimer, and Kelud volcano, respectively. We would like to acknowledge FACTS-NTU, Singapore, for the use of their analytical facilities. This research received support from the Earth Observatory
of Singapore via its funding from the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative. This project was funded by a Singapore National Research Foundation Investigator Grant (NRF-NRFI2017-06-00) awarded to FC. This work is part of Sri Budhi Utami's PhD thesis. The authors hereby declare no conflict of interest.

510

## 511 References

512

513 Albert H, Costa F, Marti J (2015) ‘Timing of Magmatic Processes and Unrest Associated with Mafic Historical Monogenetic Eruptions in Tenerife Island’. *Journal of Petrology* 56 (10): 1945–66.

514
515 Andújar J, Costa F, Scaillet B (2013) ‘Storage Conditions and Eruptive Dynamics of Central versus Flank Eruptions in Volcanic Islands: The Case of Tenerife (Canary Islands, Spain)’. *Journal of Volcanology and Geothermal Research* 260: 62–79.
516
517
518
519 https://doi.org/10.1016/j.jvolgeores.2013.05.004.

520 Annen C, Burgisser A (2021) ‘Modeling Water Exsolution from a Growing and Solidifying Felsic Magma Body’. *Lithos*, no. 402: 105799–105799.
521
522 https://doi.org/10.1016/j.lithos.2020.105799.

523 Annen C, Pinel V, Burgisser, A (2023) ‘Re-Evaluation of the Role of Volatiles in the Rupture of Magma Chambers and the Triggering of Crystal-Rich Eruptions’. *Journal of Volcanology and Geothermal Research*, January, 107755.
524
525
526 https://doi.org/10.1016/j.jvolgeores.2023.107755.

527 Bernard O, Bouvet de Maisonneuve C (2020) ‘Controls on Eruption Style at Rabaul, Papua New Guinea – Insights from Microlites, Porosity and Permeability Measurements’. *Journal of Volcanology and Geothermal Research* 406 (November): 107068.
528
529
530 https://doi.org/10.1016/j.jvolgeores.2020.107068.

531 Bernard O, Li W, Costa F, Saunders S, Itikarai I, Sindang M, Bouvet De Maisonneuve C (2022) ‘Explosive-Effusive-Explosive: The Role of Magma Ascent Rates and Paths in Modulating Caldera Eruptions’, 1–5. https://doi.org/10.1130/G50023.1.

532
533
534 Blake DM, Wilson G, Stewart C, Craig H, Hayes JL, Jenkins SF, Wilson TM, et al. (2015)
535
536 ‘Impacts of the 2014 Eruption of Kelud Volcano, Indonesia, on Infrastructure,
537
538
539 Utilities, Agriculture and Health’, no. November: 1–131.

537 Blundy J, Cashman K (2000) ‘Degassing and Crystallization of Ascending Andesite and Dacite’. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 358 (1770): 1487–1513.

540 Blundy J, Cashman K (2008) ‘Petrologic Reconstruction of Magmatic System Variables and Processes’. *Reviews in Mineralogy and Geochemistry* 69 (1): 179–239.

541
542 Boyce JW, Tomlinson SM, McCubbin FM, Greenwood JP, Treiman AH (2014) ‘The Lunar Apatite Paradox’. *Science* 344 (6182): 400–402.
543
544 https://doi.org/10.1126/science.1250398.

545 Bull, KF, Buurman H (2013) ‘An Overview of the 2009 Eruption of Redoubt Volcano, Alaska’. *Journal of Volcanology and Geothermal Research* 259 (March 2009): 2–15.
546
547 https://doi.org/10.1016/j.jvolgeores.2012.06.024.
548 Cadoux A, Scaillet B, Druitt TH, Deloule E (2014) ‘Magma Storage Conditions of Large
549 Plinian Eruptions of Santorini Volcano (Greece)’. *Journal of Petrology* 55 (6): 1129–
550 71. https://doi.org/10.1093/petrology/egu021.

551 Caricchi, L, Burlini L, Ulmer P, Gerya T, Vassalli M, Papale P (2007) ‘Non-Newtonian
552 Rheology of Crystal-Bearing Magmas and Implications for Magma Ascent
553 Dynamics’. *Earth and Planetary Science Letters* 264 (3–4): 402–19.
554 https://doi.org/10.1016/j.epsl.2007.09.032.

555 Cassidy M, Manga M, Cashman K, Bachmann O (2018) ‘Controls on Explosive-Effusive
556 Volcanic Eruption Styles’. *Nature Communications* 9 (1).
557 https://doi.org/10.1038/s41467-018-05293-3.

558 Cassidy M, Ebmeier SK, Helo C, Watt SFL, Caudron C, Odell A, Spaans K, et al (2019)
559 ‘Explosive Eruptions With Little Warning: Experimental Petrology and Volcano
560 Monitoring Observations From the 2014 Eruption of Kelud, Indonesia’.
561 *Geochemistry, Geophysics, Geosystems*, 1–30. https://doi.org/10.1029/2018gc008161.

562 Caudron C, Mazot A, Bernard A (2012) ‘Carbon Dioxide Dynamics in Kelud Volcanic
563 Lake’. *Journal of Geophysical Research: Solid Earth* 117 (B5).

564 Costa A, Caricchi L, Bagdassarov N (2009) ‘A Model for the Rheology of Particle-Bearing
565 Suspensions and Partially Molten Rocks’. *Geochemistry, Geophysics, Geosystems* 10
566 (3): 1–13. https://doi.org/10.1029/2008GC002138.

567 Costa F, Dohmen R, Chakraborty S (2008) ‘Time Scales of Magmatic Processes from
568 Modeling the Zoning Patterns of Crystals’. *Reviews in Mineralogy and Geochemistry*
569 69 (1): 545–94.

570 Costa F, Andreastuti S, Bouvet de Maisonneuve C, Pallister JS (2013) ‘Petrological Insights
571 into the Storage Conditions, and Magmatic Processes That Yielded the Centennial
572 2010 Merapi Explosive Eruption’. *Journal of Volcanology and Geothermal Research*
573 261: 209–35.

574 Couch S, Sparks RSJ, Carroll MR (2003) ‘The Kinetics of Degassing-Induced Crystallization
575 at Soufriere Hills Volcano, Montserrat’. *Journal of Petrology* 44 (8): 1477–1502.
576 https://doi.org/10.1093/petrology/44.8.1477.

577 Devine JD, Rutherford MJ, Gardner JE (1998) ‘Petrologic Determination of Ascent Rates for
578 the 1995-1997 Soufriere Hills Volcano Andesitic Magma’. *Geophysical Research
579 Letters* 25 (19): 3673–76. https://doi.org/10.1029/98GL00912.

580 Dingwell DB, Bagdassarov NS, Bussod GY, Webb SL (1993) ‘Magma Rheology’. In
581 *Experiments at High Pressure and Applications to the Earth’s Mantle: Short Course
582 Handbook*, edited by R.W. Luth, Volume 21, 131–96. Edmonton: Mineralogical
583 Association of Canada.

584 Dingwell DB (1996) ‘Volcanic Dilemma--Flow or Blow?’ *Science* 273 (5278): 1054–55.
585 https://doi.org/10.1126/science.273.5278.1054.

586 Dohmen R, Chakraborty S (2002) ‘Si and O Diffusion in Olivine and Implications for
587 Characterizing Plastic Flow in the Mantle’. *Geophysical Research Letters* 29 (21):
588 2030. https://doi.org/10.1029/2002GL015480.

589 Einstein A (1906) ‘Effect of suspended rigid spheres on viscosity’. *Annals of Physics* 19:
590 289–306.

591 Faroughi, SA, Huber C (2015) ‘A Generalized Equation for Rheology of Emulsions and
592 Suspensions of Deformable Particles Subjected to Simple Shear at Low Reynolds
593 Number’. *Rheologica Acta* 54 (2): 85–108. https://doi.org/10.1007/s00397-014-0825-8.

595 Georgeais G, Koga KT, Moussallam Y, Rose-Koga EF (2021) ‘Magma Decompression Rate
596 Calculations With EMBER: A User-Friendly Software to Model Diffusion of H2O,
597 CO2, and S in Melt Embayments'. *Geochemistry, Geophysics, Geosystems* 22 (7).

598 https://doi.org/10.1029/2020GC009542.

599 Giordano D, Russell JK, Dingwell DB (2008) 'Viscosity of Magmatic Liquids: A Model'.

600 *Earth and Planetary Science Letters* 271 (1): 123–34.

601 Global Volcanism Program, GVP (1991) 'Report on Kelut (Indonesia)'. Edited by L

602 McClelland. *Bulletin of the Global Volcanism Network* 16 (3).

603 https://doi.org/10.5479/si.GVP.BGVN199103-263280.

604 Global Volcanism Program, GVP (1996) 'Report on Kelut (Indonesia)'. Edited by R.

605 Wunderman 21 (11).

606 Global Volcanism Program, GVP (2001) 'Report on Kelut (Indonesia)'. Edited by R.

607 Wunderman. *Bulletin of the Global Volcanism Network* 26 (2).

608 Global Volcanism Program, GVP (2013) 'Kelut (Kelud)'.

609 https://doi.org/10.5479/si.GVP.VOTW4-2013.

610 Gonnermann, HM, Manga M (2013) 'Dynamics of Magma Ascent in the Volcanic Conduit'.

611 In *Modeling Volcanic Processes: The Physics and Mathematics of Volcanism*, edited

612 by S Fagents, T. K. P. Gregg, and R.M.C. Lopes, 9780521895:55–84. Cambridge

613 University Press. https://doi.org/10.1017/CBO9781139021562.004.

614 Goode LR, Handley HK, Cronin SJ, Abdurrachman M (2018) 'Insights into Eruption

615 Dynamics from the 2014 Pyroclastic Deposits of Kelut Volcano, Java, Indonesia, and

616 Implications for Future Hazards'. *Journal of Volcanology and Geothermal Research*.

617 https://doi.org/10.1016/J.JVOLGEORES.2018.02.005.

618 Harlov DE (2015) 'Apatite: A Fingerprint for Metasomatic Processes'. *Elements* 11 (3): 171–

619 76. https://doi.org/10.2113/gselements.11.3.171.

620 Hidayati S, Basuki A, Kristianto, Mulyana I (2009) 'Emergence of Lava Dome from the

621 Crater Lake of Kelud Volcano, East Java' 4 (4): 229–38.

622 Hidayati S, Triastuty H, Mulyana I, Adi S, Ishihara K, Basuki A, Kuswandarto H, Priyanto B,

623 Solikhin A (2019) 'Differences in the Seismicity Preceding the 2007 and 2014

624 Eruptions of Kelud Volcano, Indonesia'. *Journal of Volcanology and Geothermal*

625 Research 382: 50–67. https://doi.org/10.1016/j.jvolgeores.2018.10.017.

626 Humaida H (2013) 'A Geochemical Study on Merapi and Kelud Volcanic Eruptions.'.

627 Universitas Gadjah Mada, Indonesia. Unpublished PhD Thesis (in Bahasa Indonesia).

628 Humphreys MCS, Smith VC, Coumans JP, Riker JM, Stock MJ, de Hoog JCM, Brooker RA

629 (2021) 'Rapid Pre-Eruptive Mush Reorganisation and Atmospheric Volatile

630 Emissions from the 12.9 Ka Laacher See Eruption, Determined Using Apatite'. *Earth*

631 and Planetary Science Letters 576 (December): 117198–117198.

632 https://doi.org/10.1016/j.epsl.2021.117198.

633 Jeffery AJ, Gertisser R, Troll VR, Jolis EM, Dahren B, Harris C, Tindle AG, et al (2013)

634 'The Pre-Eruptive Magma Plumbing System of the 2007–2008 Dome-Forming

635 Eruption of Kelut Volcano, East Java, Indonesia'. *Contributions to Mineralogy and*

636 *Petrology* 166 (1): 275–308.

637 Ketcham, RA (2015) 'Technical Note: Calculation of Stoichiometry from EMP Data for

638 Apatite and Other Phases with Mixing on Monovalent Anion Sites'. *American*

639 *Mineralogist* 100 (7): 1620–23. https://doi.org/10.2138/am-2015-5171.

640 Kilgour, GN, Saunders KE, Blundy JD, Cashman KV, Scott BJ, Miller CA (2014)

641 'Timescales of Magmatic Processes at Ruapehu Volcano from Diffusion

642 Chronometry and Their Comparison to Monitoring Data'. *Journal of Volcanology and*

643 *Geothermal Research* 288: 62–75.

644 Lesage, P Surono (1995) 'Seismic Precursors of the February 10, 1990 Eruption of Kelut

645 Volcano, Java'. *Journal of Volcanology and Geothermal Research* 65 (1): 135–46.
646

647

648

649

Li W, Costa F (2019) ‘A Thermodynamic Model for F-Cl-OH Partitioning between Silicate Melts and Apatite Including Non-Ideal Mixing with Application to Constraining Melt Volatile Budgets’. *Geochimica et Cosmochimica Acta* 269: 203–22.
https://doi.org/10.1016/j.gca.2019.10.035.

650

651

652

Li W, Chakraborty S, Nagashima K, Costa F (2020) ‘Multicomponent Diffusion of F, Cl and OH in Apatite with Application to Magma Ascent Rates’. *Earth and Planetary Science Letters* 550: 116545–116545. https://doi.org/10.1016/j.epsl.2020.116545.

653

654

655

Li W, Costa F, Nagashima K (2021) ‘Apatite Crystals Reveal Melt Volatile Budgets and Magma Storage Depths at Merapi Volcano, Indonesia’. *Journal of Petrology* 62 (4).
https://doi.org/10.1093/petrology/egaa100.

656

657

658

Lubis, AM (2014) ‘Uplift of Kelud Volcano Prior to the November 2007 Eruption as Observed by L-Band Insar’. *Journal of Engineering and Technological Sciences* 46 (3): 245–57.

659

660

661

662

Maeno F, Nakada S, Yoshimoto M, Shimano T, Hokanishi N, Zaennudin A, Iguchi M (2019) ‘Eruption Pattern and a Long-Term Magma Discharge Rate over the Past 100 Years at Kelud Volcano, Indonesia’. *Journal of Disaster Research* 14 (1): 27–39.
https://doi.org/10.20965/JDR.2019.P0027.

663

664

Martel C, Schmidt BC (2003) ‘Decompression Experiments as an Insight into Ascent Rates of Silicic Magmas’, no. 144: 397–415. https://doi.org/10.1007/s00410-002-0404-3.

665

666

667

668

Martel, C (2012) ‘Eruption Dynamics Inferred from Microlite Crystallization Experiments: Application to Plinian and Dome-Forming Eruptions of Mt. Pelée (Martinique, Lesser Antilles)’. *Journal of Petrology* 53 (4): 699–725.
https://doi.org/10.1093/petrology/egr076.

669

670

671

Moitra P, Gonnermann HM (2015) ‘Effects of Crystal Shape- and Size-Modality on Magma Rheology’. *Geochemistry Geophysics Geosystems* 16: 1–26.
https://doi.org/10.1002/2014GC005554.

672

673

674

Moitra, P, Gonnermann HM, Houghton BG, and Tiwary CS (2018) ‘Fragmentation and Plinian Eruption of Crystallizing Basaltic Magma’. *Earth and Planetary Science Letters* 500: 97–104. https://doi.org/10.1016/j.epsl.2018.08.003.

675

676

677

678

Myers, ML., Wallace PJ, Wilson CJN (2019) ‘Inferring Magma Ascent Timescales and Reconstructing Conduit Processes in Explosive Rhyolitic Eruptions Using Diffusive Losses of Hydrogen from Melt Inclusions’. *Journal of Volcanology and Geothermal Research* 369 (January): 95–112. https://doi.org/10.1016/j.jvolgeores.2018.11.009.

679

680

681

682

Myers ML, Druitt TH, Schiavi F, Gurioli L, Flaherty T (2021) ‘Evolution of Magma Decompression and Discharge during a Plinian Event (Late Bronze-Age Eruption, Santorini) from Multiple Eruption-Intensity Proxies’. *Bulletin of Volcanology* 83 (3).
https://doi.org/10.1007/s00445-021-01438-3.

683

684

685

686

687

Nakamichi H, Iguchi M, Triastuty H, Hendrasto M, Mulyana I (2017) ‘Differences of Precursory Seismic Energy Release for the 2007 Effusive Dome-Forming and 2014 Plinian Eruptions at Kelud Volcano, Indonesia’. *Journal of Volcanology and Geothermal Research*, no. 382: 68–80.
https://doi.org/10.1016/j.jvolgeores.2017.08.004.

688

689

690

Nathwani, CL, Loader MA, Wilkinson JJ, Buret Y, Sievwright RH, Hollings P (2020) ‘Multi-Stage Arc Magma Evolution Recorded by Apatite In volcanic Rocks’. *Geology* 48 (4): 323–27. https://doi.org/10.1130/G46998.1.

691

692

693

Passarelli L, Brodsky EE (2012). ‘The Correlation between Run-up and Repose Times of Volcanic Eruptions’. *Geophysical Journal International* 188 (3): 1025–45.
https://doi.org/10.1111/j.1365-246X.2011.05298.x.

694

695

Preece K, Gertisser R, Barclay J, Berlo K, Herd R A, and Edinburgh Ion Microprobe Facility. (2014) ‘Pre- and Syn-Eruptive Degassing and Crystallisation Processes of the 2010
696

697

and 2006 Eruptions of Merapi Volcano, Indonesia'. *Contributions to Mineralogy and Petrology* 168 (4): 1-25. https://doi.org/10.1007/s00410-014-1061-z.

698

Puspito N T, Shimazaki K (1995) 'Mantle Structure and Seismotectonics of the Sunda and

699

Banda Arcs'. *Tectonophysics* 251 (1-4): 215-28. https://doi.org/10.1016/0040-1951(95)00063-1.

700

701

Putirka, KD (2008) ‘Thermometers and Barometers for Volcanic Systems’. *Reviews in Mineralogy and Geochemistry* 69 (1): 61–120.

702

703

Rasmussen DJ, Plank TA, Roman DC, Power JA, Bodnar RJ, Hauri EH (2018) ‘When Does Eruption Run-up Begin? Multidisciplinary Insight from the 1999 Eruption of Shishaldin Volcano’. *Earth and Planetary Science Letters* 486: 1–14.

704

705

706

https://doi.org/10.1016/j.eps1.2018.01.001.

707

Roman, DC, and Cashman KV (2018) ‘Top-Down Precursory Volcanic Seismicity: Implications for “Stealth” Magma Ascent and Long-Term Eruption Forecasting’.

708

709

*Frontiers in Earth Science* 6 (September): 1–18.

710

https://doi.org/10.3389/feart.2018.00124.

711

Roscoe, R (1952) ‘The Viscosity of Suspensions of Rigid Spheres’. *British Journal of Applied Physics* 3 (8): 267.

712

713

Rutherford, MJ (2003) ‘Magmatic Conditions and Magma Ascent as Indicated by Hornblende Phase Equilibria and Reactions in the 1995-2002 Soufriere Hills Magma’.

714

715

*Journal of Petrology* 44 (8): 1433–53. https://doi.org/10.1093/petrology/44.8.1433.

716

Rutherford, MJ (2008) ‘Magma Ascent Rates’. *Reviews in Mineralogy and Geochemistry* 69: 241–71. https://doi.org/10.2138/rmg.2008.69.7.

717

718

Saalfeld, MA, Myers ML, deGraffenried R, Shea T, Waelkens, CM (2022) ‘On the Rise: Using Reentrants to Extract Magma Ascent Rates in the Bandelier Tuff Caldera

719

720

Complex, New Mexico, USA’. *Bulletin of Volcanology* 84 (1).

721

https://doi.org/10.1007/s00445-021-01518-4.

722

Saunders K, Blundy J, Dohmen R, Cashman K (2012) ‘Linking Petrology and Seismology at an Active Volcano’. *Science* 336 (6084): 1023–27.

723

724

Scott JAJ, Humphreys MCS, Mather TA, Pyle DM, and Stock MJ (2015) ‘Insights into the Behaviour of S, F, and Cl at Santiaguito Volcano, Guatemala, from Apatite and Glass’. *Lithos* 232: 375–94. https://doi.org/10.1016/j.lithos.2015.07.004.

725

726

727

Stock, MJ, Humphreys MCS, Smith VC, Isaia R, Pyle DM (2016) ‘Late-Stage Volatile Saturation as a Potential Trigger for Explosive Volcanic Eruptions’. *Nature Geoscience* 9 (3): 249–54. https://doi.org/10.1038/ngeo2639.

728

729

730

Stock, MJ, Humphreys MCS, Smith VC, Isaia R, Brooker RA, Pyle DM (2018) ‘Tracking Volatile Behaviour in Sub-Volcanic Plumbing Systems Using Apatite and Glass: Insights into Pre-Eruptive Processes at Campi Flegrei, Italy’. *Journal of Petrology* 59 (12): 2463–92. https://doi.org/10.1093/petrology/egy020.

731

732

733

734

Toramaru A, Noguchi S, Oyoshihara S, Tsune A (2008) ‘MND (Microlite Number Density) Water Exsolution Rate Meter’. *Journal of Volcanology and Geothermal Research* 175 (1–2): 156–67. https://doi.org/10.1016/j.jvolgeores.2008.03.035.

735

736

737

Utami, SB, Costa F, Lesage P, Allard P, Humaida H (2021) ‘Fluid Fluxing and Accumulation Drive Decadal and Short-Lived Explosive Basaltic Andesite Eruptions Preceded by Limited Volcanic Unrest’. *Journal of Petrology* 62 (11).

738

739

740

https://doi.org/10.1093/petrology/egab086.

741

Utami SB, Andùjar J, Costa F, Scaillet B, Humaida H, Carn S (2022) ‘Pre-Eruptive Excess Volatiles and Their Relationship to Effusive and Explosive Eruption Styles in Semi-Plugged Volcanoes’. *Frontiers in Earth Science* 10 (September).

742

743

744

https://doi.org/10.3389/feart.2022.882097.
745 Vandemeulebrouck, J, Sabroux J-C, Halbwachs M, Poussielgue N, Grangeon J, Tabbagh J, et al (2000) ‘Hydroacoustic Noise Precursors of the 1990 Eruption of Kelut Volcano, Indonesia’. *Journal of Volcanology and Geothermal Research* 97 (1): 443–56.

748 Wang, L-X, Marks M A W, Wenzel T, Von Der Handt A, Keller J, Teiber H, Markl G (2014) ‘Apatites from the Kaiserstuhl Volcanic Complex, Germany: New Constraints on the Relationship between Carbonatite and Associated Silicate Rocks’. *European Journal of Mineralogy* 26 (3): 397–414. https://doi.org/10.1127/0935-1221/2014/0026-2377.

752 Whelley, PL, Newhall CG, Bradley KE (2015) ‘The Frequency of Explosive Volcanic Eruptions in Southeast Asia’. *Bulletin of Volcanology* 77 (1). https://doi.org/10.1007/s00445-014-0893-8.

755 Widiyantoro S, van der Hilst R D (1996) ‘Structure and Evolution of Lithospheric Slab’. *Science* 271 (March 1996): 1566–70.

757 Widiyantoro S., Ramdhan M, Métaxian J P, Cummins P. R., Martel C, Erdmann S, Nugraha A D, Budi-Santoso A, Laurin A, and Fahmi A A (2018) ‘Seismic Imaging and Petrology Explain Highly Explosive Eruptions of Merapi Volcano, Indonesia’. *Scientific Reports* 8 (1): 1–7. https://doi.org/10.1038/s41598-018-31293-w.

761 Williams, GT, Jenkins, SF, Biass S, Wibowo HE, Harijoko A (2020) ‘Remotely Assessing Tephra Fall Building Damage and Vulnerability: Kelud Volcano, Indonesia’. *Journal of Applied Volcanology* 9 (1): 1–18. https://doi.org/10.1186/s13617-020-00100-5.

764

List of Tables captions

765

766

767 **Table 1:** Apatite diffusion coefficients and timescales modelled using ApTimer (Li et al., 2020), along with best fit magma ascent rate ($v_{ascent}$). Input parameters considered for the 2007 dome apatite profiles span petrologic textures, temperature, acute angles between profiles and main crystallographic axes ($\alpha$ and $\gamma$), and calculated diffusion coefficients for Cl, F, and OH. The best fit timescales were calculated using the Cl or F diffusion along profiles. The Cl diffusion coefficient for the 1990 and 2014 apatite are calculated assuming the fastest diffusion parallel to the c-axis. In rare cases where the F concentration profiles show less fluctuation than the Cl, the F concentration is used instead (see Fig 4e-g for an example).

768

769

770

771

772

773

774

775

776

777

778 **Table 2:** Estimated pre-eruptive conduit geometry of the 1990, 2007, 2014 eruptions. Dense rock equivalent (DRE) volumes, erupted masses, and mean/maximum discharge rates for the eruptions were taken from Maeno *et al.* (2019).

779

780

781

782

783 **Table 3:** Input parameters and results of melt and effective viscosity calculations for 1990, 2007, and 2014 eruptions. Viscosity models used were: Einstein-Roscoe equation (Einstein, 1906; Roscoe, 1952); Caricchi *et al.* (2007); Costa *et al.* (2009); Faroughi and
784 Huber (2015); and Moitra *et al.* (2015, 2018).

785

786 **List of Figures captions**

787

788 **Fig. 1** BSE/SEI images of apatite crystals in various textural environments in the 1990, 2007, and 2014 rocks. Abbreviations are as follows: Ap = apatite; Gm = groundmass; Amph = amphibole (pseudo = pseudomorph); Cpx = clinopyroxene; Opx = orthopyroxene; Mag = magnetite; Ilm = ilmenite; MI = melt inclusion. In 1990 main- and late-stage pumices, apatite is found in: (a) the groundmass; (b) as Amph-hosted inclusion; (c) Cpx-hosted inclusion; (d) Opx-hosted inclusion; (e) inclusions in Amph pseudomorphs; (f) as inclusions in quartz-bearing enclaves; (g) glomerocrysts; (h) embedded between phenocrysts; (i) groundmass phase; (j) Cpx-hosted inclusions; (k) melt inclusion daughter phase; and (l) Amph-hosted inclusions. White scale bars are 10 µm, whereas red scale bars are 100 µm

798

799 **Fig. 2** Ternary plots of volatile compositions (in mol %) of point analysis and profiles of apatite and elemental maps from the 2007 dome. Note that the crystals show significant chemical zoning. Ternary plot of volatile compositions (in mol. %) of 2007 apatite point analyses, with white symbols are the 2007 apatite and black symbols are apatite from other eruptions. (b) Ternary plot showing cores and rims compositions of apatite profiles, where rims are relatively more enriched in Cl compared to cores. Elemental maps show spatial variations in Cl and F concentrations: (c) apatite wedged between two phenocrysts; (d) Cpx-hosted apatite, and (e) groundmass apatite

807

808 **Fig. 3** Ternary plots of volatile compositions (in mol %) of point analysis and profiles of apatite and elemental maps from the 1990 and 2014 eruptions. Note that the crystals show limited to no chemical zoning of these elements. (a) Ternary plot of volatile compositions (in mol. %) of 1990 apatite point analyses, with white symbols are the 1990 apatite and black symbols are apatite from other eruptions. Circles are main-stage apatite, and squares are late-stage apatite. (b) Ternary plot of volatile compositions (in mol. %) of 2014 apatite point analyses. Below this plot is ternary plot showing cores and rims compositions of apatite profiles. X-ray elemental maps show Cl and F distributions. (c) Cpx-hosted inclusions (2014); (d) Amph-hosted apatite (2014), and (e) groundmass

809

810

811

812

813

814

815

816
817 apatite (1990)

818

819 **Fig. 4** Two representative examples of compositional profiles of apatite crystals for the 2007 dome, with diffusion timescales calculated using ApTimer (Matlab version) from Li *et al.* (2020). Blue dashed line of each profile is the initial profile before diffusion. Red solid line shows the diffused model profile. (a) SEI image of groundmass apatite, with red arrow showing the profile. Inset shows a pole figure with the c-axis and a-axes. Orientation of the start of the profile is given as 'trv'. (b) $X_{CI}$ profile from rim to core. (c) $X_{F}$ profile from rim to core. (d) $X_{OH}$ profile from rim to core, which is a mirror image of $X_{F}$. (e) BSE image of apatite embedded in phenocrysts, with red arrow showing the profile. Inset shows pole figure with similar information as in (a). (f) $X_{CI}$ profile from rim to core. (g) $X_{F}$ profile from rim to core. (h) $X_{OH}$ profile from rim to core, which increases as in $X_{F}$

830

831 **Fig. 5** Histogram of timescales of zoned crystals from the 2007 dome of Kelud. The timescales are divided into two groups: timescales where the acute angles of the profile is closer to the c-axis (light green), and timescales where the acute angles of the profile is closer to the a-axis (dark green)

835

836 **Fig. 6** Viscosity vs. ascent rate and mass discharge rate, modified after Fig. 4.10 of Gonnerman and Manga (2013). The magma ascent rates were derived from apatite diffusion modelling and viscosity calculated for Kelud magmas were estimated using the viscosity model of Moitra *et al.* (2015, 2018). Coloured boxes represent the region outlined by magma viscosity vs. ascent velocity observed for the 1990 (grey), 2007 (red), and 2014 (blue) eruptions for a range of dissolved water content in the melt. Boxes with dashed lines are the average to maximum mass discharge rates determined by Maeno *et al.* (2019)

844

845 **List of Supplementary Materials**

846

847 **Fig. S1** Time-series of radial and tangential tilt measurements before the 2007 eruption, compared to the diffusion timescales from apatite, modified after Fig. 3b of Hidayati et al. (2019)

849
850

851

**Fig. S2** Time-series of seismicity (RSAM, tremor, low frequency (LF), shallow volcano-tectonic (VB), and volcano-tectonic (VT)) before the 2007 eruption, compared to the diffusion timescales from apatite. The figure is modified after Fig. 7a of Hidayati et al. (2019)

852

853

854

855

856

**Fig. S3** Time-series of cumulative energy and earthquake depths before the 2007 eruption compared to diffusion timescales obtained from chemical zoning in apatite. Monitoring data and figure are modified after Fig. 9 of Hidayati *et al.* (2019)

857

858

859

860

**Table S1** Compilation of unrest data for the 2007 and 2014 eruptions of Kelud, compiled from Cassidy *et al.* (2019), Caudron *et al.* (2012), Hidayati *et al.* (2009, 2019), and Lubis (2014)

861

862

863

864

**Table S2** Apatite compositions obtained from point analyses of individual grains from the 1990, 2007, and 2014 eruptions

865
Figures

<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

BSE/SEI images of apatite crystals in various textural environments in the 1990, 2007, and 2014 rocks.
Abbreviations are as follows: Ap = apatite; Gm = groundmass; Amph = amphibole (pseudo = pseudomorph); Cpx = clinopyroxene; Opx = orthopyroxene; Mag = magnetite; Ilm = ilmenite; MI = melt
inclusion. In 1990 main- and late-stage pumices, apatite is found in: (a) the groundmass; (b) as Amph-
hosted inclusion; (c) Cpx-hosted inclusion; (d) Opx-hosted inclusion; (e) inclusions in Amph
pseudomorphs; (f) as inclusions in quartz-bearing enclaves; (g) glomerocrysts; (h) embedded between
phenocrysts; (i) groundmass phase; (j) Cpx-hosted inclusions; (k) melt inclusion daughter phase; and (l)
Amph-hosted inclusions. White scale bars are 10 mm, whereas red scale bars are 100 mm

<figure><img src="image_3.png" /><figcaption>**Figure 2**</figcaption></figure>
Ternary plots of volatile compositions (in mol %) of point analysis and profiles of apatite and elemental maps from the 2007 dome. Note that the crystals show significant chemical zoning. Ternary plot of volatile compositions (in mol. %) of 2007 apatite point analyses, with white symbols are the 2007 apatite and black symbols are apatite from other eruptions. (b) Ternary plot showing cores and rims compositions of apatite profiles, where rims are relatively more enriched in Cl compared to cores. Elemental maps show spatial variations in Cl and F concentrations: (c) apatite wedged between two phenocrysts; (d) Cpx-hosted apatite, and (e) groundmass apatite

<figure><img src="image_4.png" /><figcaption>(b) 2014 apatite analyses
(points + traverse)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) 2014 clinopyroxene-hosted apatite, 14PUMI2-2-px3-apa1</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) 2014 amphibole-hosted apatite, 14PUMI1-amph-apa6</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) 1990 main-stage groundmass apatite, PL1new-px4-apa1</figcaption></figure>

<figure><img src="image_8.png" /></figure>
Figure 3

Ternary plots of volatile compositions (in mol %) of point analysis and profiles of apatite and elemental maps from the 1990 and 2014 eruptions. Note that the crystals show limited to no chemical zoning of these elements. (a) Ternary plot of volatile compositions (in mol. %) of 1990 apatite point analyses, with white symbols are the 1990 apatite and black symbols are apatite from other eruptions. Circles are main-stage apatite, and squares are late-stage apatite. (b) Ternary plot of volatile compositions (in mol. %) of 2014 apatite point analyses. Below this plot is ternary plot showing cores and rims compositions of apatite profiles. X-ray elemental maps show Cl and F distributions. (c) Cpx-hosted inclusions (2014); (d) Amph-hosted apatite (2014), and (e) groundmass apatite (1990)
<figure><img src="image_9.png" /><figcaption>Figure 4</figcaption></figure>

Two representative examples of compositional profiles of apatite crystals for the 2007 dome, with diffusion timescales calculated using ApTimer (Matlab version) from Li *et al.* (2020). Blue dashed line of each profile is the initial profile before diffusion. Red solid line shows the diffused model profile. (a) SEI image of groundmass apatite, with red arrow showing the profile. Inset shows a pole figure with the c-axis and a-axes. Orientation of the start of the profile is given as 'trv'. (b) $X_{\text{CI}}$ profile from rim to core. (c) $X_{\text{F}}$
profile from rim to core. (d) X<sub>OH</sub> profile from rim to core, which is a mirror image of X<sub>F</sub>. (e) BSE image of apatite embedded in phenocrysts, with red arrow showing the profile. Inset shows pole figure with similar information as in (a). (f) X<sub>CI</sub> profile from rim to core. (g) X<sub>F</sub> profile from rim to core. (h) X<sub>OH</sub> profile from rim to core, which increases as in X<sub>F</sub>

<figure><img src="image_10.png" /><figcaption>Figure 5</figcaption></figure>

Histogram of timescales of zoned crystals from the 2007 dome of Kelud. The timescales are divided into two groups: timescales where the acute angles of the profile is closer to the c-axis (light green), and timescales where the acute angles of the profile is closer to the a-axis (dark green)
<figure><img src="image_11.png" /><figcaption>Figure 6</figcaption></figure>

Viscosity vs. ascent rate and mass discharge rate, modified after Fig. 4.10 of Gonnerman and Manga (2013). The magma ascent rates were derived from apatite diffusion modelling and viscosity calculated for Kelud magmas were estimated using the viscosity model of Moitra *et al.* (2015, 2018). Coloured boxes represent the region outlined by magma viscosity vs. ascent velocity observed for the 1990 (grey), 2007 (red), and 2014 (blue) eruptions for a range of dissolved water content in the melt. Boxes with dashed lines are the average to maximum mass discharge rates determined by Maeno *et al.* (2019)
Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

* FigS1traverseapatiteunrest2007tilt.pdf

• FigS2traverseapatiteunrest2007seismicity.pdf

• FigS3traverseapatiteunrest2007hyp.pdf

• MS3TableS1preeruptiveunrest0714.xlsx

• MS3TableS2apatitepointAnalyse.xlsx

• Table1ApTimertimescalesv4.xlsx

• Table2MagmaDischargeRateConduitv2.xlsx

* Table3viscositycalculationsv2.xlsx