A model for permeability evolution during volcanic welding

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 2 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 2 O content. A user-friendly implementation of the model is provided. 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 , deposits degas rapidly and weld to preserve equilibrium low H 2 O concentration, (2) at , deposits weld rapidly to preserve high near-initial disequilibrium H 2 O concentration, and (3) at intermediate 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, 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 2 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.

Welding is important because it leads to substantial changes in the physical properties of the deposit, which transitions progressively from particulate to coherent , from mechanically weak to mechanically strong (Kolzenburg et al. 2012;Vasseur et al. 2013), and from permeable to impermeable 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. 2014Wadsworth et al. , 2019Kolzenburg et al. 2019;Heap et al. 2019), which may be complex.

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, volatilesupersaturated 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.

Evolution of porosity during welding
We adopt the model of Wadsworth et al. (2016aWadsworth et al. ( , 2017bWadsworth et al. ( , 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' . 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) with time as a function of material properties: where is the initial porosity, is the initial radius of the 'bubbles', is the melt viscosity, is the surface tension, and is the difference between the squeezing isotropic pressure acting on the particles and the gas pressure in the bubbles , such that . In nature, 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 : Eq. 2 where, ̅ , ̅ , and ̅ are, respectively, dimensionless time, pressure difference, and porosity. Time is normalized by the capillary relaxation timescale , pressure difference is normalized by a reference capillary pressure scale , and porosity is normalized by the initial porosity as follows Melt viscosity is a strong function of temperature , which may vary over time, and dissolved H 2 O concentration , 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 2 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 〈 〉. The spatially-averaged viscosity ( 〈 〉) is then used in the integral within the definition of ̅ , which runs from the time at which the welding starts to the time-step of interest (Eq. 3).
The equations for the rate of porosity changein dimensional form (Eq. 1) or dimensionless form (Eq. 2)can be integrated numerically to derive porosity change as a function of time ( ) for welding particles. Wadsworth et al., (2019) describe a suitable numerical scheme and provide a userfriendly implementation for download. The full model accounts for particle polydispersivity, nonisothermal 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., (2017b) and used here.
2.2 Evolution of permeability during welding Section 2.1 presents a model for the evolution of porosity with time ( ) in a welding deposit. We can extend this to derive a model for the evolution of Darcian permeability with time ( ) by finding a constitutive relationship for ( ). For sintering systems of particles, the most widely validated model is based on percolation theory (Martys et al. 1994;Wadsworth et al. 2016b): where is the specific surface area, is the percolation threshold, and is a percolation exponent. For volcanic applications, this model has been validated for sintering glass spheres (Wadsworth et al. 2017b(Wadsworth et al. , 2020bEichheimer 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, b). Whereas and are typically measured quantities, is either measured or requires estimation, and and 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 and , 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 (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 (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 directly. However, to allow Eq. 4 to be applied in situations where is not measured, we also test a constitutive model for ( ). 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 (Torquato 2013) where 〈 〉 is the th moment of the distribution of pore sizes. Wadsworth et al. (2017c) provide a method for predicting the distribution of pore sizes, and therefore 〈 〉, in sintering systems that have a distribution of initial particle sizes.
The above framework provides a method by which ( ) and ( ) 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.

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 〈 〉 m, were polydisperse in grain size, and were crystal-free. The size distribution of particles is reported in Fig. 2 along with 〈 〉 for .
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 ( ). Four runs were performed, in which the packed free-standing cylinders were heated at to temperatures and 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 voxel edge lengths. The initial H 2 O concentration was 0.14 wt.% .
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 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 (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 . Using a range of suitable segmentations, we determine the uncertainty on the porosity to be approximately 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 J o u r n a l P r e -p r o o f ( Fig. 3a-c). The pore space is initially fully connected ( Fig. 3a) but develops domains that are isolated (Fig. 3b) as welding 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 dataresults 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, b). 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 and density ) under a uniform fluid pressure gradient of . The fluid velocity is read out at every node and the average fluid velocity 〈 〉 is computed. We run the simulation until 〈 〉 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:

〈 〉
Eq. 6 The specific surface area 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 poreparticle walls in the 3D datasets. Normalizing the area of the surface by the volume of the sample or domain gives (in units of m -1 ). Using 3D datasets segmented using different segmentation threshold choices, we estimate the uncertainty on the permeability to be based on the convergence criterion used in LBflow. The uncertainty on our surface area determination is approximately . This uncertainty is determined by computing the standard error on repeat measurements of surface area on different sub-volumes selected in domain.

Numerical modelling
We solve Eq. 2 numerically, following Wadsworth et al. (2019), to obtain ̅ ( ̅ ), which we dimensionalize for each experiment via the relationships in Eq. 3, and the material parameters , , 〈 〉, and . Calculating 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 〈 〉. We use constitutive models appropriate for rhyolitic liquids: viscosity ( ) from Hess and Dingwell, (1996); H 2 O diffusivity ( ) from Zhang and Ni, (2010); and the equilibrium H 2 O solubility at the particle boundary ( ) from Liu et al., (2005), where is the partial pressure of H 2 O in the pore space. For our experiments MPa, and we use to account for standard laboratory humidity (Von Aulock et al. 2017). We use N. m -1 for nominally dry silicate melts (Parikh 1958). The average initial pore radius 〈 〉 for the experiments is determined using Wadsworth et al. (2017b), and the starting porosity is taken directly from segmented tomography data.

Results and analysis
Our results and analysis comprise three parts. First, we validate the ( ) model (Eq. 5; Section 2.2) against the values determined from the tomographic data. Second, we compare the ( ) model (Eq. 4; Section 2.2) with ( ) data from the LBflow simulations of the in situ data, for which all quantities are measured directly. Third, we validate the ( ) model (i.e. the combined ( ) and ( ) models (Eqs 1-3; Section 2.1)) against ( ) data from the LBflow simulations.

The specific surface area during welding ( )
The specific surface area measured during welding in the evolving in situ datasets varies nonlinearly with . In Fig. 4 we show the measured 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. 2017d), which uses the particle size distribution (Fig. 2) and the initial porosity as input values. In the case of the data presented here, the pore size distribution moments are 〈 〉 m, 〈 〉 , and 〈 〉 . With no fitting parameters, the data and the model agree with a coefficient of determination of .
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 taken at the initial packing porositytermed can be used to estimate the excess surface. We find that 〈 〉 , 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 , but a surface area that is 20% higher. Such a spheroid has major semi-axis and minor semi-axes . The minimum radius of curvature of the spheroid is ⁄ (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 is given by (Wadsworth et al. 2017a). For 〈 〉 m, the rounding time for our experiments at K is seconds. In our data, the porosity at which is met is approximately at a porosity of , which is consistent with the porosity at which 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 ( ) to account for angularity (e.g. by adjusting the moments of the pore size distribution by an empirical factor).

The relationship between permeability and porosity ( )
The permeability of the samples determined through fluid-flow simulations depends strongly on porosity ( Fig. 5a note the logarithmic scale for ). Permeability varies by two orders of magnitude over the porosity interval and, at the lowest porosities, simulations show that the samples are impermeable (i.e. ). Fig. 5a also plots the theoretical ( ) relationship in Eq. 4 using the model for ( ) from Eq. 5 (see Section 4.1), , and . This shows a generally good match to the data. Improvements in the fit of ( ) (Eqs 4 & 5) could be made if and were treated as fitting parameters. In the inset to Fig. 5a, we perform a minimization exercise in which we allow and 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 and are accounted for. The coefficient of determination is then the ratio of this sum of square residuals for a given and 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 and 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 and . We note that this 2% variation in the fit quality encompasses the theoretically-grounded values and , derived for monodisperse packs of spherical particles. Based on this analysis, we conclude that: 1) the quality of the model fit to the ( ) data is not sensitive to small variations (±2%) in the best-fit parameters; and 2) that the theoreticallygrounded values of the fit parameters remain valid for packs of natural volcanic particles, despite violation of the monodisperse and spherical assumptions. We therefore adopt and as empirically-validated, theoretically-grounded constants.
Due to the in situ nature of our experiments, we extract the connected porosity directly. If we then look at the connectivity of the pore space (as defined in Colombier et al. 2017) as a function of during welding, we would predict that this metric would fall to zero as during welding. The value is consistent with the observed drop in connected porosity (Fig. 5b).
The data indicate that there is no strong dependence of the ( ) 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 with respect to as , in Figure 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).

Time-evolution of permeability ( )
With a validated constitutive law for ( ) (Eqs 4-5), we can combine our model for ( ) with the vented bubble model for ( ) (Eqs 1-3), giving a method for calculating the time-dependence of permeability during welding ( ). Fluid flow simulations through the segmented 3D datasets from our in situ experiments yield ( ) 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 ( ) (Eq. 5; Fig. 3) is used to solve the model for ( ) (Eq. 4) using and (see Section 4.2; Fig. 5). This solution for ( ) is then mapped to each by first solving the dimensionless universal form of the vented bubble J o u r n a l P r e -p r o o f 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 ( ) is presented in detail elsewhere for these same experiments  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 2 O concentration in the glass throughout (Fig. 6). Neglecting diffusive degassing of the particles provides a very poor description of the data and confirms that, for our experimental conditions, accurate prediction of permeability during welding requires that degassing is accounted for.

Discussion
The results in Section 4 constitute an experimental validation of the workflow for modelling ( ) via separate models for ( ), ( ), and ( ). 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 ̅ ), 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.

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 ( ), ( ), and ( ) 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 ( ) model (Eq. 1) is already presented in dimensionless form ̅ ( ̅ ) (Eqs. 2 and 3) and has been validated extensively over a range of conditions (Wadsworth et al. 2016a(Wadsworth et al. , 2017b(Wadsworth et al. , d, 2019. We non-dimensionalize the specific surface area by the moments of the pore size distribution (see Section 4.1), giving ̅ 〈 〉 〈 〉 for the polydisperse case, and 〈 〉 for the monodisperse case. The model (Eq. 5) reduces to the dimensionless form ̅ ( ) ( ), which is shown in Figure 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 〈 〉 〈 〉 or 〈 〉 are computed following the same technique as used here (Torquato 2013;Wadsworth et al. 2017d). Data are presented in Table 1. Despite the very different pore size and particle size distributions in these studies, the collapse to the dimensionless ̅ ( ) model is generally reliable. We note that some samples from Wadsworth et al. (2017b) have ̅ 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 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 , 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 .

Journal Pre-proof
We non-dimensionalise the permeability by a reference permeability giving ̅ . The model (Eq. 4) reduces to the dimensionless form ̅ ( ) , which is shown in Figure 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., (2017a), 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 andesite particles (Kendrick et al. 2016); and natural welded systems for which is constrained (Heap et al. 2015(Heap et al. , 2019Ryan et al. 2020a). Where the original study does not give a value for , 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 .
Finally, in Fig. 7c, we test our ( ) 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 ̅ , termed ̅ , from which each specific experiment starts the process is different. For any given ̅ , there will be a different model result in ̅ ( ). Therefore, to provide a common basis that allows us to compare data and model results across all conditions, we use ̅ ̅ , which we note is equivalent to (i.e., the 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 ̅ 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 2 O, which is accounted for in the sintering obsidian particles analyzed herein . When cast in this space (Fig. 7c), the model given in Eq. 4 becomes ̅ ̅ [( ) ( )] and, following the same principles as is shown dimensionally in Fig.  6, the time dependence for ̅ ( ̅ ) is used together with ( ) to predict ̅ ( ̅ ). For the sintering of sodalime-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 ( ), it does not strongly influence the dynamics of permeability decrease. Similarly, this suggests 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).

Successes and limitations of our model and validation
The in situ nature of the data allows us to determine , , and , independently, as functions of measured experimental time ; hence to validate, independently, our models for ( ) (Fig. 4), ( ) (Fig. 5), and ( ) (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 ̅ (Eq. 3). Wadsworth et al. (2019) showed that Eqs 1-3 are valid descriptions of ( ) during welding, even when confining pressure exceeds the Laplace surface pressurei.e. when ̅ and validated up to ̅ . However, this does not necessarily confirm that the combination of ( ) and ( ), to determine ( ), would also be valid at higher confining pressures. It could be that the so-called J o u r n a l P r e -p r o o f pressure sintering regime at ̅ , induces particle-particle flattening sufficient to change the microstructure such that the analogy with static percolation theory developed for overlapping spheres, which underpins our ( ) model, is no longer applicable. For instance, as ̅ 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 ̅ , and validation of ( ) at higher ̅ should be a topic for future investigation. Nevertheless, the low ̅ regime is relevant to a range of volcanic scenarios (Gardner et al. 2018(Gardner et al. , 2019Wadsworth 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 timedependent 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 ( ) pathway taken during welding.
As a coarse test of the impact of the complexities discussed above, we show that our ( ) 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. (2020) 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 wt.% H 2 O, while the far-field hosts an equilibrium value of wt.% H 2 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 ̅ for this tuffisite. Using a solubility model (Liu et al. 2005) we find that the wall value relates to MPa and that the far-field value relates to MPa, resulting in MPa. Given the constraints on from Heap et al. (2019), we therefore estimate that ̅ for that tuffisite. Despite the high ̅ 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 from our model at any porosity (Fig. 7b). We propose that such a collapse in ( ) 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.

Applications and regimes in natural welding
J o u r n a l P r e -p r o o f

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 , initial monodisperse particle size , isothermal temperature , 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 , which represents an insulated magmatic temperature (cooling is discussed later), and an initial porosity , which is typical of loose randomly packed angular particles (e.g. Fig. 5). We vary and over a wide range representative of silicic systems: wt.%, and m, where √ ( ) is the radius above which gravitational effects dominate the particle deformation over surface tension effects, is the rhyolite particle density (Lange and Carmichael 1987), and is the acceleration due to gravity. For rhyolites, m, such that particles smaller than this value will not deform under their own mass. We assume that welding is occurring at MPa vapour 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-̅ regime (discussed above). At these conditions, the solubility model we used for our experimental validation predicts that the equilibrium water concentration is 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 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 , we find more complex welding pathways can be taken in ( ) space, depending on particle radius. We find three phenomenological regimes; to demonstrate these regimes, we use the example of ( Fig.  8d).
(1) For relatively small particles ( in this example) permeability decay is almost indistinguishable from the case where (Fig. 8a). Inspection of the evolution of the spatially-averaged H 2 O concentration 〈 〉( ) (Fig. 8e) demonstrates that this is because 〈 〉 falls rapidly to , before appreciable welding has occurred: i.e. degassing is efficient on the timescale of welding. Consequently, the particles have equilibrium H 2 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.
(2) For relatively large particles ( in this example) permeability decays much more rapidly than for the case where , and more rapidly even than the smaller particles discussed in (1) above. Inspection of 〈 〉( ) demonstrates that this is because H 2 O concentration remains close to 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 2 O concentration when welding is complete.
(3) For particles of intermediate size ( 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 and when welding is complete.

J o u r n a l P r e -p r o o f
To explore these regimes further, we consider the total time taken for welding to complete for particles of different size and initial H 2 O concentration. We take as the time elapsed when reaches . In Fig. 9 we plot as a function of the starting particle radius, showing the solution for a wide range wt.%, for rhyolitic particles welding at K and MPa. We observe that, for a given , 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 2 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 (Gardner et al. 2018), which compares the timescale for diffusive degassing with the timescale for welding in the low-̅ regime, (Eq. 3). The diffusive timescale is . For convenience, we can define for the initial conditions,

Eq. 7
where and are the properties and at the initial H 2 O concentration prior to diffusion. If , then degassing is rapid compared with welding, and early in the welding process. This manifests as the viscosity moving rapidly from an initial value determined by , toward a final value determined by . However, if , then welding is rapid compared with degassing, and does not substantially change during the welding process. This manifests as a near-constant viscosity at determined by . Then, in the limiting cases of and , we expect the time at which the permeability reaches to zero to be approximated by the welding time such that Eq. 8b and for , 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 〈 〉 and are typically of the same order of magnitude to make the approximation 〈 〉 , rendering easy to compute in these limiting cases. In Fig. 9, we plot 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 in place of 〈 〉. The limit at (Eq. 8b) is not a single curve and depends on (Fig. 9).
The regime boundary between and is not necessarily precisely at To test this, we find the point at which for each , which, by definition, is the point at which , and plot this continuously in Fig. 9. The curve cuts the ( ) curves within their transitional zones (i.e. the inflected portion) in the majority of cases, indicating that is a reasonable approximation of the regime boundary. The transitional zone is more sharply defined for high that for low , and the absolute difference between the low and high solutions is greater high .

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 J o u r n a l P r e -p r o o f 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 . While there is evidence that depends on the cooling rate,-termed , as (Gottsmann et al. 2002), we make the simplifying assumption that is the temperature at which 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 over the interval between an initial temperature and , which is usually a few hundred kelvin, and therefore relatively small compared with the total cooling path from to ambient conditions. During linear cooling, the viscosity of the particles increases and the diffusivity of H 2 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 . At the fastest cooling rates, little welding occurs before is reached, and the initial permeability is preserved. At the slowest cooling rates, welding goes to completion, and permeability drops to zero, by the time 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 , hence no degassing occurs, we can see that the final permeability is a simple, smooth function of , whereas when , the final permeability may be a more complex function of .
There is a critical , which we term , below which the permeability will always fall to zero before (Fig. 10). We take 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 as a function of for the same range of as given for isothermal conditions in Fig. 9. This shows that, for relatively small particles, the critical cooling rate collapses to a limiting ( ) for any ; for relatively large particles, ( ) evolves toward a second limiting solution at relatively higher . Following the analysis in Section 6.1, we find that the limiting solution ( ) at small is the solution for the . Conversely, the high ( ) limit at large is the solution for the . This implies that there are again three regimes for a given (1) At small and , welding completes and the resulting deposit is degassed to equilibrium H 2 O concentration.
(2) At large and , welding completes but the resulting deposit retains disequilibrium H 2 O concentration.
(3) At any and , welding is incomplete, , and 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 Eq. 9b J o u r n a l P r e -p r o o f and as before these regimes are separated by (solid curve on Fig. 11).

Direct volcanic applications
Welding occurs in volcanic eruptions in distinct scenarios: on the walls of open conduits during Plinian eruptions (Gardner et al. 2017(Gardner et al. , 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 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 , 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 2 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, to , 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 viscoelastic 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 ( ) 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 〈 〉 ( ) where is the thermal diffusivity and is the system lengthscale in the dominant direction of heat flow. Using data from Bagdassarov and Dingwell (1994) for rhyolitic melts, we can approximate as independent of temperature over the range K and take . For ignimbrites, we assume that the lower limit on is the depositional boundary layer thickness, estimated to be on the order of (Andrews and Branney 2011) and the upper limit is taken to be m as a reasonable estimate of half-thickness of an entire ignimbrite. For tuffisites, we take the lower limit on to be and the upper limit to be based on available observations (Castro et al. 2012;Saubin et al. 2016;Heap et al. 2019). In all cases we take to be K. This results in 〈 〉 for ignimbrites, consistent with geospeedometry measurements from glassy welded ignimbrites (Gottsmann and Dingwell 2001), and 〈 〉 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. fiammé), the dominant particle size range is Streck and Grunder 1995;Andrews and Branney 2011) and a maximum initial H 2 O concentration estimate is (Sumner J o u r n a l P r e -p r o o f . For tuffisites, the particle size range is similar and (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 as a function of radius of the welding particles (as in Figure 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 magmatic 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 , 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 , 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 2 O concentrations.
Both the ignimbrite and tuffisite scenarios investigated in Fig. 12 are in the low-̅ 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-̅ regime. An increase in ̅ moves the ( ) 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 and the 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 tuffisites is km and for a 10 cm tuffisite is m.
Referring back to Fig. 9, isothermal tuffisites injected in hot magma are likely to weld rapidly ) and to record degassed H 2 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 . 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 bombinformation 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).

Outlook: permeability-limited welding and compaction
In a welding deposit, there is a competition between the timescale over which welding occurs (in the low-̅ regime), and the timescale over which gas escapes from the welding material , 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 , then 〈 〉 ( ) where is the J o u r n a l P r e -p r o o f distance the gas has to travel to escape the welding system. Using we can define a Darcy number as (Wadsworth et al. 2016a) Eq. 10 such that at , outgassing is efficient and welding will complete in the way described by our model. However for , 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 regime, our analysis shows that many volcanic systems weld to completion and permeability drops to zero. Since as (Eq. 10); it's clear that 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 behaviour (Castro et al. 2014). This regime remains under-explored.

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 ( ) and ( ). 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 2 O. In turn, this diffusion solution is used to solve the 'vented bubble model' for welding and output ( ) and ( ). The Excel™ implementation (Supplementary Information), neglects the diffusion of H 2 O, and therefore is suited to conditions where H 2 O movement can be neglected on the timescale of welding ( and ; 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 ( ) and find the mean particle size 〈 〉. Both the VolcWeld program and the Excel™ implementation (Supplementary Information) use Wadsworth et al. (2017c) to convert ( ) into a pore size distribution ( ) and find the mean pore size 〈 〉. 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 〈 〉.
J o u r n a l P r e -p r o o f (2) Measure the initial H 2 O concentration of the particles and determine whether or not this is in equilibrium by calculating the H 2 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 2 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 (e.g. Hess and Dingwell 1996), use , and use 〈 〉 from step (1) to compute a welding timescale.
(4) If the system is not in H 2 O equilibrium (see step (2) For non-isothermal cooling conditions: (5) As with isothermal conditions, the Excel™ implementation can be used for or , where the first sheet in the implementation permits input of a temperature rate If, however, , then the VolcWeld program must be used.

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 , deposits degas rapidly and weld to preserve equilibrium low H 2 O concentration, (2) at , deposits weld rapidly to preserve high near-initial disequilibrium H 2 O concentration, and (3) at intermediate 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, 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 2 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.   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). (gi) 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. and (Feng et al. 1987). The grey band represents the range of model solutions that encompass 2% uncertainty around the best fit and pair. Inset: the goodness of fit (fit in logarithmic space) given by a Deming (see text) of Eq. 4 for different and value pairs, for which the best-fit value is indicated with the dashed line. The region of and value pairs given by the grey band in (a) is enclosed by the white dashed line. (b) The relationship between the pore connectivity and .

Figure 8.
Model results for permeability as a function of time ( ) and the average particle H 2 O concentration as a function of time 〈 〉( ) for different particle radius and initial water concentration : (a-b) for ; (c-d) for and for which is determined at 0.1 MPa. In (b) and (d), the filled circle against which the curves terminate represents the time at which welding completes . All results are using the ( ) relationship for rhyolites (Hess and Dingwell 1996) and at K. c) Schematic cartoon showing the simultaneous evolution of pore-space geometry and H 2 O distribution in the particle pack, as welding and degassing proceed.