# A model for permeability evolution during volcanic welding

Fabian Wadsworth, Jérémie Vasseur, Edward Llewellin, Richard Brown, Hugh Tuffen, James Gardner, Jackie Kendrick, Yan Lavallée, Katherine Dobson, Michael Heap, et al.

► To cite this version:

Fabian Wadsworth, Jérémie Vasseur, Edward Llewellin, Richard Brown, Hugh Tuffen, et al. A model for permeability evolution during volcanic welding. Journal of Volcanology and Geothermal Research, 2021, 409, pp.107118. 10.1016/j.jvolgeores.2020.107118. hal-03547493

HAL Id: hal-03547493

https://hal.science/hal-03547493

Submitted on 29 Jan 2022

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
<figure><img src="image_2.png" /></figure>

# A model for permeability evolution during volcanic welding

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

Check for
updates

Fabian B. Wadsworth a,*, Jérémie Vasseur b, Edward W. Llewellin a, Richard J. Brown a,
Hugh Tuffen c, James E. Gardner d, Jackie E. Kendrick e, Yan Lavallée f, Katherine J. Dobson g,
Michael J. Heap h, Donald B. Dingwell b, Kai-Uwe Hess b, Jenny Schauroth f, Felix W. von Aulock f,
Alexandra R.L. Kushnir i, Federica Marone j

## ARTICLE INFO

**Article history:**

Received 15 September 2020

Received in revised form 5 November 2020

Accepted 6 November 2020

Available online 10 November 2020

**Keywords:**

Tuffisite

Compaction

Sintering

Ignimbrite

Synchrotron

Experiment

Magma

## ABSTRACT

Volcanic ash and pyroclasts can weld when deposited hot by pyroclastic density currents, in near-vent fall deposits, or in fractures in volcano interiors. Welding progressively decreases the permeability of the particle packs, influencing a range of magmatic and volcanic processes, including magma outgassing, which is an important control on eruption dynamics. Consequently, there is a need for a quantitative model for permeability evolution during welding of ash and pyroclasts under the range of conditions encountered in nature. Here we present in situ experiments in which hydrous, crystal-free, glassy pyroclasts are imaged via x-ray tomography during welding at high temperature. For each 3D dataset acquired, we determine the porosity, Darcian gas permeability, specific surface area, and pore connectivity. We find that all of these quantities decrease as a critical percolation threshold is approached. We develop a constitutive mathematical model for the evolution of permeability in welding volcanic systems based on percolation theory, and validate the model against our experimental data. Importantly, our model accounts for polydispersivity of the grainsize in the particle pack, the pressures acting on the pack, and changes in particle viscosity arising from degassing of dissolved H₂O during welding. Our model is theoretically grounded and has no fitting parameters, hence it should be valid across all magma compositions. The model can be used to predict whether a cooling pyroclast pack will have sufficient time to weld and to degas, the scenarios under which a final deposit will retain a permeable network, the timescales over which sealing occurs, and whether a welded deposit will have disequilibrium or equilibrium H₂O content. A user-friendly implementation of the model is provided.

© 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Magma may fragment during ascent, forming volcanic ash and pyroclasts, which are subsequently transported in a gas phase. Following deposition the particles may weld – i.e. coalesce to form a coherent material – if they remain hot enough to deform viscously. Welded deposits are observed in a number of volcanic settings (Fig. 1). Subaerially, welded spatter deposits (Sumner et al., 2005; Sánchez et al., 2012) and welded fall deposits (Houghton and Carey, 2015) may form proximal to

the vent, and welding may occur in the depositional boundary layer at the base of a pyroclastic density current (Branney and Kokelaar, 2002) or after deposition and cessation of flow (Riehle, 1973; Heap et al., 2014; Lavallée et al., 2015). In the volcanic plumbing system, welding may occur in vent-filling ignimbrites (Almond, 1971), at conduit margins (Gardner et al., 2017; Wadsworth et al., 2020a), and in tuffisite veins (Tuffen et al., 2003; Tuffen and Dingwell, 2005; Stasiuk et al., 1996; Kendrick et al., 2016; Kolzenburg et al., 2019) which are produced when gas-particle mixtures hydraulically fracture and intrude the

<footnote>a Department of Earth Sciences, Durham University, Durham DH1 3LE, United Kingdom</footnote>

<footnote>b Department of Earth and Environmental Sciences, Ludwig-Maximilians-Universität, Munich 80333, Germany</footnote>

<footnote>c Lancaster Environment Centre, Lancaster University, Lancaster LA1 4YQ, United Kingdom</footnote>

<footnote>d Department of Geological Sciences, Jackson School of Geosciences, 2305 Speedway Stop C1160, Austin, TX 78712-1692, United States</footnote>

<footnote>e School of Geosciences, University of Edinburgh, The King's Buildings, James Hutton Road, Edinburgh EH9 3FE, Scotland, United Kingdom</footnote>

<footnote>f School of Environmental Sciences, University of Liverpool, Liverpool L69 3BX, United Kingdom</footnote>

<footnote>g Civil and Environmental Engineering, University of Strathclyde, Glasgow G1 1XJ, United Kingdom</footnote>

<footnote>h Université de Strasbourg, CNRS, Institut de Physique du Globe de Strasbourg UMR 7516, F-67000 Strasbourg, France</footnote>

<footnote>i Université de Strasbourg, CNRS, École et Observatoire des Sciences de la Terre UMS 830, F-67000 Strasbourg, France</footnote>

<footnote>j Swiss Light Source, Paul Scherrer Institute, Villigen PSI 5232, Switzerland</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: fabian.b.wadsworth@durham.ac.uk (F.B. Wadsworth).</footnote>
country rock or overlying/surrounding magma (Heiken et al., 1988; Stasiuk et al., 1996; Tuffen et al., 2003; Owen et al., 2019).

Welding is important because it leads to substantial changes in the physical properties of the deposit, which transitions progressively from particulate to coherent (Branney and Kokelaar, 1992), from mechanically weak to mechanically strong (Kolzenburg et al., 2012; Vasseur et al., 2013), and from permeable to impermeable (Colombier et al., 2017; Farquharson et al., 2017; Wadsworth et al., 2017b; Kolzenburg et al., 2019; Heap et al., 2019). The effect on permeability is of particular interest because it has been proposed that closure of outgassing pathways through tuffisites (Farquharson et al., 2017; Heap et al., 2019) and clastic vent-filling deposits (Quane et al., 2009) can cause vent plugging and build-up of gas pressure, leading to explosive eruption.

The physical process of welding of silicate melts is synonymous with viscous sintering as it is understood in materials science research (e.g. Mackenzie and Shuttleworth, 1949; Prado et al., 2001; Wadsworth et al., 2016a); however, welding in natural environments presents added complications not usually encountered in materials research. These include: (1) the particles may be supersaturated in volatiles and therefore may degas during welding (Sparks et al., 1999; Wadsworth et al., 2019; Heap et al., 2019); (2) the particles may be undersaturated in volatiles and undergo resorption during welding (Sparks et al., 1999; Gardner et al., 2018); (3) the particles may be variably crystalline, or internally porous (Quane and Russell, 2005a; Wright and Cashman, 2014; Heap et al., 2015; Kendrick et al., 2016); and (4) the sintering pack of particles may be subjected to local stresses that promote shearing (e.g. rheomorphism; Andrews and Branney, 2011), or to substantial confining pressure (Sparks et al., 1999). Furthermore, welding is strongly affected by the thermal pathway that particles follow after fragmentation (Branney and Kokelaar, 2002; Wadsworth et al., 2014, 2019; Kolzenburg et al., 2019; Heap et al., 2019), which may be complex.

We develop a general model for evolution of permeability with time during welding of deposits of crystal-poor volcanic particles, which we validate against experimental data. The important novelty over previous work on welding dynamics (Riehle, 1973; Riehle et al., 1995; Sparks et al., 1999; Quane and Russell, 2005a; Russell and Quane, 2005; Lavallée et al., 2015) is that the model accounts for a range of factors that complicate welding in volcanic systems: particle polydispersivity, non-isothermal conditions, degassing of dissolved H₂O during welding, gas pressure, and confining pressure. The model is therefore suitable for application to welding in a wide range of volcanic scenarios. We apply our model to welding of rhyolites in a range of regimes, with the principal aim of exploring the conditions under which a welded deposit is expected to preserve equilibrium degassed volatile concentrations, the conditions under which disequilibrium volatile concentrations would be preserved, and the conditions under which welding is arrested before permeability decreases to zero.

An implementation of the model is available via Vhub (https://vhub.org/resources/4568) and a simplified version implemented in Excel™ is provided as Supplementary Material A. A practical guide to its use is presented in Section 7.

## 2. Theoretical development

The welding of hot pyroclasts (above their glass transition temperature) is a specific case of viscous sintering. In physical volcanology, 'sintered' has commonly been used as a descriptor for poorly or incipiently welded deposits (Quane and Russell, 2005b; Wright and Cashman, 2014), but in materials science, viscous sintering and volcanic welding are the same fundamental process, describing the progressive coalescence of droplets toward a low-porosity, dense end-state. In this work we therefore use the terms sintering and welding interchangeably. Welding is driven by stresses arising from surface tension and confining stress (if applicable) and is opposed by the viscosity of the droplets and by pressure in the interstitial gas phase. As welding progresses, the geometry of the pore space and the degree of pore

connectivity evolve, therefore causing the permeability to evolve. In this section, we develop a mathematical model for the evolution of the permeability of welding, volatile-supersaturated magmatic particles, building on previous models for the evolution of porosity of welding ash, and for the relationship between permeability and porosity for particulate and percolating materials.

### 2.1. Evolution of porosity during welding

We adopt the model of Wadsworth et al. (2016a, 2017b, 2019), who treat the welding system as a collection of unit cells that abstract the pore-particle system as a spherical 'bubble' of gas surrounded by a shell of viscous melt; the gas is allowed to escape freely from the bubble, giving the model its name: the 'vented bubble model' (Wadsworth et al., 2016a, 2016b). An overview of the main equations from Wadsworth et al. (2019) is presented in this section; the reader is referred to the original work for their derivation, solution, and experimental validation. The model gives the rate of change of the porosity (i.e. gas volume fraction) $\phi$ with time $t$ as a function of material properties:

$$ \frac{d\phi}{dt} = -\frac{3\Delta P}{4\mu}\phi - \frac{3\Gamma}{2\mu a_i}\left(\frac{\phi_i}{1-\phi_i}\right)^{1/3}\phi^{2/3}(1-\phi)^{1/3} \quad \text{Eq. 1} $$

where $\phi_i$ is the initial porosity, $a_i$ is the initial radius of the 'bubbles', $\mu$ is the melt viscosity, $\Gamma$ is the surface tension, and $\Delta P$ is the difference between the squeezing isotropic pressure $P$ acting on the particles and the gas pressure in the bubbles $P_g$, such that $\Delta P = P - P_g$. In nature, $\Delta P$ is controlled by a balance between the overburden pressure and the pore fluid pressure.

Eq. 1 can be presented in dimensionless form to facilitate collapse of experimental data over a range of conditions, and to identify regimes of behavior (Wadsworth et al., 2019):

$$ \frac{d\bar{\phi}}{dt} = -\frac{3}{2} \left[ \bar{P}\phi + \left( \frac{1-\phi_i\bar{\phi}}{1-\phi_i} \right)^{1/3} \phi^{2/3} \right] \quad \text{Eq. 2} $$

where, $\bar{t}$, $\bar{P}$, and $\bar{\phi}$ are, respectively, dimensionless time, pressure difference, and porosity. Time is normalized by the capillary relaxation time-scale $\lambda = \mu a_i/\Gamma$, pressure difference is normalized by a reference capillary pressure scale $P_L = 2\Gamma/a_i$, and porosity is normalized by the initial porosity $\phi_i$ as follows:

$$ \bar{t} = \frac{t}{\lambda} = \frac{\Gamma}{a_i} \int_{t_i}^{t} \frac{1}{\mu} dt; \quad \bar{P} = \frac{\Delta P}{P_L} = \frac{(P-P_g)a_i}{2\Gamma}; \quad \bar{\phi} = \frac{\phi}{\phi_i} \quad \text{Eq. 3} $$

Melt viscosity is a strong function of temperature $T$, which may vary over time, and dissolved H₂O concentration $C$, which may vary with both time and spatial position within the particles if they are degassing during welding. This variation is accounted for by determining the spatial distribution of dissolved H₂O within the particles at each time via Fick's second law of diffusion, cast in spherical coordinates (Crank, 1975), and computing the spatial average $\langle C \rangle$. The spatially-averaged viscosity $\mu = f(T, \langle C \rangle)$ is then used in the integral within the definition of $\bar{t}$, which runs from the time at which the welding starts $t_i$ to the time-step of interest (Eq. 3).

The equations for the rate of porosity change – in dimensional form (Eq. 1) or dimensionless form (Eq. 2) – can be integrated numerically to derive porosity change as a function of time $\phi(t)$ for welding particles. Wadsworth et al. (2019) describe a suitable numerical scheme and provide a user-friendly implementation for download. The full model accounts for particle polydispersivity, non-isothermal conditions, degassing during welding, and confining pressure. Polydispersivity of particles in the welding deposit manifests as a distribution of initial particle radii and interstitial pore radii. This is accounted for via convolution techniques discussed in Wadsworth et al. (2017c) and used here.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Settings in which welding of volcanic particles occurs. (a) A schematic of a generic explosive eruption (re-drawn and adapted from Harris et al., 2011). Welding may occur in tuffisites (black veins) which hydraulically fracture country rock or overlying magma, in the deposits of pyroclastic density currents on the volcano flanks, and in aircraft jet engines. Here, we focus on tuffisites and ignimbrites. (b) A tuffisite hosted in obsidian from the 2008 eruption of Volcán Chaiten (image from Heap et al., 2019). (c) Detail of a partially welded ignimbrite from the 668 ka phonolitic Arico eruption of Las Cañadas volcano, Tenerife (image courtesy of Rich Brown).</figcaption></figure>

## 2.2. Evolution of permeability during welding

Section 2.1 presents a model for the evolution of porosity with time $\phi(t)$ in a welding deposit. We can extend this to derive a model for the evolution of Darcian permeability with time $k(t)$ by finding a constitutive relationship for $k(\phi)$. For sintering systems of particles, the most widely validated model is based on percolation theory (Martys et al., 1994; Wadsworth et al., 2016b):

$$k = \frac{2[1-(\phi-\phi_c)]}{s^2} (\phi-\phi_c)^b \qquad \text{Eq.4}$$

where $s$ is the specific surface area, $\phi_c$ is the percolation threshold, and $b$ is a percolation exponent. For volcanic applications, this model has been validated for sintering glass spheres (Wadsworth et al., 2017b, 2020b; Eichheimer et al., 2020), volcanic tuffisites (Heap et al., 2019), welded impact breccias (Heap et al., 2020), and for sintered volcanic fault gouge (Ryan et al., 2020a, 2020b). Whereas $k$ and $\phi$ are typically measured quantities, $s$ is either measured or requires estimation, and $\phi_c$ and $b$ are usually treated as fit parameters (e.g. Mueller et al., 2005; Colombier et al., 2017).

Wadsworth et al. (2016a) proposed that, from a micro-structural perspective, sintering systems of particles are a dynamic analogue of static overlapping sphere geometries, which are frequently used to consider the properties of low-porosity porous media (e.g. Martys et al., 1994). By drawing this analogy, we can use continuum percolation theory and simulation results to constrain $b$ and $\phi_c$, removing the need to fit for an unknown parameter. Based on theoretical arguments, the percolation exponent for overlapping sphere systems can be constrained to be $b = 4.4$ (Feng et al., 1987), which has been found to match both experimental and simulation data for overlapping spheres and systems of sintering particles (Vasseur and Wadsworth, 2017). For monodisperse overlapping spheres, simulations show that $\phi_c \approx 0.03$ (Kertesz, 1981; Elam et al., 1984; Vasseur and Wadsworth, 2017). These constraints are for idealized sphere systems and require validation for sintering of non-spherical particles.

In our experimental campaign presented here (see Section 3), we can measure $s$ directly. However, to allow Eq. 4 to be applied in situations where $s$ is not measured, we also test a constitutive model for $s(\phi)$. By again drawing an analogy between dynamic sintering systems

of particles and static simulations of overlapping sphere domains, we can use a theoretical constraint for $s$ (Torquato, 2013)

$$s = -3(1-\phi) \ln(1-\phi) \frac{\langle a^2 \rangle}{\langle a^3 \rangle} \qquad \text{Eq.5}$$

where $\langle a^n \rangle$ is the nth moment of the distribution of pore sizes. Wadsworth et al. (2017c) provide a method for predicting the distribution of pore sizes, and therefore $\langle a^n \rangle$, in sintering systems that have a distribution of initial particle sizes.

The above framework provides a method by which $\phi(t)$ and $k(\phi)$ can be predicted. In regimes where the permeability does not limit the gas escape and closure of the pore network by sintering (discussed in Wadsworth et al., 2019), the time-dependence of permeability can then be computed a priori by combining these two models, resulting in model predictions that can be tested against experiment.

<figure><img src="image_5.png" /><figcaption>Fig. 2. The particle size distribution used in the experiments presented in Wadsworth et al. (2019) and re-analyzed herein. The values for the first, second, and third moments of the distribution are given on the figure and the first moment (mean) is indicated with a vertical dashed line. The data are binned at 0.04 log units in R.</figcaption></figure>
# 3. Methods

## 3.1. *Experiments and data analysis*

We analyze the results from high temperature welding experiments using volatile supersaturated rhyolitic obsidian particles, presented in Wadsworth et al. (2019). The rhyolite particles had a mean size $\langle R \rangle = 2.33 \times 10^{-5}$ m, were polydisperse in grain size, and were crystal-free. The size distribution of particles is reported in Fig. 2 along with $\langle R^n \rangle$ for $n \le 3$.

The rhyolite particles were packed into free-standing cylinders 3 mm high using a compression forming method (Wadsworth et al., 2016b) and were heated using a laser heating system (Fife et al., 2012). Temperature was monitored using a calibrated pyrometer and additional calibration steps (see Wadsworth et al., 2017b). Temperature fluctuations were greater than the absolute uncertainty on temperature and were maximally ~5 K (i.e. ± 2.5 K). Four runs were performed, in which the packed free-standing cylinders were heated at 0.067 K. $s^{-1}$ to temperatures $T = 1053, 1144, 1217$ and $1248$ K. The experiments were performed in situ at the TOMCAT beamline of the Swiss Light Source and imaged simultaneously via x-ray computed tomography,

resulting in a time-series of 3D reconstructions of the evolving geometry of the particle pack throughout the welding process. Full 3D scans were collected at intervals between 1 and 107 s with a spatial resolution of 1.7 – 2.0 µm voxel edge lengths. The initial H₂O concentration was 0.14 wt% (Wadsworth et al., 2019).

A sub-set of the 3D datasets in each run were segmented to reconstruct the non-pore volume (i.e. the volume of the welding particle phase) and the pore volume between the particles, using Avizo™ using an auto-thresholding algorithm. A time-series of representative images of segmented volumes are shown in Fig. 3. The porosity $\phi$ is calculated directly from the tomography data as the ratio of the measured pore volume to the total volume (non-pore volume plus pore volume). The pore volume is divided into connected porosity $\phi_p$ (i.e. pores and pore networks that connect across at least one opposite face of the measured domain) and isolated porosity, such that the isolated porosity is $\phi - \phi_p$. Using a range of suitable segmentations, we determine the uncertainty on the porosity to be approximately 0.5% of the stated values. This is a conservative estimate based on variations in the segmentation thresholds used. During welding, the pore volume decreases and the sample shrinks (Fig. 3a-c). The pore space is initially fully connected (Fig. 3a) but develops domains that are isolated (Fig. 3b) as welding

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Internal textures in packs of volcanic particles that are welding at high temperature, imaged using in situ x-ray tomography. (a-c) 3D renderings of the pore spaces using grey to denote pores that are connected from edge-to-edge and green to denote pores that are isolated from through-going connections (Insets: the full exterior sample shape in 2D showing the evolution of relative 2D area of the side of the free-standing cylinders. (d-f) 2D slices of the same sample time-steps given in (a-c). (g-i) Collapsed 2D views of the steady-state fluid flow vector distribution computed using LBflow and visualized using Paraview. Red and blue colours respectively indicate relatively high and low flow speeds; the absolute values of the colour scale are arbitrary and for comparison only. Note that the scale of the cube edge in (a-c) and the square edge in (d-i) is 250 microns, and that the square 2D slices given in (d-f) are taken at the mid-point of the z-axis from the boxes in (a-c).</figcaption></figure>
progresses. The stable end-state is a melt continuum containing isolated bubbles (Fig. 3c); the topological inverse of the initial state (Fig. 3a).

In this study we use these tomography data – results from which are presented in Wadsworth et al. (2019) – to determine the permeability of the samples through the welding process. The segmented pore volume is used as the input for LBflow, a numerical package for simulating fluid flow using the lattice-Boltzmann method (Llewellin, 2010a, 2010b). LBflow discretizes the pore volume into a cubic lattice of fluid nodes, to which we add a periodic boundary and a buffer zone to allow flow to connect across faces, and simulates fluid flow under an applied body force, equivalent to applying a pressure gradient across the sample. We simulate the flow of air at ambient conditions (fluid viscosity $\mu_f = 1.8205 \times 10^{-5}$ Pa. s and density $\rho_f = 1.2047$ kg. m$^{-3}$) under a uniform fluid pressure gradient of $\nabla P_f = 10^{-2}$ Pa. m$^{-1}$. The fluid velocity $u$ is read out at every node and the average fluid velocity $\langle u \rangle$ is computed. We run the simulation until $\langle u \rangle$ stabilizes to a steady state value, and ensure the flow is always in the low-Reynolds number (i.e. creeping flow) regime (for more details, see Llewellin, 2010b). Visualizations of computed fluid flow lines through representative segmented volumes are shown in Fig. 3. The permeability is calculated using Darcy’s law for each 3D dataset:

$$k = - \frac{\mu_f}{\nabla P_f} \langle u \rangle \qquad \text{Eq.6}$$

The specific surface area $s$ for each the 3D dataset is determined using a ‘marching cubes’ algorithm (Lorensen and Cline, 1987; Lewiner et al., 2003). This algorithm extracts an isosurface for the pore-particle walls in the 3D datasets. Normalizing the area of the surface by the volume of the sample or domain gives $s$ (in units of m$^{-1}$). Using 3D datasets segmented using different segmentation threshold choices, we estimate the uncertainty on the permeability to be $10^{-17}$ m$^2$ based on the convergence criterion used in LBflow. The uncertainty on our surface area determination is approximately 0.8%. This uncertainty is determined by computing the standard error on repeat measurements of surface area on different sub-volumes selected in domain.

### 3.2. Numerical modelling

We solve Eq. 2 numerically, following Wadsworth et al. (2019), to obtain $\bar{\phi}(t)$, which we dimensionalize for each experiment via the relationships in Eq. 3, and the material parameters $\mu$, $\Gamma$, $\langle a_i \rangle$, and $\phi_i$. Calculating $\mu$ requires the spatial distribution of dissolved H$_2$O in the particles at every timestep throughout the welding process, which we compute via Fick’s second law, following Wadsworth et al. (2019). We then integrate across the measured particle radii to find a spatially averaged viscosity and convolve this with the initial particle size distribution (Fig. 2) to arrive at a global average of the viscosity for the population of particles used – termed $\langle u \rangle$. We use constitutive models appropriate for rhyolitic liquids: viscosity $\mu(T, C)$ from Hess and Dingwell (1996); H$_2$O diffusivity $D(T, P_g, C)$ from Zhang and Ni (2010); and the equilibrium H$_2$O solubility at the particle boundary $C_c(T, P_{H2O})$ from Liu et al. (2005), where $P_{H2O}$ is the partial pressure of H$_2$O in the pore space. For our experiments $P = 0.1$ MPa, and we use $P_{H2O} = 0.1 P_g$ to account for standard laboratory humidity (Von Aulock et al., 2017). We use $\Gamma = 0.3$ N.m$^{-1}$ for nominally dry silicate melts (Parikh, 1958). The average initial pore radius $\langle a_i \rangle$ for the experiments is determined using Wadsworth et al. (2017c), and the starting porosity $\phi_i$ is taken directly from segmented tomography data.

## 4. Results and analysis

Our results and analysis comprise three parts. First, we validate the $s$ ($\phi$) model (Eq. 5; Section 2.2) against the values determined from the tomographic data. Second, we compare the $k(\phi)$ model (Eq. 4; Section 2.2) with $k(\phi)$ data from the LBflow simulations of the in situ

<figure><img src="image_7.png" /><figcaption>Fig. 4. The relationship between specific surface area $s$ and porosity $\phi$ during high-temperature welding of a pack of volcanic particles. Both $s$ and $\phi$ are determined from analysis of tomographic data (Section 3.1).</figcaption></figure>

data, for which all quantities are measured directly. Third, we validate the $k(t)$ model (i.e. the combined $k(\phi)$ and $\phi(t)$ models (Eqs. 1–3; Section 2.1)) against $k(t)$ data from the LBflow simulations.

### 4.1. The specific surface area during welding $s(\phi)$

The specific surface area $s$ measured during welding in the evolving in situ datasets varies non-linearly with $\phi$. In Fig. 4 we show the measured $s$ values compared with the model (Eq. 5; Section 2.2). When solving the model we use the moments of a polydisperse pore size distribution calculated from a polydisperse overlapping sphere model used for sintering problems (Wadsworth et al., 2017c), which uses the particle size distribution (Fig. 2) and the initial porosity $\phi_i$ as input values. In the case of the data presented here, the pore size distribution moments are $\langle a \rangle = 8.01 \times 10^{-6}$ m, $\langle a^2 \rangle = 9.40 \times 10^{-11}$ m$^2$, and $\langle a^3 \rangle = 1.35 \times 10^{-15}$ m$^3$. With no fitting parameters, the data and the model agree with a coefficient of determination of $r^2 = 0.8$.

We hypothesize that the deviation between the model and the data at high porosity is attributable to the initial angularity of the particles, consistent with the observation that initially-spherical particles do not show this discrepancy during sintering (Wadsworth et al., 2017b), and that the data approach the model more closely as the particles progressively round during welding. The magnitude of the deviation between the measured and the predicted $s$ taken at the initial packing porosity – termed $\Delta s_i$ – can be used to estimate the excess surface. We find that $\Delta s_i \langle a_i \rangle \approx 0.2$, which implies that the surface is approximately 20% greater than it would be for a pack of spherical particles with the same size distribution. A sphere has a single radius of curvature and, for a given volume, a sphere has the lowest surface area of any object. Any solid object of the same volume that has a higher surface area must necessarily have non-uniform surface curvature which includes regions with smaller radius of curvature than a sphere. In order to find a conservative estimate for the minimum radius of curvature for such an object, we consider a prolate spheroid with the same volume as a sphere of radius $a$, but a surface area that is 20% higher. Such a spheroid has major semi-axis $c \approx 2.16a$ and minor semi-axes $m \approx 0.68a$. The minimum radius of curvature of the spheroid is $m^2/c \approx 0.21a$ (at each pointed end); hence an object with a surface area that is 20% larger than a sphere of the same volume must have at least one radius of curvature that is no larger than 21% of the sphere’s radius.

We can use this information to estimate the maximum time taken for the irregular particles in our experiments to round under the action of surface tension. The timescale for relaxation of a surface irregularity
with radius of curvature $\xi$ is given by $\lambda_{\xi} = \xi\mu/\Gamma$ (Wadsworth et al., 2017a). For $\xi = 0.21(a_i) = 1.68 \times 10^{-6}$ m, the rounding time for our experiments at 1248 K is $\lambda_{\xi} = 126$ seconds. In our data, the porosity at which $t = \lambda_{\xi}$ is met is approximately at a porosity of $\phi \approx 0.4$, which is consistent with the porosity at which $\Delta s_i$ approaches zero and the model and data coincide (Fig. 4). This implies that particle irregularity may only be important during the very initial stage of welding, and becomes irrelevant quickly. We propose that this simple test supports the hypothesis that initial surface angularity of the particles is responsible for the deviation from the model for specific surface at high porosities. We conclude that, for the case of initially angular particles, our model becomes more valid as welding progresses. Empirical adjustments could be made to the model for $s(\phi)$ to account for angularity (e.g. by adjusting the moments of the pore size distribution by an empirical factor).

## 4.2. The relationship between permeability and porosity $k(\phi)$

The permeability of the samples determined through fluid-flow simulations depends strongly on porosity (Fig. 5a – note the logarithmic scale for $k$). Permeability varies by two orders of magnitude over the porosity interval $0.1 \le \phi \le 0.6$ and, at the lowest porosities, simulations show that the samples are impermeable (i.e. $k = 0$). Fig. 5a also plots the theoretical $k(\phi)$ relationship in Eq. 4 using the model for $s(\phi)$ from Eq. 5 (see Section 4.1), $b = 4.4$, and $\phi_c = 0.03$. This shows a generally good match to the data.

Improvements in the fit of $k(\phi)$ (Eqs. 4 & 5) could be made if $\phi_c$ and $b$ were treated as fitting parameters. In the inset to Fig. 5a, we perform a minimization exercise in which we allow $\phi_c$ and $b$ to vary. Our minimization parameter is the sum of the square residuals between the data and the model result, where we define the residuals on the basis of the Deming approach in which both the mismatch in $k$ and $\phi$ are accounted for. The coefficient of determination is then the ratio of this sum of square residuals for a given $\phi_c$ and $b$ combination, and the total variance in the data. In the inset to Fig. 5a, we show the results of this fitting exercise, demonstrating that the best-fit combinations of $\phi_c$ and $b$ co-vary, such that the residual function has a trench-like minimum. In Fig. 5a, we show a grey band, which bounds all solutions to Eq. 4 for which the coefficient of determination is within 2% of the absolute best fit pair of $\phi_c$ and $b$. We note that this 2% variation in the fit quality encompasses the theoretically-grounded values $\phi_c = 0.03$ and $b = 4.4$, derived for monodisperse packs of spherical particles. Based on this analysis, we conclude that: 1) the quality of the model fit to the $k(\phi)$ data is not sensitive to small variations ($\pm 2\%$) in the best-fit parameters; and 2) that the theoretically-grounded values of the fit

parameters remain valid for packs of natural volcanic particles, despite violation of the monodisperse and spherical assumptions. We therefore adopt $\phi_c = 0.03$ and $b = 4.4$ as empirically-validated, theoretically-grounded constants.

Due to the in situ nature of our experiments, we extract the connected porosity $\phi_p$ directly. If we then look at the connectivity of the pore space $\phi_p/\phi$ (as defined in Colombier et al., 2017) as a function of $\phi$ during welding, we would predict that this metric would fall to zero as $\phi \to \phi_c$ during welding. The value $\phi_c = 0.03$ is consistent with the observed drop in connected porosity (Fig. 5b).

The data indicate that there is no strong dependence of the $k(\phi)$ trend with temperature, or the thermal path taken to reach the peak temperature. This is consistent with Eqs. 1–4 which predict that, while temperature controls the time-evolution of permeability, it does not influence the path of the permeability–porosity relationship. Minor anisotropy in the permeability at the beginning of the experiment (e.g. at high porosity) appears to grow at low porosities (observe the divergence of $k_x$ with respect to $k_z$ as $\phi \to \phi_c$, in Fig. 5a), which could be attributed to minor effects of the compression-packing method of sample preparation. No anisotropy is predicted by the model approach (Section 2).

## 4.3. Time-evolution of permeability $k(t)$

With a validated constitutive law for $k(\phi)$ (Eqs. 4–5), we can combine our model for $k(\phi)$ with the vented bubble model for $\phi(t)$ (Eqs. 1–3), giving a method for calculating the time-dependence of permeability during welding $k(t)$. Fluid flow simulations through the segmented 3D datasets from our in situ experiments yield $k(t)$ data, which we plot in Fig. 6. The data show that permeability drops during welding, and drops more rapidly at higher temperature.

In Fig. 6 we also show our model results in which the model for $s(\phi)$ (Eq. 5; Fig. 3) is used to solve the model for $k(\phi)$ (Eq. 4) using $\phi_c = 0.03$ and $b = 4.4$ (see Section 4.2; Fig. 5). This solution for $k(\phi)$ is then mapped to each $t$ by first solving the dimensionless universal form of the vented bubble model (Eq. 2) and then rendering that solution dimensional via Eq. 3 and the constitutive equations for the material used here (Section 3.2). The validation of this model for $\phi(t)$ is presented in detail elsewhere for these same experiments (Wadsworth et al., 2019) and not restated here. The model and the data agree well. For comparison, we also show what the model solution would be if we did not account for the diffusive degassing of the particles during sintering, and had instead assumed either the initial or equilibrium H₂O concentration in the glass throughout (Fig. 6). Neglecting diffusive degassing of the particles provides a very poor description of the data

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** The relationship between porosity $\phi$ and permeability $k$ of welding systems. (a) The $k(\phi)$ data for the in situ experiments analyzed here. The solution to the $k(\phi)$ model using Eq. 4 is shown for $\phi_c = 0.03$ and $b = 4.4$ (Feng et al., 1987). The grey band represents the range of model solutions that encompass 2% uncertainty around the best fit $\phi_c$ and $b$ pair. Inset: the goodness of fit (fit in logarithmic space) given by a Deming $r^2$ (see text) of Eq. 4 for different $\phi_c$ and $b$ value pairs, for which the best-fit value is indicated with the dashed line. The region of $\phi_c$ and $b$ value pairs given by the grey band in (a) is enclosed by the white dashed line. (b) The relationship between the pore connectivity $\phi_p/\phi$ and $\phi$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The time-evolution of permeability for different isothermal dwell temperatures: (a-b) 1053 K, (c-d) 1144 K, (e-f) 1217 K, and (g-h) 1248 K; symbols as for Fig. 5. The top panels show the temperature-time paths; the bottom panels show the associated permeability-time paths. The range of the time axis shown is kept constant for comparison of the timescales involved for different temperatures. In the bottom panels we also show the solution to the *k(t)* model we develop in Section 2. The solid curves are the coupled model for which the syn-welding degassing of the particles is accounted for, while the dashed and dash-dot curves are for the case where the particles are assumed to be at the initial or equilibrium H₂O concentration throughout, respectively.</figcaption></figure>

and confirms that, for our experimental conditions, accurate prediction of permeability during welding requires that degassing is accounted for.

**5. Discussion**

The results in Section 4 constitute an experimental validation of the workflow for modelling *k(t)* via separate models for *φ(t)*, *s(φ)*, and *k(s, φ)*. The model accounts for syn-welding pyroclast degassing and the associated change in pyroclast viscosity, contains no empirical adjustments, and is grounded in percolation and sintering theory. This means that, once validated in a given regime (e.g. low *P*), it can be extrapolated and applied to conditions or compositions, within that regime, beyond those for which it was specifically validated. In this section, we test the extent to which we have indeed validated this model and state the conditions that would require further validation.

**5.1. Toward a general model for permeability and permeability evolution in welding systems**

The data collected here are for a specific composition, particle size distribution, and thermal treatment. In order to test the extent to

which our model results for *s(φ)*, *k(s, φ)*, and *k(t)* can be applied more generally, we can compare our model with other available data for those quantities. The data used for this comparison are summarized in Table 1. To compare across a range of material properties and experimental or simulation conditions, it is convenient to use the dimensionless form of our models.

The *φ(t)* model (Eq. 1) is already presented in dimensionless form $\bar{\phi}(\bar{t})$ (Eqs. 2 and 3) and has been validated extensively over a range of conditions (Wadsworth et al., 2016a, 2017b, 2017c, 2019).

We non-dimensionalize the specific surface area *s* by the moments of the pore size distribution (see Section 4.1), giving $\bar{s} = s(a^3)/(a^2)$ for the polydisperse case, and $s\langle a \rangle$ for the monodisperse case. The model (Eq. 5) reduces to the dimensionless form $\bar{s} = -3(1-\phi)\ln(1-\phi)$, which is shown in Fig. 7a, along with data from this study, literature data from sintering anhydrous spherical glass beads using the same in situ technique (Wadsworth et al., 2017b), and literature data from simulations of overlapping sphere packs (Vasseur and Wadsworth, 2017; Vasseur et al., 2020) where $\langle a^3 \rangle/\langle a^2 \rangle$ or $\langle a \rangle$ are computed following the same technique as used here (Torquato, 2013; Wadsworth et al., 2017c). Data are presented in Table 1. Despite the very different pore size and particle size distributions in these studies, the collapse to the
<table><caption>**Table 1**
Parameters used to normalize data.</caption><thead><tr><th rowspan="3">Composition</th><th rowspan="3"></th><th colspan="6">Normalization for s(φ)</th><th colspan="2">Normalization for k(φ)</th><th rowspan="2">Normalization for k(t)</th><th rowspan="2">Reference for data</th></tr><tr><th>φ<sub>i</sub></th><th>⟨R⟩</th><th>⟨R<sup>2</sup>⟩</th><th>⟨R<sup>3</sup>⟩</th><th>⟨a<sub>i</sub>⟩</th><th>⟨a<sub>i</sub><sup>2</sup>⟩</th><th>⟨a<sub>i</sub><sup>3</sup>⟩</th><th>φ<sub>c</sub></th><th>s(φ) from data?*</th></tr><tr><th>m</th><th>m<sup>2</sup></th><th>m<sup>3</sup></th><th>m</th><th>m<sup>2</sup></th><th>m<sup>3</sup></th><th></th><th>μ(T) model</th><th>Γ</th><th>N/m</th></tr></thead><tbody><tr><td>Rhyolite</td><td>Polydisperse</td><td>0.55</td><td>2.33 x 10<sup>-5</sup></td><td>7.66 x 10<sup>-10</sup></td><td>2.96 x 10<sup>-14</sup></td><td>8.01 x 10<sup>-6</sup></td><td>9.40 x 10<sup>-11</sup></td><td>1.35 x 10<sup>-15</sup></td><td>0.038</td><td>Yes</td><td>Hess and Dingwell (1996)</td><td>0.3</td><td>This study</td></tr><tr><td>Soda-lime-silica</td><td>Monodisperse</td><td>0.39</td><td>9.29 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>1.16 x 10<sup>-5</sup></td><td>2.33 x 10<sup>-10</sup></td><td>6.59 x 10<sup>-15</sup></td><td>0.020</td><td>Yes</td><td>–</td><td>–</td><td>Blair et al. (1993)</td></tr><tr><td>Soda-lime-silica</td><td>Monodisperse</td><td>0.46</td><td>4.09 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>1.88 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>0.013</td><td>Yes</td><td>Wadsworth et al. (2016a)</td><td>0.3</td><td>Wadsworth et al. (2017b)</td></tr><tr><td>Soda-lime-silica</td><td>Monodisperse</td><td>0.50</td><td>9.41 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>1.73 x 10<sup>-5</sup></td><td>4.45 x 10<sup>-10</sup></td><td>1.41 x 10<sup>-14</sup></td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Wadsworth et al. (2020a, 2020b)</td></tr><tr><td>Soda-lime-silica</td><td>Monodisperse</td><td>0.50</td><td>6.02 x 10<sup>-4</sup></td><td>3.64 x 10<sup>-7</sup></td><td>2.21 x 10<sup>-10</sup></td><td>1.75 x 10<sup>-4</sup></td><td>1.65 x 10<sup>-14</sup></td><td>7.11 x 10<sup>-31</sup></td><td>0.030</td><td>Yes</td><td>–</td><td>–</td><td>Eichheimer et al. (2020)</td></tr><tr><td>Rhyolite</td><td>Polydisperse</td><td>0.35</td><td>7.50 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>8.02 x 10<sup>-6</sup></td><td>9.78 x 10<sup>-11</sup></td><td>1.49 x 10<sup>-25</sup></td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Okumura and Sasaki (2014)</td></tr><tr><td>Steel and brass</td><td>Monodisperse</td><td>–</td><td>–</td><td>–</td><td>–</td><td>4.5 - 218 x 10<sup>-6</sup></td><td>–</td><td>–</td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>https://amespore.com/en/Vasseur et al. (2020)</td></tr><tr><td>–</td><td>Polydisperse</td><td>–</td><td>–</td><td>–</td><td>–</td><td>1.12 - 1.6 x 10<sup>-6</sup></td><td>1.28 - 3.76 x 10<sup>-10</sup></td><td>0.15 - 1.64 x 10<sup>-14</sup></td><td>–</td><td>Yes</td><td>–</td><td>–</td><td>Vasseur and Wadsworth (2017)</td></tr><tr><td>–</td><td>Monodisperse</td><td>–</td><td>1 - 2 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>0.030</td><td>Yes</td><td>–</td><td>–</td><td>Vasseur and Wadsworth (2017)</td></tr><tr><td>Dacite</td><td>Polydisperse</td><td>0.26</td><td>2.5 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>1.80 x 10<sup>-5</sup></td><td>5.01 x 10<sup>-12</sup></td><td>1.74 x 10<sup>-17</sup></td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Heap et al. (2015)</td></tr><tr><td>Dacite</td><td>Polydisperse</td><td>0.43</td><td>1 x 10<sup>-4</sup></td><td>–</td><td>–</td><td>1.46 x 10<sup>-5</sup></td><td>3.19 x 10<sup>-10</sup></td><td>8.64 x 10<sup>-15</sup></td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Kendrick et al. (2016)</td></tr><tr><td>Rhyolite</td><td>Polydisperse</td><td>0.50</td><td>2.5 - 15 x 10<sup>-6</sup></td><td>–</td><td>–</td><td>0.46 - 2.75 x 10<sup>-6</sup></td><td>0.03 - 1.13 x 10<sup>-11</sup></td><td>0.03 - 5.72 x 10<sup>-17</sup></td><td>0.030</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Heap et al. (2019)</td></tr><tr><td>Dacite</td><td>Polydisperse</td><td>0.50</td><td>1 x 10<sup>-5</sup></td><td>–</td><td>–</td><td>1.84 x 10<sup>-6</sup></td><td>5.02 x 10<sup>-12</sup></td><td>1.69 x 10<sup>-17</sup></td><td>–</td><td>No (monodisperse)</td><td>–</td><td>–</td><td>Ryan et al. (2020a, 2020b)</td></tr></tbody></table>

dimensionless $\bar{s}(\phi)$ model is generally reliable. We note that some samples from Wadsworth et al. (2017b) have $\bar{s}$ that deviate from the model, which is likely to be due to heterogeneity in the initial polydisperse packing. The discrepancy between our data and the model at high $\phi$ was noted in Section 4.1, and attributed to particle angularity. The literature data, which are for systems of spheres, show closer agreement at high $\phi$, which is consistent with that hypothesis. This additionally supports the conclusion that welding systems are indeed microstructural analogues of overlapping particle domains. Across all datasets, we find reasonable agreement with our model (Fig. 7a) with a global goodness of fit of $r^2 = 0.80$.

We non-dimensionalise the permeability $k$ by a reference permeability $k_r = 2(1 - (\phi - \phi_c)/s^2$, giving $\bar{k} = k/k_r$. The model (Eq. 4) reduces to the dimensionless form $\bar{k} = (\phi - \phi_c)^b$, which is shown in Fig. 7b, along with data from this study, and literature data from studies using samples which were (or are inferred to have been) initially particulate and which have undergone welding: sintered soda-lime-silica glass beads from Wadsworth et al. (2017b), which were also welded using the same in situ technique; glass beads welded ex situ (Blair et al., 1993; Wadsworth et al., 2020b; Eichheimer et al., 2020); obsidian welded ex situ (Okumura and Sasaki, 2014); welded crystal-bearing anes desire particles (Kendrick et al., 2016); and natural welded systems for which $s$ is constrained (Heap et al., 2015, 2019; Ryan et al., 2020a). Where the original study does not give a value for $s$, we calculate it using Eq. 5. Data are presented in Table 1. Across this range of experimental-to-natural variability, we find reasonable agreement with our model (Fig. 7b) with a global goodness of fit of $r^2 = 0.94$.

Finally, in Fig. 7c, we test our $k(t)$ model against previous data. We compare the results from this study with those from Wadsworth

et al. (2017a) in which soda-lime-silica glass beads were sintering in situ, yielding a similar time-dependent description of the permeability. Due to subtle differences in the initial packing and the random close packing difference between spherical and angular particles, the initial value of $\bar{k}$, termed $\bar{k}_i$, from which each specific experiment starts the process is different. For any given $\bar{k}_i$, there will be a different model result in $\bar{k}(t)$. Therefore, to provide a common basis that allows us to compare data and model results across all conditions, we use $\bar{k}/\bar{k}_i$, which we note is equivalent to $k/k_i$ (i.e., the $k_r$ terms cancel). Note that this normalization step does not affect the generality of Eq. 4, but is simply a device to allow all of the data to be compared to a single, master curve.

The dimensionless time $\bar{t}$ is defined in Eq. 3 and is specific to the temperature-time pathway taken. For the anhydrous soda-lime-silica glass beads, there is not the additional requirement to account for the degassing of $H_2O$, which is accounted for in the sintering obsidian particles analyzed herein (Wadsworth et al., 2019). When cast in this space (Fig. 7c), the model given in Eq. 4 becomes $\bar{k}/\bar{k}_i = [(\phi - \phi_c)/(\phi_i - \phi_c)]^b$ and, following the same principles as is shown dimensionally in Fig. 6, the time dependence for $\bar{\phi}(\bar{t})$ is used together with $k(\phi)$ to predict $\bar{k}(\bar{t})$. For the sintering of soda-lime-silica glass beads, we use the monodisperse form of this model, because the beads were sieved to a single size fraction, while for the data presented herein, we use the polydisperse model presented here. This difference results in two different model curves. The agreement between the models and these two datasets leads us to conclude that, while the particle angularity played a minor but measurable role in controlling $s(\phi)$, it does not strongly influence the dynamics of permeability reduction. Similarly, this suggests

<footnote>a If 'yes' then s was a measured quantity, if 'no' then s is computed from <Rn>.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** The relationships from Figs. (4–6) rendered dimensionless and compared with available published data. In all panels the data presented in colour are the experimental data from this study (refer to panel (c) for the temperatures for each colour). The normalization values for each published dataset are given in Table 1. (a) The dimensionless specific surface area as a function of porosity $\bar{s}(\phi)$. Additional data shown are for glass bead packs sintered in situ (Wadsworth et al., 2017b) and ex situ (Eichheimer et al., 2020), and for sphere pack simulations (Vasseur et al., 2020). (b) The dimensionless permeability as a function of the scaled porosity $k(\phi - \phi_c)$. Additional experimental data shown are for glass bead packs sintered in situ (Wadsworth et al., 2017b) and ex situ (Blair et al., 1993; Wadsworth et al., 2020b; Eichheimer et al., 2020), for experimentally welded obsidian particles (Okumura and Sasaki, 2014), and for commercially available AmesPore® sintered metal filters (data from: https://amespore.com/en/). Natural data shown are for welded crystal-poor (Heap et al., 2019) and crystal-rich (Kendrick et al., 2016) tuffisites, welded block-and-ash flow deposits (Heap et al., 2015), and dacitic volcanic fault gouge (Ryan et al., 2020a). (c) The dimensionless permeability as a function of dimensionless time $\bar{k}(t)$. The curves represent the time-dependent model given in Eq. 2 for $\bar{\phi}(\bar{t})$ coupled with the $\bar{k}(\phi)$ model. The difference between the solid and the dashed curves arises from the definition of $a_i$ for monodisperse (dashed curve) and polydisperse (solid curve) sintering particles. We compare our data with previous data for monodisperse sintering glass beads (Wadsworth et al., 2017b).</figcaption></figure>

that degassing and the associated evolution of viscosity gradient in the particles does not substantially alter the evolution of the microstructure and so does not appear to impact the permeability decay when accounted for via the approach of Wadsworth et al. (2019).

## 5.2. Successes and limitations of our model and validation

The in situ nature of the data allows us to determine $s$, $\phi$, and $k$, independently, as functions of measured experimental time $t$; hence to validate, independently, our models for $s(\phi)$ (Fig. 4), $k(\phi)$ (Fig. 5), and $k(t)$ (Fig. 6). In each case we find reasonable agreement between model and experiment.

The principal limitation to the model validation is associated with the parameter $\bar{P}$ (Eq. 3). Wadsworth et al. (2019) showed that Eqs. 1–3 are valid descriptions of $\phi(t)$ during welding, even when confining pressure exceeds the Laplace surface pressure – i.e. when $\bar{P}>1$ and validated up to $\bar{P} = 100$. However, this does not necessarily confirm that the combination of $\phi(t)$ and $k(\phi)$, to determine $k(t)$, would also be valid at higher confining pressures. It could be that the so-called pressure sintering regime at $\bar{P}>1$, induces particle-particle flattening sufficient to change the microstructure such that the analogy with static percolation theory developed for overlapping spheres, which underpins our $k(\phi)$ model, is no longer applicable. For instance, as $\bar{P}$ increases above unity, the sintering rate increases (Rahaman and De Jonghe, 1990; Ryan et al., 2018; Wadsworth et al., 2019), and it would be reasonable to assume that the microstructural evolution would be different as slow, capillary-driven processes, such as neck formation and pore geometry smoothing, become negligible on the timescale of porosity decay. We therefore consider that our model is validated only for $0<\bar{P}<1$, and validation of $k(\phi, t)$ at higher $\bar{P}$ should be a topic for future investigation. Nevertheless, the low $\bar{P}$ regime is relevant to a range of volcanic scenarios (Gardner et al., 2018, 2019; Wadsworth et al., 2019; Heap et al., 2019). Similarly, we have not provided a tensorial description of Eq. 4 that would allow permeability anisotropy to be predicted under stresses that give rise to shear deformation. However, syn-welding shear is clearly preserved by textures in volcanic deposits (Fig. 1; c.f. Andrews and Branney, 2011), so the development of anisotropy of permeability during welding under shear is another priority for future investigation.

There are further complexities that may arise during welding in nature that are not accounted for in our model. For example, in Fig. 1 we show that clasts in tuffisites or ignimbrites may be internally porous (e.g. pumice and fiammé), or variably crystalline (e.g. lithics). Additionally, welding packs may be heterogeneous, composed of particles with a range of compositions, representing a range of provenance pre-deposition (Saubin et al., 2016) or variable alteration during transport or gas flow (Berlo et al., 2013; Casas et al., 2019; Paisley et al., 2019), and the pore space can be in-filled with vapor-phase deposition or expanded by leaching of pore walls (Schipper et al., 2015). This kind of diversity of clast types and evolution of pore geometries in natural systems could be accounted for in a model that convolves Eqs. 1–3 with a distribution of clast rheology, which in turn could be rendered time-dependent in a similar manner to our quantitative account of degassing which leads to a time-dependent particle rheology. However, we have not validated those effects here, nor has the impact of these effects on the permeability been tested, and it is possible that a wide distribution of particle rheology has a non-trivial impact on the $k(t)$ pathway taken during welding.

As a coarse test of the impact of the complexities discussed above, we show that our $k(\phi)$ model is effective across a wide range of natural samples collected from systems that sintered under a wide range of conditions (Fig. 7b). For example, the sintered tuffisite materials from Heap et al. (2019) and the volcanic gouge sintered in the solid-state from Ryan et al. (2020a, 2020b) are both shown to have formed at appreciable confining pressure. Heap et al. (2019) show that the wall of their studied tuffisites host an equilibrium value of 0.46 wt.% H₂O, while the far-field hosts an equilibrium value of 1.04 wt.% H₂O, and that these limits are bridged by a diffusion profile toward the tuffisite. Assuming
that the wall value would be in equilibrium with the gas at the vapor pressure in the tuffisite (and we assume that this value was approximately constant), and that the far-field value would be in equilibrium with the liquid at magmastatic pressure, we can estimate P̄ for this tuffisite. Using a solubility model (Liu et al., 2005) we find that the wall value relates to P<sub>g</sub> ≈ 2 MPa and that the far-field value relates to P ≈ 9 MPa, resulting in ΔP ≈ 7 MPa. Given the constraints on P<sub>L</sub> from Heap et al. (2019), we therefore estimate that P̄≈116 for that tuffisite. Despite the high P̄ regime and the polymictic nature of the tuffisite fill, the permeability–porosity relationship collapses to our general description to within less than an order of magnitude deviation in *k* from our model at any porosity (Fig. 7b). We propose that such a collapse in *k* (φ) space is a reasonable indication that, while the effects discussed above may influence pore geometries, this does not appear to translate into a substantial difference in the hydraulic properties of these rocks.

# 6. Applications and regimes in natural welding

## 6.1. Degassing regimes in natural rhyolite welding

We can apply our experimentally validated model to investigate volcanic welding with syn-welding degassing. Here, we investigate welding in a rhyolite ash pack with an initial water concentration Cᵢ, initial monodisperse particle size R, isothermal temperature T, and initial packing porosity φᵢ, using our model to predict the evolution of permeability with time, and the time it takes for the permeability reach zero (i.e. the time at which φ < φᵢ). We explore the scenario where T = 1100 K, which represents an insulated magmatic temperature (cooling is discussed later), and an initial porosity φᵢ = 0.5, which is typical of loose randomly packed angular particles (e.g. Fig. 5). We vary Cᵢ and R over a wide range representative of silicic systems: 0.1 < Cᵢ < 5 wt.%, and 10⁻⁶ < R < R_c m, where R_c = √Γ/(ρg) is the radius above which gravitational effects dominate the particle

deformation over surface tension effects, ρ ≈ 2000 kg.m⁻³ is the rhyolite particle density (Lange and Carmichael, 1987), and g is the acceleration due to gravity. For rhyolites, R_c ≈ 4 × 10⁻³ m, such that particles smaller than this value will not deform under their own mass. We assume that welding is occurring at P_g = 0.1 MPa vapor pressure at the Earth's surface, simulating welding in very shallow tuffisites or at the base of hot pyroclastic density currents. We confine our analysis to the low-P̄ regime (discussed above). At these conditions, the solubility model we used for our experimental validation predicts that the equilibrium water concentration is C_e = 0.1 wt.% (Liu et al., 2005), which, given sufficient time, is the value to which the particles will degas by diffusion.

In Fig. 8, we show indicative model results. We find that when the initial water concentration Cᵢ is close to or equal to the equilibrium value (e.g. initial supersaturation is low or zero), the permeability during welding decays smoothly. In this case packs with larger particle radius will both (1) start from a higher initial permeability, and (2) take longer to weld to zero permeability (Fig. 8a). When the initial water concentration is much higher than the equilibrium value, such that Cᵢ ≫ C_e, we find more complex welding pathways can be taken in k(t) space, depending on particle radius. We find three phenomenological regimes; to demonstrate these regimes, we use the example of Cᵢ = 1 wt. % (Fig. 8d).

(1) For relatively small particles (R = 10⁻⁵ m in this example) permeability decay is almost indistinguishable from the case where Cᵢ = C_e (Fig. 8a). Inspection of the evolution of the spatially-averaged H₂O concentration 〈C〉(t) (Fig. 8e) demonstrates that this is because 〈C〉 falls rapidly to C_e, before appreciable welding has occurred: i.e. degassing is efficient on the timescale of welding. Consequently, the particles have equilibrium H₂O concentration for the majority of the welding process, and the final outcome is expected to be a deposit that is fully degassed when welding is complete.

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Model results for permeability as a function of time k(t) and the average particle H₂O concentration as a function of time 〈C〉(t) for different particle radius R and initial water concentration Cᵢ; (a-b) for Cᵢ = C_e = 0.1 wt.%; (d-e) for Cᵢ = 1 wt.% and for which C_e is determined at 0.1 MPa. In (b) and (e), the filled circle against which the curves terminate represents the time at which welding completes λ'. All results are using the μ(T,C) relationship for rhyolites (Hess and Dingwell, 1996) and at T = 1100 K. c) Schematic cartoon showing the simultaneous evolution of pore-space geometry and H₂O distribution in the particle pack, as welding and degassing proceed.</figcaption></figure>
(2) For relatively large particles ($R = 10^{-3}$ m in this example) permeability decays much more rapidly than for the case where $C_i = C_e$, and more rapidly even than the smaller particles discussed in (1) above. Inspection of $\langle C \rangle(t)$ demonstrates that this is because H₂O concentration remains close to $C_i$ throughout most of the welding process: i.e. degassing is inefficient on the timescale of welding. Consequently, the viscosity of the particles is much lower than for the equilibrium case, leading to faster welding. The final outcome is expected to be a deposit that retains disequilibrium H₂O concentration when welding is complete.

(3) For particles of intermediate size ($R = 10^{-4}$ m in this example) the permeability decay is more complex, with a pronounced inflection during welding. This is because degassing and welding operate over similar timescales, and the increase in viscosity that occurs as the particles degas leads to an appreciable decrease in welding rate as welding progresses. Depending on the exact conditions, the final deposit will have a water content between $C_i$ and $C_e$ when welding is complete.

To explore these regimes further, we consider the total time taken for welding to complete $\lambda'$ for particles of different size and initial H₂O concentration. We take $\lambda'$ as the time elapsed when $\phi$ reaches $\phi_c$. In Fig. 9 we plot $\lambda'$ as a function of the starting particle radius, showing the solution for a wide range $C_e \le C_i \le 5$ wt.%, for rhyolitic particles welding at 1100 K and $P_g = 0.1$ MPa. We observe that, for a given $C_i$, the welding time changes with particle size, and passes through three regimes, corresponding to regimes (1)–(3) above. In regimes (1) and (2), welding time increases with particle size, but welding time in regime (2) may be lower than welding time in regime (1), despite larger particle size, if H₂O concentration is sufficiently high. In regime (3), at intermediate particle size, welding time decreases as particle size increases.

These 3 regimes can be understood by considering the capillary Peclet number $P_c$ (Gardner et al., 2018), which compares the timescale for diffusive degassing $\lambda_D$ with the timescale for welding in the low-$\bar{P}$ regime, $\lambda$ (Eq. 3). The diffusive timescale is $\lambda_D \approx R^2/D$. For convenience, we can define $P_c$ for the initial conditions,

$$P_c = \frac{\lambda_D}{\lambda} = \frac{\Gamma R^2}{D_i \mu_i \langle a_i \rangle} \qquad \text{Eq. 7}$$

<figure><img src="image_12.png" /><figcaption>Fig. 9. The welding time $\lambda'$ as a function of the radius $R$ of the welding particles. Model curves are plotted for a wide range of initial H₂O concentration 0.1 < $C_i$ < 5 wt.%. Models are plotted only for $R < R_c$, such that gravitational effects are unimportant. The low-Pc solution for $\lambda'(R)$ is shown in the dashed line, and the solid line connects the $P_c = 1$ point for each $C_i$, (i.e. the point where $\lambda = \lambda_D$).</figcaption></figure>

where $D_i$ and $\mu_i$ are the properties $D$ and $\mu$ at the initial H₂O concentration prior to diffusion. If $P_c \ll 1$, then degassing is rapid compared with welding, and $C_i \to C_e$ early in the welding process. This manifests as the viscosity moving rapidly from an initial value $\mu_i$ determined by $C_i$, toward a final value $\mu_e$ determined by $C_e$. However, if $P_c \gg 1$, then welding is rapid compared with degassing, and $C_i$ does not substantially change during the welding process. This manifests as a near-constant viscosity at $\mu_i$ determined by $C_i$. Then, in the limiting cases of $P_c \ll 1$ and $P_c \gg 1$, we expect the time at which the permeability reaches to zero to be approximated by the welding time such that

$$\lambda' = \frac{\mu_e \langle a_i \rangle}{\Gamma} \approx \frac{\mu_e R}{\Gamma} \qquad \text{Pc} \ll 1 \qquad \text{Eq. 8a}$$

$$\lambda' = \frac{\mu_i \langle a_i \rangle}{\Gamma} \approx \frac{\mu_i R}{\Gamma} \qquad \text{Pc} \gg 1 \qquad \text{Eq. 8b}$$

and for $P_c \ll 1$, we expect that the full model is required to predict the welding time, with no simple limiting solution available. In Eq. 8 we use the fact that $\langle a_i \rangle$ and $R$ are typically of the same order of magnitude to make the approximation $\langle a_i \rangle \approx R$, rendering $\lambda'$ easy to compute in these limiting cases. In Fig. 9, we plot $\lambda' = \mu_e R/\Gamma$ as the dashed curve (Eq. 8a), and find a minor offset between the limit that is approached by the solutions and this approximation, which relates to the use of $R$ in place of $\langle a_i \rangle$. The limit at $P_c \gg 1$ (Eq. 8b) is not a single curve and depends on $C_i$ (Fig. 9).

The regime boundary between $P_c \ll 1$ and $P_c \gg 1$ is not necessarily precisely at $P_c = 1$. To test this, we find the point at which $\lambda = \lambda_D$ for each $C_i$, which, by definition, is the point at which $P_c = 1$, and plot this continuously in Fig. 9. The $P_c = 1$ curve cuts the $\lambda'(R)$ curves within their transitional zones (i.e. the inflected portion) in the majority of cases, indicating that $P_c = 1$ is a reasonable approximation of the regime boundary. The transitional zone is more sharply defined for high $C_i$ that for low $C_i$, and the absolute difference between the low and high $P_c$ solutions is greater high $C_i$.

## 6.2. Cooling during welding and the preservation of permeable networks

In volcanic environments, ignimbrites are emplaced onto a relatively low-temperature substrate, cooling dynamically during emplacement (Branney and Kokelaar, 2002). Similarly, while tuffisites are sourced in the magmatic conduit itself, they are often emplaced in hydraulic fractures in cold country rock (Stasiuk et al., 1996). In both cases, dynamic cooling occurs and even simple models show that this can result in incomplete welding if the material quenches before the welding time is reached (Wadsworth et al., 2014; Kolzenburg et al., 2019). Indeed, the value that the welding time takes (e.g. Eq. 3; c.f. Eq. 8) is a function of time as the temperature drops.

We can apply the approach taken in Section 6.1 for isothermal conditions, to investigate welding permeability under non-isothermal conditions of continuous cooling. During cooling, the limiting temperature below which we assume the system is no longer able to weld viscously is the glass transition temperature $T_g$. While there is evidence that $T_g$ depends on the cooling rate, – termed $q$, as $T \to T_g$ (Gottsmann et al., 2002), we make the simplifying assumption that $T_g$ is the temperature at which $\mu = 10^{12}$ Pa.s (Gottsmann et al., 2002; Wadsworth et al., 2017a). Using Eq. 3, our full model can be mapped to any non-isothermal path. However, for simplicity and illustrative purposes, we assume a linear cooling rate $q$ over the interval between an initial temperature $T_i$ and $T_g$, which is usually a few hundred kelvin, and therefore relatively small compared with the total cooling path from $T_i$ to ambient conditions. During linear cooling, the viscosity of the particles increases and the diffusivity of H₂O decreases, causing a non-trivial competition between continued degassing of the particles and the welding rate.

In Fig. 10 we show the results of our model runs for non-isothermal environments, which we cast as the permeability preserved when $T = T_g$. At the fastest cooling rates, little welding occurs before $T_g$ is reached, and the initial permeability is preserved. At the slowest cooling
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** The permeability that is preserved in a welding deposit when temperature falls below the glass transition temperature (i.e. when $T = T_g$) as a function of the cooling rate $q$. There is a critical cooling rate $q'$, below which the permeability falls to zero before $T_g$ is reached (vertical arrows). Results are shown for (a) low $C_i = 0.1$ wt.% and (b) relatively high $C_i = 1$ wt.%. In all cases the equilibrium H₂O concentration is $C_e = 0.1$ wt.%.</figcaption></figure>

rates, welding goes to completion, and permeability drops to zero, by the time $T_g$ is reached. At intermediate cooling rates, welding is arrested before completion, and the final permeability has a value intermediate between the initial value and zero. When $C_i = C_e$, hence no degassing occurs, we can see that the final permeability is a simple, smooth function of $q$, whereas when $C_i > C_e$, the final permeability may be a more complex function of $q$.

There is a critical $q$, which we term $q'$, below which the permeability will always fall to zero before $T = T_g$ (Fig. 10). We take $q'$ as a cooling rate of first-order interest, as it separates the regime of complete welding from the regime of incomplete welding. In Fig. 11 we show the critical $q'$ as a function of $R$ for the same range of $C_i$ as given for isothermal conditions in Fig. 9. This shows that, for relatively small particles, the critical cooling rate collapses to a limiting $q'(R)$ for any $C_i$; for relatively large particles, $q'(R)$ evolves toward a second limiting solution at relatively higher $q'$. Following the analysis in Section 6.1, we find that the limiting solution $q'(R)$ at small $R$ is the solution for the $Pc \ll 1$. Conversely, the high $q'(R)$ limit at large $R$ is the solution for the $Pc \gg 1$. This implies that there are again three regimes for a given $C_i$

(1) At small $R$ and $q < q'$, welding completes and the resulting deposit is degassed to equilibrium H₂O concentration.

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** The critical cooling rate $q'$ above which welding will not complete and $k > 0$ will be preserved in deposits, as a function of $R$. Model results are shown in colour for a wide range of $C_i$ (see Fig. 9). The low-Pc solution for $q'(R)$ is shown in the dashed line, and the solid line connects the $Pc = 1$ point for each $C_i$.</figcaption></figure>

(2) At large $R$ and $q < q'$, welding completes but the resulting deposit retains disequilibrium H₂O concentration.

(3) At any $R$ and $q > q'$, welding is incomplete, $k > 0$, and $t/\lambda_D$ can be used to determine if the deposit is degassed or not, independent of the welding model.

In regimes (1) and (2), we can find an analytical solution for the critical cooling rate as

$$q' = \frac{T_i - T_g}{\lambda'} = \frac{(T_i - T_g)\Gamma}{\mu_e \langle a_i \rangle} \approx \frac{(T_i - T_g)\Gamma}{\mu_e R} \quad \text{Pc} \ll 1 \qquad \text{Eq.9a}$$

$$q' = \frac{T_i - T_g}{\lambda'} = \frac{(T_i - T_g)\Gamma}{\mu_i \langle a_i \rangle} \approx \frac{(T_i - T_g)\Gamma}{\mu_i R} \quad \text{Pc} \gg 1 \qquad \text{Eq.9b}$$

and as before these regimes are separated by $Pc = 1$ (solid curve on Fig. 11).

### 6.3. Direct volcanic applications

Welding occurs in volcanic eruptions in distinct scenarios: on the walls of open conduits during Plinian eruptions (Gardner et al., 2017, 2019); in the hot components of jet engines (Song et al., 2014; Giehl et al., 2016), in tuffisites (Tuffen and Dingwell, 2005; Castro et al., 2012; Heap et al., 2019); at the depositional base of moving pyroclastic density currents (Branney and Kokelaar, 1992, 2002); and during lightning strikes in-air or air-to-ground (Cimarelli et al., 2017; Mueller et al., 2018). While the dynamics of welding is important in controlling the bulk dynamics in all of these scenarios (Wadsworth et al., 2019), we focus on tuffisites and the sedimentation of ignimbrites as the two scenarios where permeability evolution may be most important for understanding outgassing and pressure evolution. Here, we use our model to predict whether a cooling ignimbrite or tuffisite will complete welding before cooling finishes, and whether the final deposit will preserve equilibrium or disequilibrium H₂O.

As discussed in Section 6.2, ignimbrites are sedimented onto a substrate that is usually at a lower temperature than the particles deposited; therefore, any initial welding will be strongly non-isothermal and indeed the quench against the substrate may inhibit welding. However, as an ignimbrite accumulates, the cooling rate at any given point may decrease. Depending on the sedimentation rate, the cooling rate may end up being low, especially far from the substrate (Gottsmann and Dingwell, 2001). Therefore, the average cooling rate from deposition temperature, $T_i$ to $T_g$, may vary in time and space in welding ignimbrites. By contrast, the particles in tuffisites are emplaced sub-surface
and generally involve smaller transport distances; they may therefore be expected to have a higher emplacement temperature compared with ignimbrites. For tuffisites, from a cooling perspective, there are two main types: (1) those that open into hydraulic fractures in visco-elastic magma (Tuffen and Dingwell, 2005) and (2) those that open into hydraulic fractures in country rock surrounding the conduit (Stasiuk et al., 1996). In both cases, the expected pressure drop due to fracture propagation may induce adiabatic cooling, but in general, the cooling in case (1) can be expected to be much less than in case (2) (Kolzenburg et al., 2019) and in general a single isothermal temperature is usually estimated when considering the welding dynamics of tuffisites in regime (1) (Castro et al., 2012; Wadsworth et al., 2019; Heap et al., 2019).

Any real cooling path will be non-linear and calculation of $T(t)$ requires solution of Fourier's law for conduction. Given that this is geometry and condition dependent, we use a scaling approach, and estimate average cooling rate using $\langle q \rangle \approx (T_i - T_g)D_T/L^2$ where $D_T$ is the thermal diffusivity and $L$ is the system lengthscale in the dominant direction of heat flow. Using data from Bagdassarov and Dingwell (1994) for rhyolitic melts, we can approximate $D_T$ as independent of temperature over the range $700 \le T \le 1200$ K and take $D_T = 1 \times 10^{-6} \text{ m}^2/\text{s}$. For ignimbrites, we assume that the lower limit on $L$ is the depositional boundary layer thickness, estimated to be on the order of 1 m (Andrews and Branney, 2011) and the upper limit is taken to be 10 m as a reasonable estimate of half-thickness of an entire ignimbrite. For tuffisites, we take the lower limit on $L$ to be 1 mm and the upper limit to be 10 cm based on available observations (Castro et al., 2012; Saubin et al., 2016; Heap et al., 2019). In all cases we take $T_i$ to be 1100 K. This results in $10^{-6} \le \langle q \rangle \le 10^{-4}$ K. s$^{-1}$ for ignimbrites, consistent with geospeedometry measurements from glassy welded ignimbrites (Gottsmann and Dingwell, 2001), and $10^{-2} \le \langle q \rangle \le 10^2$ K. s$^{-1}$ for tuffisites injected into country rock. For tuffisites injected into viscoelastic magma, the cooling rate will be significantly lower and cooling may not occur. For ignimbrites, while there are clearly larger particles involved (e.g.

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Application of our model to natural conditions (see Fig. 1). The model curves show the critical cooling rate $q'$ as a function of particle radius $R$ for $C_i = 0.7$ and 0.9 wt.%, typical of the far-field value found in obsidian close to tuffisites (Castro et al., 2012; Saubin et al., 2016) and consistent with the upper value estimated for natural welded ignimbrites (Sumner and Branney, 2002). For a given curve, points below the curve represent conditions under which welding completes before the deposit quenches; for points above the curve, welding does not complete before the deposit quenches. The box shown for tuffisites is bounded by $10^{-6} \le R \le 10^{-3}$ m and by $q'$ scaled for tuffisite widths $0.1 \le L \le 10$ cm. The box for ignimbrites is bounded by $10^{-5} \le R \le 10^{-3}$ m and by $q'$ scaled for thickness from a lower limit of the rising deposition boundary layer in ignimbrites $L = 1$ m to an ignimbrite deposit thickness $L = 10$ m. Tuffisites are expected to quench before welding completes, whereas ignimbrites may weld fully before quenching.</figcaption></figure>

fiammé), the dominant particle size range is $10^{-5} \le R \le 10^{-3}$ m (Branney et al., 1992; Streck and Grunder, 1995; Andrews and Branney, 2011) and a maximum initial H₂O concentration estimate is $C_i = 0.87$ wt. % (Sumner and Branney, 2002). For tuffisites, the particle size range is similar $10^{-6} \le R \le 10^{-3}$ m and $0.7 < C_i < 0.9$ wt. % (Castro et al., 2012; Saubin et al., 2016).

In Fig. 12 we use the estimates given above to constrain welding and degassing regimes for tuffisites and ignimbrites, plotting the critical cooling rate $q'$ as a function of radius $R$ of the welding particles (as in Fig. 11). Tuffisites plot well above the curves, indicating that they are always expected to quench before they weld under the conditions investigated. This implies that they would remain permeable and potentially act as efficient propped-fractures and outgassing pathways for magnetic gases from the conduit. Ignimbrites can straddle the regimes where welding will or will not occur on the timescale of cooling, consistent with the observation of welded and unwelded ignimbrites in nature. Moreover, for ignimbrites dominated by particles with radius $R \le 2.5 \times 10^{-4}$ m, for which the cooling is occurring slowly (e.g. fast accumulation rate), our model predicts that the deposit will be degassed and will weld over the approximate timescale given in Eq. 8b. For ignimbrites dominated by larger particles $R \ge 2.5 \times 10^{-4}$ m, our model predicts that welding will occur even for high cooling rates (e.g. slow accumulation rates or basal quenched contacts), but that the deposit will preserve disequilibrium H₂O concentrations.

Both the ignimbrite and tuffisite scenarios investigated in Fig. 12 are in the low-$\bar{P}$ regime, where confining pressure is small compared with the Laplace surface pressure (Eq. 3). We note that if the accumulation rate of an ignimbrite is particularly fast, or if the tuffisite is injected at depth, then welding will instead occur in the high-$\bar{P}$ regime. An increase in $\bar{P}$ moves the $q'(R)$ curves upwards (to faster cooling rates) compared with their position in Fig. 12. In that regime, welding in ignimbrites is always expected to occur even under the fastest cooling rates used here, and tuffisites may move into the regime where welding goes to completion before quench. For example, if we assume that the pressure is dominantly lithostatic with a mean country rock density of $\rho = 2000$ kg. m$^{-3}$ and the $P_g \approx 2$ MPa described earlier, we can estimate the depth below which even rapidly cooled tuffisites will weld shut on the timescale of cooling. This depth for a 1 mm tuffisite is 23.2 km and for a 10 cm tuffisite is 475.3 m.

Referring back to Fig. 9, isothermal tuffisites injected in hot magma are likely to weld rapidly (Wadsworth et al., 2019) and to record degassed H₂O contents in the particle fill. This is consistent with the densely welded nature of tuffisites observed in situ in Icelandic rhyolitic conduits (Tuffen et al., 2003; Tuffen and Dingwell, 2005). Tuffisites in obsidian bombs undergo cooling at rates intermediate between those injected into country rock (Fig. 12) and those welding near-isothermally in conduit interiors (Fig. 9), and likely cool at rates proportional to the square of the bomb size. The fact that many bombs preserve tuffisites that are incompletely welded (Heap et al., 2019) suggests that this bomb cooling rate is greater than $q'$. However, we note that a model prediction of the end state of permeability in a bomb would require knowledge of the timing of formation of the tuffisite relative to the time it was ejected as a bomb, as well as the subsequent cooling path of the bomb – information that at present is not available for any known example. Nevertheless, our model would predict that these tuffisites span the range between permeable and impermeable, consistent with observations (Castro et al., 2012; Saubin et al., 2016; Heap et al., 2019).

## 6.4. Outlook: permeability-limited welding and compaction

In a welding deposit, there is a competition between the timescale over which welding occurs $\lambda$ (in the low-$\bar{P}$ regime), and the timescale over which gas escapes from the welding material $\lambda_o$, which is a Darcy timescale that depends on the permeability and the length of the system (Michaut
et al., 2009; Wadsworth et al., 2016a; Kennedy et al., 2016; Heap and Wadsworth, 2016). If we assume the stress driving the flow of gas out from between the particles is $P_L$, then $\lambda_0 = \mu_g L_0^2 \langle a \rangle / (k\Gamma)$ where $L_0$ is the distance the gas has to travel to escape the welding system. Using $\lambda_0$ we can define a Darcy number as $\lambda_0/\lambda$ (Wadsworth et al., 2016a)

$$Da = \frac{\lambda_0}{\lambda} = \frac{\mu_g L_0^2}{k\mu} \qquad \text{Eq.10}$$

such that at Da << 1, outgassing is efficient and welding will complete in the way described by our model. However for Da >> 1, the welding is rapid compared with outgassing, and therefore the gas cannot be driven from the welding system sufficiently quickly, and gas pressure may rise in response to the Laplace surface pressure that is driving welding. In turn, this can lead to non-linear welding behavior in a regime termed 'compaction' by Wadsworth et al. (2019) on the basis that it is expected to be governed by non-linear compaction dynamics described elsewhere (Michaut et al., 2009). While Wadsworth et al. (2019) showed that most volcanic systems are initially in the Da << 1 regime, our analysis shows that many volcanic systems weld to completion and permeability drops to zero. Since Da → ∞ as k → 0 (Eq. 10); it's clear that Da >> 1 must be met during welding of such systems before welding completes.

The implication of Eq. 10 is that, at relatively high porosity in a highly permeable system, the gas phase is truly geochemically open. However, close to the percolation threshold, as the system becomes permeability limited, the gas phase could behave as if it were geochemically closed even in a permeable system. This is qualitatively consistent with the interpretation that the hydrogen-isotope systematics in a range of partially degassed tuffisites show evidence of batched open-closed system behavior (Castro et al., 2014). This regime remains underexplored.

## 7. Using this model: A practical guide

We note that the full model presented here is general, such that any pressure-temperature pathway could be incorporated, and the equations solved for $\phi(t)$ and $k(t)$. In this section, we present a short guide for navigating the solutions to our model framework that may be useful in various scenarios of practical interest.

To use this model as presented, we provide a VHub resource 'VolcWeld' (available at https://vhub.org/resources/4568) and a simplified version implemented in Excel™ is provided as Supplementary Material. The full model provided in the 'VolcWeld' program solves the full system of equations presented herein and therefore solves for the diffusive movement of H₂O. In turn, this diffusion solution is used to solve the 'vented bubble model' for welding and output $\phi(t)$ and $k(t)$. The Excel™ implementation (Supplementary Information), neglects the diffusion of H₂O, and therefore is suited to conditions where H₂O movement can be neglected on the timescale of welding (Pc << 1 and Pe >> 1; see Section 6), but still solves for the 'vented bubble model' using the Euler method for solving ordinary differential equations. The following steps can be followed to determine which solution tool is required:

Initial conditions:

(1) Constrain the particle size distribution $f(R)$ and find the mean particle size $\langle R \rangle$. Both the VolcWeld program and the Excel™ implementation (Supplementary Information) use Wadsworth et al. (2017c) to convert $f(R)$ into a pore size distribution $f(a_i)$ and find the mean pore size $\langle a_i \rangle$. If no information about the full distribution of particle sizes is known, then take an estimate for the typical particle size and define this as $\langle R \rangle$.

(2) Measure the initial H₂O concentration of the particles and determine whether or not this is in equilibrium by calculating the H₂O solubility using a constitutive law (e.g. Liu et al., 2005). Estimate the gas and liquid pressures and initial temperature of the system.

Isothermal welding

(3) If the system is in H₂O equilibrium and is isothermal, then no diffusion will occur and our Excel™ implementation can be used (Supplementary Information). In this situation, for scaling purposes, one can compute the melt viscosity $\mu_i$ (e.g. Hess and Dingwell, 1996), use $\Gamma = 0.3$ N. m⁻¹, and use $\langle a_i \rangle$ from step (1) to compute a welding timescale.

(4) If the system is not in H₂O equilibrium (see step (2)), then calculate the H₂O diffusivity $D_i$ (e.g. Zhang and Ni, 2010) and melt viscosity $\mu_i$ (e.g. Hess and Dingwell, 1996) under initial conditions. Using $\Gamma = 0.3$ N. m⁻¹, calculate the capillary Peclet number $P_c = \Gamma R^2 / (D_i \mu_i \langle a_i \rangle)$ for these initial conditions. If $P_c \le 0.1$ or $P_c \ge 10$, use the initial or equilibrium H₂O concentration to compute the melt viscosity, respectively. In either case, our Excel™ implementation can be used (Supplementary Information). If instead $0.1 < P_c < 10$, then the diffusive movement of H₂O renders the problem more complex and the VolcWeld program can be used (available at Vhub).

(5) For non-isothermal cooling conditions:

As with isothermal conditions, the Excel™ implementation can be used for $P_c \le 0.1$ or $P_c \ge 10$, where the first sheet in the implementation permits input of a temperature rate $q$. If, however, $0.1 < P_c < 10$, then the VolcWeld program must be used.

## 8. Summary and conclusions

We have experimentally validated a mathematical model that can predict how the Darcian permeability of welding deposits evolves with time, including the effect of syn-welding volatile diffusion and non-isothermal cooling paths, which are typical of natural scenarios. We use our new experimental and simulation data for the permeability of welded obsidian particles collected in situ to validate our model across a range of temperatures. We take an additional validation step by comparing the component parts of our model to natural and experimental data for other compositions, temperatures, particle sizes, and conditions, and conclude that our model is generally applicable for small systems. Across all isothermal conditions, we identify three dynamic regimes captured by the capillary Peclet number, (1) at $P_c \ll 1$, deposits degas rapidly and weld to preserve equilibrium low H₂O concentration, (2) at $P_c \gg 1$, deposits weld rapidly to preserve high near-initial disequilibrium H₂O concentration, and (3) at intermediate $P_c = 1$, welding produces partially degassed deposits. In cooling systems, we find the additional effect that deposits may cool before they weld or degas, preserving permeable systems that did not completely seal.

We show that, for all reasonable conditions, tuffisites injected in colder country rock are likely to quench before welding completely, preserving permeable pathways. By contrast, we show that ignimbrites straddle the critical cooling rates below which welding can complete and above which it cannot, suggesting that the details of the emplacement temperature, cooling rates, particle radii, and initial H₂O concentration are critical controls on whether an ignimbrite will be welded or not. This same approach shows that if ignimbrites are welded, then whether or not they are fully degassed depends primarily on the particle size. We anticipate that the model framework that we develop will find application in studies of the dynamics of outgassing, and implications for eruptive processes.

## Author contributions

Wadsworth, Vasseur and Llewellin derived the mathematical model and led the data analysis. Vasseur and Wadsworth developed the online resource VolcWeld, which is used and updated herein (Vasseur and Wadsworth, 2019). Brown, Tuffen, and Gardner provided context and images from welded deposits. Schauroth and Dobson led the Avizo™ segmentation work. Wadsworth, von Aulock, Dobson, Schauroth, and
Vasseur assisted at the beamline. Lavallée, Dingwell, and Hess provided laboratory, logistical and financial support. Kendrick, Heap, and Kushnir, contributed additional data and context for application of the model. Marone is the beamline scientist supporting this work. All authors contributed to the manuscript.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

We thank Holly Unwin for useful input in the field in the summer of 2019, and Mathieu Colombier, Sebastian Wiesmaier, and Julie Fife for beamtime support. We acknowledge the Paul Scherrer Institut, Villigen, Switzerland for provision of synchrotron radiation beamtime at the TOMCAT beamline of the SLS. Funding was provided by the European Research Council (Advanced Grant EAVESDROP 834225; Starting Grant SLIM 306488), the Leverhulme Trust (ECF-2016-325 and RF-2019-526 \4), the Natural Environment Research Council grants (NE/N002954/1, NE/M018687/1, and NE/M018687/2), and the National Science Foundation (NSF EAR-1852449). All raw data associated with the experiments presented herein (c.f. Figs. 4-6) and required to reproduce our analysis is available in the Supplementary Information. We are grateful to Kelly Russell for editorial oversight and to two reviewers for careful comments on our work which led to improvements.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jvolgeores.2020.107118.

## References

* Almond, D.C., 1971. Ignimbrite vents in the Sabaloka cauldron, Sudan. Geol. Mag. 108, 159-176. https://doi.org/10.1017/S0016756800051190.

* Andrews, G.D.M., Branney, M.J., 2011. Emplacement and rheomorphic deformation of a large, lava-like rhyolitic ignimbrite: Grey's Landing, southern Idaho. Geol. Soc. Am. Bull. 123, 725-743. https://doi.org/10.1130/b30167.1.

* Bagdassarov, N., Dingwell, D., 1994. Thermal properties of vesicular rhyolite. J. Volcanol. Geotherm. Res. 60, 179-191.

* Berlo, K., Tuffen, H., Smith, V.C., et al., 2013. Element variations in rhyolitic magma resulting from gas transport. Geochim. Cosmochim. Acta 121, 436-451.

* Blair, S.C., Berge, P.A., Berryman, J.G., 1993. Two-Point Correlation Functions to Characterize Microgeometry and Estimate Permeabilities of Synthetic and Natural Sandstones. Lawrence Livermore National Lab, CA (United States).

* Branney, M.J., Kokelaar, B.P., 2002. Pyroclastic Density Currents and the Sedimentation of Ignimbrites. The Geological Society of London, London.

* Branney, M.J., Kokelaar, P., 1992. A reappraisal of ignimbrite emplacement: progressive aggradation and changes from particulate to non-particulate flow during emplacement of high-grade ignimbrite. Bull. Volcanol. 54, 504-520. https://doi.org/10.1007/BF00301396.

* Branney, M.J., Kokelaar, BpP, McConnell, B., 1992. The Bad Step Tuff: a lava-like rheomorphic ignimbrite in a calc-alkaline piecemeal caldera, English Lake District. Bull. Volcanol. 54, 187-199. https://doi.org/10.1007/BF00278388.

* Casas, A.S., Wadsworth, F.B., Ayris, P.M., et al., 2019. SO2 scrubbing during percolation through rhyolitic volcanic domes. Geochim. Cosmochim. Acta https://doi.org/10.1016/j.gca.2019.04.013.

* Castro, J.M., Cordonnier, B., Tuffen, H., et al., 2012. The role of melt-fracture degassing in defusing explosive rhyolite eruptions at volcán Chaitén. Earth Planet. Sci. Lett. 333, 63-69. https://doi.org/10.1016/j.epsl.2012.04.024.

* Castro, J.M., Bindeman, I.N., Tuffen, H., Ian Schipper, C., 2014. Explosive origin of silicic lava: textural and δD-H2O evidence for pyroclastic degassing during rhyolite effusion. Earth Planet. Sci. Lett. 405, 52-61. https://doi.org/10.1016/j.epsl.2014.08.012.

* Cimarelli, C., Yilmaz, T., Colombier, M., et al., 2017. Micro-and nano-CT textural analysis of an experimental volcanic fulgurite. EGU. Gen. Assem. Conf. Abstr. 19, 17982.

* Colombier, M., Wadsworth, F.F.B., Gurioli, L., et al., 2017. The evolution of pore connectivity in volcanic rocks. Earth Planet. Sci. Lett. 462. https://doi.org/10.1016/j.epsl.2017.01.011.

* Crank, J., 1975. The Mathematics of Diffusion. Second Edition. Oxford University Press, New York pp 411.

* Eichheimer, P., Thielmann, M., Fujita, W., et al., 2020. Combined numerical and experimental study of microstructure and permeability in porous granular media. Solid Earth 11, 1079-1095. https://doi.org/10.5194/se-11-1079-2020.

* Elam, W.T., Kerstein, A.R., Rehr, J.J., 1984. Critical properties of the void percolation problem for spheres. Phys. Rev. Lett. 52, 1516.

* Farquharson, J.L., Wadsworth, F.B., Heap, M.J., Baud, P., 2017. Time-dependent permeability evolution in compacting volcanic fracture systems and implications for gas overpressure. J. Volcanol. Geotherm. Res. 339, 81-97. https://doi.org/10.1016/j.jvolgeores.2017.04.025.

* Feng, S., Halperin, B.B.I., Sen, P.P.N., 1987. Transport properties of continuum systems near the percolation threshold. Phys. Rev. B 35, 197.

* Fife, J.L., Rappaz, M., Pistone, M., et al., 2012. Development of a laser-based heating system for in situ synchrotron-based X-ray tomographic microscopy. J. Synchrotron Radiat. 19, 352-358.

* Gardner, J.E., Llewellin, E.W., Watkins, J.M., Befus, K.S., 2017. Formation of obsidian pyroclasts by sintering of ash particles in the volcanic conduit. Earth Planet. Sci. Lett. 459, 252-263. https://doi.org/10.1016/j.EPSL.2016.11.037.

* Gardner, J.E., Wadsworth, F.B., Llewellin, E.W., et al., 2018. Experimental sintering of ash at conduit conditions and implications for the longevity of tuffisites. Bull. Volcanol. 80, 23. https://doi.org/10.1007/s00445-018-1202-8.

* Gardner, J.E., Wadsworth, F.B., Llewellin, E.W., et al., 2019. Experimental constraints on the textures and origin of obsidian pyroclasts. Bull. Volcanol. 81, 22. https://doi.org/10.1007/s00445-019-1283-z.

* Giehl, C., Brooker, R., Marxer, H., Nowak, M., 2016. An experimental simulation of volcanic ash deposition in gas turbines and implications for jet engine safety. Chem. Geol. 461, 160-170.

* Gottsmann, J., Dingwell, D.B., 2001. Cooling dynamics of spatter-fed phonolite obsidian flows on Tenerife, Canary Islands. J. Volcanol. Geotherm. Res. 105, 323-342. https://doi.org/10.1016/S0377-0273(00)00262-6.

* Gottsmann, J., Giordano, D., Dingwell, D.B., 2002. Predicting shear viscosity during volcanic processes at the glass transition: a calorimetric calibration. Earth Planet. Sci. Lett. 198, 417-427.

* Harris, P.D., Branney, M.J., Storey, M., 2011. Large eruption-triggered ocean-island landslide at Tenerife: onshore record and long-term effects on hazardous pyroclastic dispersal. Geology 39, 951-954. https://doi.org/10.1130/G31994.1.

* Heap, M.J., Wadsworth, F.B., 2016. Closing an open system: pore pressure changes in permeable edifice rock at high strain rates. J. Volcanol. Geotherm. Res. 315, 40-50. https://doi.org/10.1016/j.jvolgeores.2016.02.011.

* Heap, M.J., Kolzenburg, S., Russell, J.K., et al., 2014. Conditions and timescales for welding block-and-ash flow deposits. J. Volcanol. Geotherm. Res. 289, 202-209. https://doi.org/10.1016/j.jvolgeores.2014.11.010.

* Heap, M.J., Farquharson, J.L., Wadsworth, F.B., et al., 2015. Timescales for permeability reduction and strength recovery in densifying magma. Earth Planet. Sci. Lett. 429, 223-233. https://doi.org/10.1016/j.epsl.2015.07.053.

* Heap, M.J., Tuffen, H., Wadsworth, F.B., et al., 2019. The permeability evolution of tuffisites and implications for outgassing through dense rhyolitic magma. J. Geophys. Res. Solid Earth 124, 8281-8299. https://doi.org/10.1029/2018JB017035.

* Heap, M.J., Gilg, H.A., Byrne, P.K., et al., 2020. Petrophysical properties, mechanical behavior, and failure modes of impact melt-bearing breccia (suevite) from the Ries impact crater (Germany). Icarus 349, 113873. https://doi.org/10.1016/j.icarus.2020.113873.

* Heiken, G., Wohletz, K., Eichelberger, J., 1988. Fracture fillings and intrusive pyroclasts, Inyo Domes, California. J. Geophys. Res. 93, 4335-4350. https://doi.org/10.1029/JB093iB05p04335.

* Hess, K.U., Dingwell, D.B., 1996. Viscosities of hydrous leucogranitic melts: a non-Arrhenian model. Am. Mineral. 81, 1297-1300.

* Houghton, B., Carey, R.J., 2015. Pyroclastic fall deposits. The Encyclopedia of Volcanoes, pp. 599-616.

* Kendrick, J.E., Lavallée, Y., Varley, N.R., et al., 2016. Blowing off steam: tuffisite formation as a regulator for lava dome eruptions. Front. Earth Sci. 4, 41. https://doi.org/10.3389/feart.2016.00041.

* Kennedy, B.M., Wadsworth, F.B., Vasseur, J., et al., 2016. Surface tension driven processes densify and retain permeability in magma and lava. Earth Planet. Sci. Lett. 433, 116-124. https://doi.org/10.1016/j.epsl.2015.10.031.

* Kertesz, J., 1981. Percolation of holes between overlapping spheres: Monte Carlo calculation of the critical volume fraction. J. Phys. Lett. 42, 393-395.

* Kolzenburg, S., Heap, M.J., Lavallée, Y., et al., 2012. Strength and permeability recovery of tuffisite-bearing andesite. Solid Earth 3, 191-198. https://doi.org/10.5194/se-3-191-2012.

* Kolzenburg, S., Ryan, A.G., Russell, J.K., 2019. Permeability evolution during non-isothermal compaction in volcanic conduits and tuffisite veins: implications for pressure monitoring of volcanic edifices. Earth Planet. Sci. Lett. 527, 115783. https://doi.org/10.1016/j.epsl.2019.115783.

* Lange, R.A., Carmichael, I.S.E., 1987. Densities of Na 2O-K 20-MgO-MgO-FeO-Fe 20 3-Al 30 3-TiO 2-SiO 2 liquids: new measurements and derived partial molar properties. Geochim. Cosmochim. Acta 51, 2931-2946.

* Lavallée, Y., Wadsworth, F.B., Vasseur, J.J., et al., 2015. Eruption and emplacement timescales of ignimbrite super-eruptions from thermo-kinetics of glass shards. Front. Earth Sci. 3, 2. https://doi.org/10.3389/feart.2015.00002.

* Lewiner, T., Lopes, H., Vieira, A.A.W., Tavares, G., 2003. Efficient implementation of marching cubes' cases with topological guarantees. J. Gr. Tools 8, 1-15.

* Liu, Y., Zhang, Y., Behrens, H., 2005. Solubility of H < sub>2</sub> 2<sub>sub</sub> O in rhyolitic melts at low pressures and a new empirical model for mixed H < sub>2</sub> 2<sub>sub</sub> O < sub>CO</sub> < sub>sub</sub> 2<sub>sub</sub> O<sub>2</sub> solubility in rhyolitic melts. J. Volcanol. Geotherm. Res. 143, 219-235.
Llewellin, E.W., 2010a. LBflow: an extensible lattice Boltzmann framework for the simulation of geophysical flows. Part II: usage and validation. Comput. Geosci. 36, 123-132. https://doi.org/10.1016/j.cageo.2009.08.003.

Llewellin, E.W., 2010b. LBflow: an extensible lattice Boltzmann framework for the simulation of geophysical flows. Part I: theory and implementation. Comput. Geosci. 36, 115-122. https://doi.org/10.1016/j.cageo.2009.08.004.

Lorensen, W., Cline, H., 1987. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. ACM siggraph Comput Graph.

Mackenzie, J.K., Shuttleworth, R., 1949. A phenomenological theory of sintering. Proc. Phys. Soc. Sect. B. 62, 833.

Martys, N.S., Torquato, S., Bentz, D.P., 1994. Universal scaling of fluid permeability for sphere packings. Phys. Rev. E 50, 403.

Michaut, C., Bercovici, D., Sparks, R.S.J., 2009. Ascent and compaction of gas rich magma and the effects of hysteretic permeability. Earth Planet. Sci. Lett. 282, 258-267.

Mueller, S., Melnik, O., Spieler, O., et al., 2005. Permeability and degassing of dome lavas undergoing rapid decompression: an experimental determination. Bull. Volcanol. 67, 526-538.

Mueller, S.P., Helo, C., Keller, F., et al., 2018. First experimental observations on melting and chemical modification of volcanic ash during lightning interaction. Sci. Rep. 8, 1389. https://doi.org/10.1038/s41598-018-19608-3.

Okumura, S., Sasaki, O., 2014. Permeability reduction of fractured rhyolite in volcanic conduits and its control on eruption cyclicity. Geology 42, 843-846.

Owen, J., Shea, T., Tuffen, H., 2019. Basalt, unveiling fluid-filled fractures, inducing sediment intra-void transport, ephemerally: examples from Katla 1918. J. Volcanol. Geotherm. Res. 369, 121-144. https://doi.org/10.1016/j.jvolgeores.2018.11.002.

Paisley, R., Berlo, K., Whattam, J., et al., 2019. Degassing-induced chemical heterogeneity at the 2011-2012 Cordón Caulle eruption. Volcanica 2, 211-237. https://doi.org/10.30909/vol.02.02.211237.

Parikh, N.M., 1958. Effect of atmosphere on surface tension of glass. J. Am. Ceram. Soc. 41, 18-22.

Prado, M., Dutra Zanotto, E., Müller, R., 2001. Model for sintering polydispersed glass particles. J. Non-Cryst. Solids 279, 169-178.

Quane, S.L., Russell, J.K., 2005a. Welding: insights from high-temperature analogue experiments. J. Volcanol. Geotherm. Res. 142, 67-87. https://doi.org/10.1016/j.jvolgeores.2004.10.014.

Quane, S.L., Russell, J.K., 2005b. Ranking welding intensity in pyroclastic deposits. Bull. Volcanol. 67, 129-143.

Quane, S.L., Russell, J.K., Friedlander, E.A., 2009. Time scales of compaction in volcanic systems. Geology 37, 471-474.

Rahaman, M.N., De Jonghe, L.C., 1990. Sintering of spherical glass powder under a uniaxial stress. J. Am. Ceram. Soc. 73, 707-712. https://doi.org/10.1111/j.1151-2916.1990.tb06576.x.

Riehle, J.R., 1973. Calculated compaction profiles of rhyolitic ash-flow tufts. Bull. Geol. Soc. Am. 84, 2193-2216. https://doi.org/10.1130/0016-7606(1973)84 <2193:CCPOR>A.2.0.CO;2

Riehle, J.R., Miller, T.F., Bailey, R.A., 1995. Cooling, degassing and compaction of rhyolitic ash flow tufts: a computational model. Bull. Volcanol. 57, 319-336. https://doi.org/10.1007/BF00301291.

Russell, J.K., Quane, S.L., 2005. Rheology of welding: inversion of field constraints. J. Volcanol. Geotherm. Res. 142, 173-191. https://doi.org/10.1016/j.jvolgeores.2004.10.017.

Ryan, A.G., Friedlander, E.A., Russell, J.K., et al., 2018. Hot pressing in conduit faults during lava dome extrusion: insights from Mount St. Helens 2004-2008. Earth Planet. Sci. Lett. 482, 171-180. https://doi.org/10.1016/j.epsl.2017.11.010.

Ryan, A.G., Heap, M.J., Russell, J.K., et al., 2020a. Cyclic shear zone cataclasis and sintering during lava dome extrusion: insights from Chaos Crags, Lassen Volcanic Center (USA). J. Volcanol. Geotherm. Res. 401, 106935. https://doi.org/10.1016/j.jvolgeores.2020.106935.

Ryan, A.G., Russell, J.K., Heap, M.J., et al., 2020b. Timescales of porosity and permeability loss by solid-state sintering. Earth Planet. Sci. Lett. 549, 116533. https://doi.org/10.1016/j.epsl.2020.116533.

Sánchez, M.C., Sarrionandia, F., Arostegui, J., et al., 2012. The transition of spatter to lava-like body in lava fountain deposits: features and examples from the Cabezo Segura volcano (Calatrava, Spain). J. Volcanol. Geotherm. Res. 227, 1-14.

Saubin, E., Tuffen, H., Gurioli, L., et al., 2016. Conduit dynamics in transitional rhyolitic activity recorded by tuffisite vein textures from the 2008-2009 Chaitén Eruption. Front. Earth Sci. 4. https://doi.org/10.3389/feart.2016.00059.

Schipper, C.I., Castro, J.M., Tuffen, H., et al., 2015. Cristobalite in the 2011-2012 Cordón Caulle eruption (Chile). Bull. Volcanol. 77. https://doi.org/10.1007/s00445-015-0925-z.

Song, W., Hess, K.K.-U., Damby, D.E.D., et al., 2014. Fusion characteristics of volcanic ash relevant to aviation hazards. Geophys. Res. Lett. 41, 2326-2333. https://doi.org/10.1002/2013GL059182.

Sparks, R.S.J., Tait, S.R., Yanev, Y., 1999. Dense welding caused by volatile resorption. J. Geol. Soc. Lond. 156, 217-225. https://doi.org/10.1144/gsjgs.156.2.0217.

Stasiuk, M.V., Barclay, J., Carroll, M.R., et al., 1996. Degassing during magma ascent in the Mule Creek vent (USA). Bull. Volcanol. 58, 117-130.

Streck, M.J., Grunder, A.L., 1995. Crystallization and welding variations in a widespread ignimbrite sheet; the Rattlesnake Tuff, eastern Oregon, USA. Bull. Volcanol. 57, 151-169.

Sumner, J.M., Branney, M.J., 2002. The emplacement history of a remarkable heterogeneous, chemically zoned, rheomorphic and locally lava-like ignimbrite: 'TL' on Gran Canaria. J. Volcanol. Geotherm. Res. 115, 109-138. https://doi.org/10.1016/S0377-0273(01)00311-0.

Sumner, J.M., Blake, S., Matela, R.J., Wolff, J.A., 2005. Spatter. J. Volcanol. Geotherm. Res. 142, 49-65.

Torquato, S., 2013. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer Science & Business Media.

Tuffen, H., Dingwell, D., 2005. Fault textures in volcanic conduits: evidence for seismic trigger mechanisms during silicic eruptions. Bull. Volcanol. 67, 370-387. https://doi.org/10.1007/s00445-004-0383-5.

Tuffen, H., Dingwell, D.B., Pinkerton, H., 2003. Repeated fracture and healing of silicic magma generate flow banding and earthquakes? Geology 31, 1089-1092. https://doi.org/10.1130/g19777.1.

Vasseur, J., Wadsworth, F., 2019. Volcanic Welding Model: VolcWeld. https://vhub.org/resources/4568.

Vasseur, J., Wadsworth, F.B., 2017. Sphere models for pore geometry and fluid permeability in heterogeneous magmas. Bull. Volcanol. 79, 77. https://doi.org/10.1007/s00445-017-1165-1.

Vasseur, J., Wadsworth, F.B., Lavallée, Y., et al., 2013. Volcanic sintering: timescales of viscous densification and strength recovery. Geophys. Res. Lett. 40, 5658-5664. https://doi.org/10.1002/2013GL058105.

Vasseur, J., Wadsworth, F.B., Dingwell, D.B., 2020. Permeability of polydisperse magma foam. Geology. https://doi.org/10.1130/G47094.1.

Von Aulock, F.W., Kennedy, B.M., Maksimenko, A., et al., 2017. Outgassing from open and closed magma foams. Front. Earth Sci. 5. https://doi.org/10.3389/feart.2017.00046.

Wadsworth, F.B., Vasseur, J., Aulock, F.W., et al., 2014. Nonisotopic viscous sintering of volcanic ash. J. Geophys. Res. Solid Earth 119, 8792-8804. https://doi.org/10.1002/2014JB011453.

Wadsworth, F.B., Vasseur, J., Llewellin, E.W., et al., 2016a. Sintering of viscous droplets under surface tension. Proc. R. Soc. A. Math. Phys. Eng. Sci. 472, 20150780. https://doi.org/10.1098/rspa.2015.0780.

Wadsworth, F.B., Vasseur, J., Scheu, B., et al., 2016b. Universal scaling of fluid permeability during volcanic welding and sediment diagenesis. Geology 44, 219-222. https://doi.org/10.1130/G37559.1.

Wadsworth, F.B., Vasseur, J., Llewellin, E.W., et al., 2017a. Size limits for rounding of volcanic ash particles heated by lightning. J. Geophys. Res. Solid Earth 122, 1977-1989. https://doi.org/10.1002/2016JB013864.

Wadsworth, F.B., Vasseur, J., Llewellin, E.W., et al., 2017b. Topological inversions in coalescing granular media control fluid-flow regimes. Phys. Rev. E 96, 033113. https://doi.org/10.1103/PhysRevE.96.033113.

Wadsworth, F.B., Vasseur, J., Llewellin, E.W., Dingwell, D.B., 2017c. Sintering of polydisperse viscous droplets. Phys. Rev. E 95. https://doi.org/10.1103/PhysRevE.95.033114.

Wadsworth, F.B., Vasseur, J., Schauroth, J., et al., 2019. A general model for welding of ash particles in volcanic systems validated using in situ X-ray tomography. Earth Planet. Sci. Lett. 525, 115726. https://doi.org/10.1016/j.epsl.2019.115726.

Wadsworth, F.B., Llewellin, E.W., Vasseur, J., et al., 2020a. Explosive-effusive volcanic eruption transitions caused by sintering. Sci. Adv. 6, 7940-7963. https://doi.org/10.1126/sciadv.aba7940.

Wadsworth, F.B., Vossen, C.E.J., Schmid, D., et al., 2020b. Determination of permeability using a classic Darcy water column. Am. J. Phys. 88, 20-24. https://doi.org/10.1119/10.000296.

Wright, H.M., Cashman, K.V., 2014. Compaction and gas loss in welded pyroclastic deposits as revealed by porosity, permeability, and electrical conductivity measurements of the Shevlin Park Tuff. Geol. Soc. Am. Bull. 126, 234-247. https://doi.org/10.1130/b30668.1.

Zhang, Y., Ni, H., 2010. Diffusion of H, C, and O components in silicate melts. Rev. Miner. Geochem. 72, 171-225.