Heat Transfer in Pyroclastic Density Current‐Ice Interactions: Insights From Experimental and Numerical Simulations

Stratovolcanoes are common globally, with high‐altitude summit regions that are often glacier‐clad and intersect the seasonal and perennial snow line. During an eruption, interaction between snow/ice and hot, pyroclastic deposits will potentially lead to extensive melt and steam production. This is particularly pertinent when pyroclastic density currents (PDCs) are emplaced onto and propagate over glacierised substrates. Generated melt and steam are incorporated into the flow, which can cause a transformation from a hot, dry granular flow, to a water‐saturated, sediment‐laden flow, termed a lahar. Both PDCs and ice‐melt lahars are highly hazardous due to their high energy during flow and long runout distances. Knowledge of the physics that underpin these interactions and the transformation to ice‐melt lahar is extremely limited, preventing accurate descriptions within hazard models. To physically constrain the thermal interactions we conduct static melting experiments, where a hot granular layer was emplaced onto an ice substrate. The rate of heat transfer through the particle layer, melt and steam generation were quantified. Experiments revealed systematic increases in melt and steam with increasing particle layer thicknesses and temperatures. We also present a one‐dimensional numerical model for heat transfer, calibrated against experimental data, capable of accurately predicting temperature and associated melting. Furthermore, similarity solutions are presented for early‐time melting which are used to benchmark our numerical scheme, and to provide rapid estimates for meltwater flux hydrographs. These data are vital for predicting melt volume and incorporation into PDCs required to facilitate the transformation to and evolution of ice‐melt lahars.

When PDCs propagate over and are emplaced onto snow or ice they mechanically and thermally scour the substrate (Pierson et al., 1990;Thouret et al., 2007;Walder, 2000b).These processes can generate steam and meltwater that are incorporated into the flow, which can fundamentally affect its dynamics and cause transformations in both flow mobility and character (Figure 1).Generation and escape of steam can increase pore fluid pressures within the granular layer and cause fluidization of the flow, enhancing its overall mobility (Roche et al., 2002;Rowley et al., 2014).Incorporation of meltwater can affect the friction and cohesion properties of the bulk particle layer and can transform the flow into an ice-melt lahar if sufficient melt is generated (Ahmed et al., 2012;Branney & Gilbert, 1995;Huggel et al., 2007;Walding et al., 2023).
PDC interactions with frozen substrates are difficult to study in situ due to their unpredictable and hazardous nature and poor preservation potential due to the susceptibility of snow and ice to melt out of deposits, reworking them in the process (Breard, Calder, & Ruth, 2020).Few evidence-based field studies of PDC-ice interactions and subsequent lahar generation exist.For example, highly detailed investigations were conducted following the catastrophic 1985 eruption of Nevado del Ruiz.These studies provided constraints on (a) total ice loss and melt volume (Thouret, 1990), and (b) PDC, tephra fall and lahar events and deposits (Naranjo et al., 1986;Pierson et al., 1990).Kilgour et al. (2010) also provides a detailed study following the 25 September 2007 eruption of Ruapehu, which generated small-volume lahars.These studies provide context for the geophysical scale modeling described in Section 4.2.Experiments and theoretical modeling complement detailed field studies and can offer additional insights into the microphysical interactions between PDCs and frozen substrates.
In order to comprehensively investigate the physics underpinning PDC-ice interactions, the thermal and mechanical effects must be disentangled.Our study builds on previous work where the thermal effects are studied in isolation by considering the scenario where particles are rapidly placed on a horizontal frozen surface, so that there is no mechanical scour.Few theoretical and experimental studies of hot particle-ice interactions exist offering insights into the heat transfer from particle to ice layers.Walder (2000aWalder ( , 2000b) ) developed a theory for pyroclast-snow interactions and thermally-driven slurry formation based on vertical thermal transfer between a porous, hot, monodisperse, and polydisperse particle layer and a snow substrate.Experiments where hot particles were released onto shaved ice revealed a continuum of behaviors.Where no convective bubbling occurred within the granular layer, the particles melted into the snow as a wetting front rose upwards through the particle layer.In other cases rising vapor bubbles caused complete convective overturning of the particle layer by fluidization, thermally scouring, and incorporating the snow, facilitating the transformation from a dry non-cohesive mass of particles into a slurry.The latter regimes were favored by higher initial particle temperatures and smaller grain diameters.Cowlyn (2016) conducted complementary experiments to determine the amount of melting that could be generated by an individual pyroclast.These experiments provided numerical constraints for the rate of melt and steam production to inform predictive models for PDC-ice interaction.In this paper we focus on the thermal interactions between a layer of hot particles and ice.The previous experimental work of Walder (2000b) and Cowlyn (2016) are extended to include the interactions of volcanic and non-volcanic particles with ice across the expected thermal range for PDC-ice interactions.We report on the time evolution of temperature through the particle layer and the products generated by the interactions between hot particles and ice.Numerical simulations of this heat transfer are also presented, along with mathematical analysis at geophysical scales, highlighting the implications and importance of these simulations for natural PDC-ice interactions, and ice-melt lahar generation.

Experiments
A series of static melting experiments were conducted, where hot particles were poured onto a horizontal ice substrate to (a) investigate heat transfer between the particle and ice layers, and (b) quantify melt and steam generation.These experiments were designed as an analog to thermally-driven pyroclast-ice interactions and do not investigate the role of mechanical scouring in incorporating water and steam into hot density currents.Our experiment data are freely available in an online repository (Vale et al., 2023).

Materials
Three particle types were used in the experiments: glass ballotini, crushed pumice (acquired from: Specialist Aggregates, product code: 7803), and an andesitic ash sample from Ruapehu Volcano.Glass ballotini have frequently been used in granular flow experiments because of their highly regular shape and packing structure, making them a good particle type for comparison purposes (Roche et al., 2004;Rowley et al., 2014).The natural samples used in experiments were selected to encompass a range of compositions, from felsic, vesicular volcanic glass (pumice) to more mafic and heterogeneous volcanic samples (Ruapehu ash).Grain characteristics were constrained for the three particle types through image acquisition and analysis techniques, described below (Figure 2, Table 1).
The natural experiment samples were sieved to be within the 500-2,000 μm grain size fractions, while Ballotini particles were pre-sorted into 1,000-1,400 μm sieve fractions.Particle grain size and shape distributions were determined for all three particle types using Dynamic image analysis using a CAMSIZER X2 (Figure 2, Table 1) (Buckland et al., 2021;Microtrac MRB, 2024).Particle sphericity is a measure of shape, computed for an imaged particle as the ratio of the area to the perimeter squared, multiplied by a factor of 4π (Retsch Technology, 2023), which has a maximum value of 1 for a perfect sphere (perfect circle in image, Figure 2).Ballotini particles are spherical, and regular in shape, with sphericities consistently close to 1. Ruapehu ash samples have sphericities ranging between 0.8 and 0.9.Pumice are the most variable in shape, with sphericities between 0.6 and 0.9.
A Hitachi S3500-N scanning electron microscope (SEM) operating in backscattered electron (BSE) mode at the University of Bristol was used to image and characterize Ruapehu ash and pumice samples particle componentry and vesicularity.The samples were mounted in epoxy resin, manually ground to grades PSI 240 to 1200, polished using a Buehler AutoMet 250 autopolisher to grades 9, 3 and 1 μm, and then carbon coated.Pumice samples were composed almost entirely of glass (>95%), with few crystals present.Ruapehu ash samples were more varied in composition and texture, with the presence of microlites, phenocrysts and glass (Figure 2).SEM BSE images were analyzed using ImageJ, an open-source, Java-based image processing software (Schneider et al., 2012).Images were manually edited (i.e., colored black) to remove particle fragments trapped in epoxy in vesicle interiors.Image thresholiding was applied to distinguish between the groundmass and vesicles.The percentage of vesiculated area in the image was calculated using the "analyze particles" function within ImageJ, similar to Liu et al. (2017).
A glass pycnometer was filled with water, and subsequently one type of particles were added to the vessel.The volume of particles was determined by measuring the volume of displaced water (Flint & Flint, 2002).Ballotini particles were the densest, and pumice particles were the least dense.This is consistent with particle vesicularity measurements, where ballotini particles had no vesicles, meanwhile pumice particles had a vesiculated area up to 78%.However, for porous samples these densities are maxima because water may enter some of the vesicles.
The packing fraction of the granular layer in a cylindrical vessel was calculated by dividing the bulk density of the particle layer by the grain density of each particle type.The bulk density was determined by dividing the mass of Note.We also present the corresponding dimensionless parameters for each material using Equations 7a-7d and taking d = 20 mm (b).particles by the volume of the bed.For these packing calculations the grain density was obtained by multiplying the solid density (from He pycnometry on powdered samples) by the average solid fraction of the particles based on vesicle area versus solid fraction analysis of SEM images of polished grains.Analysis of packing fraction showed that packing densities for pumice and ballotini particles were equivalent, while packing density for Ruapehu ash was slightly higher.
Bulk particle thermal properties in a water-saturated state were measured at room temperature using a Portable Electronic Divided Bar (PEDB, product code: Hot Dry Rocks HDR01), at GNS Science, Taupo, New Zealand.
The PEDB determines a ratio between the thermal gradient across the sample and a known material.With this method, thermal conductivity measurements are accurate to within ±3.5% (Antriasian, 2009).Specific heat capacities were also determined by introducing a temperature perturbation and comparing the net thermal energy absorbed by the sample during thermal re-equilibration from one steady-state temperature to another (Antriasian & Beardsmore, 2014).Thermal conductivity measurements of water-saturated experiment particle samples ranged from 0.75 to 1.56 W m 1 K 1 , with pumice particles characterized by the lowest thermal conductivities, and Ruapehu ash particles the highest.

Experiment Configuration and Procedure
Experiments were initiated by releasing hot particles onto a horizontal layer of ice contained within a cylindrical alumina beaker whose initial temperature was approximately 20°C (Figure 3).The particle release lasted ∼2 s, and over that time a layer of particles of uniform horizontal thickness was formed.The particle mass, and hence the layer thickness (5-45 mm), was varied.The nominal initial temperature of particles was varied over the range 200-700°C, though the thermocouples recorded the actual initial temperature following particle layer emplacement.The maximum temperature recorded by each internal thermocouple was used to calibrate our model.The experiment temperature range was informed by PDC temperatures estimated by direct and proxy evidence (Banks & Hoblitt, 1996;Cole et al., 2002;Druitt et al., 2002;Lerner et al., 2019;Scott & Glasspool, 2005).
The evolution of temperature through the particle layer was recorded every second using eight ring-mounted type-K thermocouples.These thermocouples have a temperature recording range from 40 to 1,100 ± 1.5°C, and an experimentally-determined response time of <3 s (RS product code: 397-1236).For all experiments, the thermocouples were fixed at several heights from the ice-particle interface (0 mm) up to 45 mm (surface of the thickest particle layer) (Figure 3).After 10 min of particle-ice contact, the particles were separated from the ice, weighed with the melt, dried, and then reweighed.Melt was calculated as the mass difference between the wet particles (particles plus meltwater) and dry particles.The mass of steam generated was inferred from the difference between the mass loss of ice and total meltwater generated.From each experiment, a single measurement of the total amount of melt and steam produced in the 10 min after the particles first come into contact with the ice was obtained.
A limited set of experiments were also conducted at lower temperatures (20-200°C) and smaller particle depths (10-30 mm) using a larger, rectangular cross-section (330 × 230 mm) apparatus to confirm that the smaller-scale (75 mm diameter) static experiments were scalable to larger systems.In the larger scale experiments, particles were heated in a laboratory oven with a lower thermal range than the furnace used for the small scale experiments.
The oven was used due to capacity constraints.Scaling between experiment configurations is presented in Supporting Information S1 (Vale et al., 2023).

Numerical Simulations
We analyze a one-dimensional domain comprised of hot ash, of thickness d, overlying ice, with initial thickness H 0 (Figure 3).The coordinate system is upwards positive, with an origin defined such that z = 0 marks the initial position of the ash-ice interface.Melting of the ice will move the positions of the ash-ice interface z = s(t), and airash interface z = d + s(t).
In each solid phase, thermal diffusion is described by where j = [A, I] denotes the solid phase for ash and ice respectively; ρ j , c p,j , and k j are the density, specific heat capacity, and thermal conductivity of the jth solid phase respectively.For simplicity, we consider thermal transport in the solid phases only.That is, higher-order thermal effects associated with the imbibition of meltwater are neglected.In Section 4.1, we calibrate our model and demonstrate that this is a reasonable assumption.
During melting, conservation of energy requires that the latent heat of melting (L) is balanced by the difference in heat flux across the ash-ice interface: where T m is the melting temperature of ice and q = k j (∂T/∂z)| z=s is the heat flux, given by Fourier's law, evaluated at the ash-ice interface.This equation-which is often referred to as the Stefan condition-is widely used in moving-boundary problems to describe the velocity of a phase-change interface (e.g., Meirmanov, 2011).Note that due to our sign convention, melting occurs when ds/dt < 0. Therefore, the cumulative melt at time t is given by (ρ w /ρ I )s(t), where ρ w is the water density.
A Dirichlet condition is required to couple the ash and ice subdomains.The magnitude of the interface temperature must account for occurrence, or absence, of melting.Melting is an isothermal process; meaning that the interface temperature is pinned at the melting temperature when Equation 2 is negative.When melting terminates, the heat flux either side of the interface is continuous.Solving Equation 2 with ds/dt = 0 yields a Dirichlet condition for the interface temperature that is a weighted arithmetic mean of ash and ice temperatures either side of the interface.
At the top boundary, heat losses to the air are likely dominated by convection, which are approximated by equating the surface heat flux to a linear constitutive function that is proportional to the temperature difference at the surface: where T air is the ambient air temperature, and γ is a dimensional parameter that control the strength of convective heat losses (e.g., Vollmer, 2009).We assume that the basal boundary is perfectly insulating: This is likely a reasonable assumption for most geophysical settings, where the heat capacity of ice greatly exceeds that of ash (ρ . This is confirmed in Section 4.1, where we demonstrate that our model calibration improves as d/H decreases.In all cases, a uniform initial temperature distribution in the ice (T I (t = 0, z)) is assumed.When calibrating our model at the laboratory scale, we use an initial temperature distribution (T A (t = 0, z)) that is determined experimentally (see Section 4.1).At geophysical-scales (Section 4.2), we assume for simplicity that T A (t = 0, z) is uniform.

Non-Dimensionalization and Rescaling
To better understand the key controls of volcanically-induced ice melting, we reduce the number of parameters in our model by non-dimensionalizing the governing equations and boundary conditions using rescaled variables in combination with characteristic scales Note that the characteristic timescale is diffusive, whereas the characteristic temperature scale is a ratio between the volumetric latent heat of ice to the volumetric heat capacity of ash.
The remaining dimensionless parameters are where α j is the thermal diffusivity of the jth phase.To avoid numerical complexities associated with solving diffusion in a shrinking (ice) domain, we transform our dimensionless coordinate system onto a fixed domain using a bilinear mapping (see Figure 3 and Appendix A).Dropping the tilde notation used above, the resulting remapped non-dimensional governing equations are: Note that rescaling Equation 1 introduces advective terms to account for the motion of the ash-ice interface, such that thermal transport is now described by two coupled advection-diffusion equations with four principle parameters: R α , which compares the strength of thermal diffusion (α j ) in the ice and ash; R k compares the strength of thermal conduction in the solid phases; H compares the initial thicknesses of the solid phases; and the Nusselt number Nu is the ratio of convective to conductive heat transfer at the surface of the ash (η = 1).Note that the range of dimensionless parameters for each material are presented in Table 1.
The coupled system of ordinary and partial differential equations is solved using the method of lines (e.g., Schiesser, 2012).We use a standard first-order finite-volume scheme to discretize our remapped spatial domains; allowing for the resulting system of equations to be expressed as a series of coupled ordinary differential equations (ODEs), which are integrated in time using MATLAB's stiff ODE solver ODE15s (Shampine & Reichelt, 1997).The ODE solver is provided with the pattern of the Jacobian matrix.This significantly reduces the computation time by allowing the solver to only evaluate the Jacobian's non-sparse elements (e.g., Goudarzi et al., 2016).Our model is validated in Section 4.2, where we demonstrate that our numerical scheme agrees with similarity solutions that describe early-time diffusive melting.The convergence of our numerical scheme is verified by exploring the effect of grid resolution on the melting end state s ∞ ≡ s(t → ∞) in the reference case used in Section 4.2.For our simulations, we use 2,500 grid cells in each solid phase.At this resolution, simulations take ∼9 s on a single i7-6500U processor, compared with a run time of ∼65 s for simulations that do not utilize the Jacobian pattern.Further reducing the grid spacing results in variations to |s ∞ | by less than 0.8%.Our numerical solver and associated plotting scripts are freely available in an online repository (Vale et al., 2023).

Heat Transfer From Particle to Ice
The time evolution of temperature at set heights through the particle layer was obtained using eight type-K thermocouples.For all experiments, the thermocouples captured the initial spike in temperature following emplacement and the subsequent cooling of the particles as heat was transferred from the particle layer to the underlying ice and air above.Figure 4 provides examples from 200°C Ruapehu ash experiments as our numerical model is calibrated against this data, however all experiments exhibited the same broad pattern of thermal decay, but over different timescales.The temperature data for all other experiments can be found in Supporting Information S1 (Vale et al., 2023).
The actual maximum particle temperatures recorded by the thermocouples were lower than the furnace temperature because the particles cooled over a period of seconds as they were transported 1.5 m from the furnace, and then poured ∼10 cm into the ice container.Thicker particle layers and particles characterized by lower thermal diffusivities retained more heat in transit and release and so attained higher maximum temperatures.The time at which the maximum temperature was observed varied according to particle type and initial temperature.At higher maximum temperatures the thermocouples took longer to equilibrate with the particles.Ballotini particles reached maximum temperatures fastest, meanwhile pumice particles reached maximum temperatures slowest.This results from the particle's thermal diffusivity.In line with this, ballotini particles also reached thermal equilibrium fastest, and pumice the slowest.Within individual experiments the maximum temperatures obtained by thermocouples varied with position through the particle layer.Typically, the highest temperatures were recorded in the mid-particle region, with cooler maximum temperatures recorded closer to the ice and the particle surface.
The temperature data recorded steam generation signals in two forms, (a) as noise in the thermocouple data, and (b) as a period of stability around 100°C (see Supporting Information S1, Vale et al., 2023).These signals were used to determine the duration and intensity of steam generation in experiments.Of the three particle types examined, pumice particles produced the least steam.This is supported by a general lack of noise or thermocouple stability around 100°C in the temperature data.Ballotini and Ruapehu ash particles produced significant amounts of steam in some experiments.The experiment that produced the most steam, B_M250_T700 (34.81 g), generated steam for over 300 s, or half the total experiment duration, based on duration of the noise signal.
Stepped reductions in the particle temperature followed by stabilizing of the temperature curves were observed in some experiments (e.g., Figure 4c).These steps were initially recorded close to the particle-ice interface, but were subsequently recorded by sequentially higher thermocouples in the particle layer.In the experiments where this behavior was recorded, these stepped reductions in temperature occurred at different rates between successive thermocouples, dependent on the initial experiment temperature conditions.We observed these stepped features across all particle types, most notably in experiments with greater particle layer thicknesses.

Melt and Steam
Systematic increases in melt with increasing particle layer thickness and temperature were observed across all particle types (Figure 5).The amount of melt generated was more sensitive to increasing particle temperature for pumice and Ruapehu ash than for ballotini.The melt data also show a reduction in the rate of increase in melt with increasing particle layer thickness for all particle types.Steam generation also systematically increased with particle mass and temperature for ballotini and Ruapehu ash particles, but was negligible for pumice (Figure 5d).Ballotini particles produced the most steam for any given layer thickness and were more sensitive to particle temperature than the other particle types.
The generation and escape of steam through the upper particle layer also offered insights into the thermodynamics of particle-ice interactions.As expected, steam generation was most productive in the first seconds to minutes following emplacement, where the duration of steam production was dependent on the initial experiment conditions, with the hottest temperatures and thickest particle layers producing steam for the longest durations.In some experiments where steam was produced the particles could be fluidized for up to 10 s.In fluidized experiments, steam escape occurred across the whole surface area of the particle layer.Fluidization of Ruapehu ash particles also resulted in the elutriation of fines from the particle layer in experiments.In some other experiments steam escape via the surface was localized and temporally sporadic.The escape of steam via the upper particle surface also brought melt to the surface with it (see Figure S1 in Supporting Information S1, Vale et al., 2023).

Model Calibration
We calibrate our model using experiments with Ruapehu ash particles heated to 200°C.These experiments provide two sets of measurements that must be approximated by a well-calibrated model: (a) thermal evolution of the ash measured by internal thermocouples; (b) a single measurement of ice melting, recorded after 10 min.In addition to these experimental constraints, our model calibration is further aided by experimental measurements of the specific heat capacity and thermal conductivity of the ash.Moreover, the physical properties of ice are well constrained in the literature.The parameters used in our model calibration are summarized in Table 2.Note that γ is the only free parameter that is undetermined by experimental measurement or literature values.However, we fix γ = 1 W m 2 K 1 as our results are weakly sensitive to typical variations in γ.We explain the physical mechanisms related to this weak sensitivity further in Section 4.2.
Analysis of the experimental thermocouple data demonstrates that a significant amount of heat is lost during transfer of the ash from the oven to the experimental apparatus.This heat loss, which increases for thinner ash layers (see Figure 6), imparts an initial thermal profile in the ash that must be accounted for in an accurate calibration.We implement this in our model using a quadratic thermal initial condition in the ash based on the maximum temperature measured by each internal thermocouple.Specifically, we use an unconstrained multidimensional nonlinear minimization algorithm (Nelder & Mead, 1965) to find the quadratic coefficients that correspond to a minimized total residual between the initial condition and the maximum thermocouple temperatures.This quadratic profile accounts for both preexperiment heat loss and the initial thermal gradients that redistribute heat throughout the system.
We compare our calibrated model with experimental data in Figure 4, which shows the thermal evolution within the ash, and in Table 3, which lists the proportion of melting after 10 min.Note that we present percentage melting values to avoid introducing arbitrary length scales into our 1-D model results.
We find that for all ash thicknesses, our model can accurately simulate the measured thermal behavior of the system.Our simulations accurately recover the magnitude and timescale of heat loss within the ash.Moreover, our calibrated model predicts melting values consistent with those observed in the laboratory.Our calibration performs best for progressively thinner ash layers.The growth of these small errors with ash thickness can be attributed to several effects: (a) maximum experimental temperatures are measured at progressively later times for thicker ash layers.Therefore, our assumed initial condition in the ash (a quadratic fit of maximum thermocouple temperatures) is more appropriate as d decreases.(b) The basal insulating boundary condition (Equation 4) is typically valid provided that (ρ I c I H 0 )/(ρ A c A d) ≫ 1.For thicker ash layers the total thermal capacity of each phase becomes comparable, meaning that boundary effects begin to impact the dynamics of the system.(c) The reduced thermal capacity of thinner ash layers induces less melting.Therefore, we expect that our model, which does not incorporate dynamic effects related to meltwater, to be more valid for thinner layers.(d) Due to the increased total thermal capacity of ash, steam generation increases with d.For ash temperatures below 400°C, steam generation is negligible in Ruapehu ash samples.However at higher temperatures, and when ρ I c I H 0 )/ ρ A c A d) = O(1) or smaller, the latent heat of vapourization becomes non-negligible to the total thermal balance of the system.
Note that the overall melting results in our experiments represent a "snapshot" of the of the melting dynamics.In the model, the end state (i.e. when ds/dt = 0) is defined as s ∞ = lim t→∞ s(t) .For melting to terminate, the ash must lose sufficient heat to balance the flux terms in the Stefan condition.This is expected at the geophysical scale where typically d ≪ H 0 .In this regime, where melting is expected to be negligible relative to the initial thickness of ice (i.e.s ∞ /H ≪ 1), our model calibration performs best.Naturally, this motivates the use of our calibrated   11) are plotted to demonstrate its dependence on R α and R k .In both panels T A T m ≈ 1.73 and T I T m ≈ 0.17.model to evaluate the evolution of potentially hazardous melt generation following the emplacement of hot ash onto ice for natural-scale flows.

Predictions for Melting During Geophysical Scale Processes
We explore a reference scenario using depths of ash deposited on ice as reported by Pierson et al. (1990) and Kilgour et al. (2010).Our reference scenario uses length scales of d = 0.1 m and H 0 = 1 m in conjunction with the calibrated parameters in Table 2 and an initial uniform ash temperature of 200°C.The corresponding dimensionless parameters are listed in Table 4.
Melting occurs in two distinct regimes (Figure 6): early-time diffusive motion of the ash-interface, before a transition to a late-time regime when melting terminates.Both of these regimes are relevant to volcanic hazards: the former governs the rate of meltwater supply; while the latter describes the magnitude of melt generation.In our reference case melting occurs over 0.65-1.56characteristic (diffusive) time scales (dimensional equivalent of 3.48-8.36hr).Note that changing the thickness of ice and/or ash will modify the magnitude and duration of melting, but not the qualitative behavior or existence of these regimes.Before considering the late time regime, we describe the early-time transient behavior.In this regime, melting obeys the classical Stefan problem (e.g., Meirmanov, 2011) and permits the derivation of a similarity solution for the motion of the ash-ice interface: where λ is a constant that determines the early-time melting rate.At early-times, the motion of the ash-ice interface is invariant to the length scales d and H 0 , which allows for the derivation of analytic expressions for temperature in the ash and ice region respectively: where T A and T I are the initial ash and ice temperatures respectively, ξ = z/ ̅ ̅ t √ is a diffusive coordinate transform, and erf(⋅) is the error function.Differentiating these analytical expressions before substituting into the Stefan condition (Equation 8c) yields an expression for λ, which is given by the solution to We solve Equation 11 using MATLAB's nonlinear root finding algorithm fzero.We demonstrate the accuracy of our numerical scheme by overlaying these similarity solutions for different values of R α in Figure 6a.The similarity solution's independence of length scales means that the early-time melting rate is determined by the interplay between R α , R k , T m , T I , and T A only.Solving Equation 11 allows for efficient exploration of this parameter space.Given that the principal temperatures trivially modulate the initial melting rate (e.g., λ monotonically increases with T A and T I ), we consider the impact of R k and R α in Figure 6b.Here we show that λ monotonically increases and decreases with R α and R k respectively.Increasing R k linearly increases the first term in Equation 8c, which reduces the strength of melting by decreasing the heat flux differential across the ash-ice interface.The melting rate increases with R α due to the increased relative strength of thermal diffusion in the ice.Furthermore, by reducing the relative strength of thermal diffusion in the ash, convective heat losses to the atmosphere are transmitted Note.Labels (a), (b) and (c) correspond to the subplots in Figure 4. to the ash-interface at a slower rate, thus delaying the transition from early-time self similar melt propagation to late-time termination of melting (see below).
At intermediate-times, the motion of the ash-ice interface deviates from the early-time s ∝ ̅ ̅ t √ scaling toward a late-time state (s ∞ ) when melting stops.This transition develops as convective heat losses at the upper boundary (η = 1) begin to impact thermal diffusion across the ash-ice interface.Eventually, the ash has lost sufficient heat that the heat fluxes across the ash-ice interface balance and melting terminates.Accordingly the ash-ice interface will remain motionless unless provided with latent heat to resume melting.At intermediate-and late-times, thermal transport and associated motion of the ash-ice interface, is coupled to the inherent length scales of the system, and also to the parameters that govern the rate of heat loss within the system.The impacts of select parameters (R α , H, T A T m , and Nu) on the total (i.e., end state) melt production are illustrated in Figure 7.As in the early-time regime, increasing R α , which strengthens thermal diffusion in the ice, results in increased total melting (Figure 7a).For relevant geophysical settings where H ≫ 1, |s ∞ | is invariant to the ice thickness (Figure 7a).This is because, when melting is small relative to the initial ice thickness, heating of the ice remains confined away from the basal boundary which therefore does not impose any scale effects on the total melt production.In Figure 7b we highlight the impact of the ash temperature and Nusselt number on s ∞ .Intuitively, |s ∞ | monotonically increases with the ash temperature above a threshold temperature which is required to supply latent heat across the ash-ice interface.This threshold is derived by setting λ = 0 in Equation 11, which yields We also demonstrate that, relative to other dimensionless parameters, the Nusselt number has a limited impact on the total melt production.By strengthening convective heat losses to the atmosphere (increasing Nu), melt generation reduces, but only by a small fraction.This relative insensitivity to Nu is because thermal transport in the ash is dominated by conduction into the ice.This is expected for 0 < Nu < 1, and further explains the weak sensitivity to Nu in our model calibration.Note also that as the experiments have not reached end state melting, and therefore as Equation 11is invariant to Nu we expect our calibration to be insensitive to variations in the Nusselt number.In Figure 7b our reference case (Nu ∼ 0.073) will lie between the dashed black and uppermost orange curves.The proximity of our reference case to the convection free (Nu = 0) limit demonstrates that convective heat losses to the atmosphere are essentially negligible over the time scale of melting for typical geophysical parameters and temperatures considered.

Discussion
When PDCs are emplaced onto snow and ice substrates, they rapidly transfer heat from the particle layer into the substrate due to large temperature gradients between the two media.This heat transfer generates melt and steam  12) respectively.
which can be incorporated into the flow transforming its mobility and character.The role of melt and steam in PDC-ice interactions differ, as do their production timescales.Melting is a continuous process for as long as (a) the hot ash can supply latent heat (via a difference in heat flux across the ash-ice interface), and (b) there is a supply of ice to melt.The production and incorporation of melt can cause a PDC to transform from a dry granular flow into a saturated, sediment-laden flow, termed an ice-melt lahar (Thouret et al., 2007).Both melt and steam production depend on the initial temperature gradient and the thermal diffusivity of the particles in contact with the ice but steam production requires heating of water above its vaporization temperature and therefore terminates earlier than melting.The production of steam in PDC-ice interactions can result in elevated pore fluid pressures, fluidization and convective overturning of the particle layer.Fluidization increases the mobility of the layer, while convective overturning sustains the heat transfer from the particle layer to the ice as hot particles from the interior of the particle layer are brought into contact with the ice, essentially "resetting" the thermal boundary condition at the ice-particle interface (Walder, 2000a).Fluidization and convective overturning of the particle layer ceases when the temperature of the particles in contact with the ice drops below the temperature required for steam production.
In natural volcanic settings PDCs scour the ice thermally and mechanically, but before the physical coupling between the thermal and mechanical mechanisms can be considered, the thermal interactions must first be isolated.In the previous sections we presented a series of systematic static melting experiments, along with a calibrated 1-D numerical model and mathematical analysis to resolve the rate of heat transfer from a static hot particle layer to an ice substrate, and to quantify melt.The model can be used to derive a time-series of melt generation (see Section 5.3).From this, ice-melt lahar source hydrographs can be generated, which can be used as an input in surface flow hazard models.From here on, we discuss heat transfer in particle-ice interactions, including the generation and role of melt and steam.Experimental insights and how the 1-D model can be applied at geophysical scales, including constraining the ice-melt lahar hazard are also considered.Finally, the limitations of this investigation are reviewed and recommendations for further work are suggested.

Heat Transfer in Particle-Ice Interactions
Static melting experiments, where a layer of hot particles were emplaced onto a horizontal ice layer were conducted to investigate heat transfer from particle to ice, as an analog to the thermal interactions between hot pyroclasts and ice.Our 1-D model simulates heat transfer between two solid phases, (a) a hot ash layer, and (b) an underlying ice layer.The model is calibrated with experiments using Ruapehu ash heated to 200°C.The experiment data provided two constraints for the model, (a) the time evolution of temperature through the particle layer, and (b) a mass of melt at the end of the 10 min experiment.In all experiments the thermocouples recorded an initial spike in the temperature as particles were emplaced onto the ice, followed by subsequent gradual cooling as heat was transferred from the particle layer into the ice substrate.Our model captures the magnitude and timescale of the heat loss from the ash.With our geophysical extension this model can be used to determine the cooling timescales of deposited pyroclastic material.

Melt and Steam Generation in Particle-Ice Interactions
The heat transfer from the particle layer to the ice substrate initiates the production of melt and steam.The quantities of melt and steam produced in experiments were determined by the particle type, initial particle temperature, and layer thickness.The differing sensitivities of particles to layer thickness and temperature in generating melt and steam are caused by variations in both bulk particle and individual grain characteristics, specifically grain shape, vesicularity, and thermal diffusivity.Ballotini and Ruapehu ash experiments produced comparable quantities of melt, but ballotini particles produced more steam.Pumice experiments produced the least melt and negligible steam at all temperatures investigated.Steam generation was not present in all experiments, and was negligible for all cases where the particles were initially cooler than 400°C.For natural samples with T A (t = 0, z) > 400°C the bulk of ice that is melted or evaporated is measured to be in the liquid phase (Figure 5d); meaning that, melting is the dominant phase transition over the range of temperatures experienced during PDC emplacement.Therefore for simplicity, higher-order terms related to vaporization are not included in our 1-D model.Our numerical model successfully predicts melt generation to c. 5% where steam production is negligible (Table 3).With our geophysical extension, we can predict melt generation at geophysical scales and provide a time-series of this melt generation.This informs ice-melt lahar genesis, which is discussed further in Section 5.3.

Melt
In experiments, melt systematically increased with initial particle temperature and layer thickness.The melt data showed a reduction in the rate of increase in melt with increasing particle layer thickness for all particle types (Figure 5).Two potential explanations for this melt curve flattening are proposed.First, the observed flattening could result from the limited timescale of the experiment such that heat from particles in the upper region of a thick particle layer did not have time to transfer heat to the particle-ice interface.A second possibility is that the observed flattening relates to changes in the partitioning of energy within the system with increasing particle layer thickness.With thinner particle layers the thermal energy melted the underlying ice layer and caused a single phase change from ice to meltwater.With thicker particle layers the increased thermal energy can be expended through (a) heating of meltwater to higher temperatures, and (b) meltwater vaporization.The ratio of this energy partitioning will depend on the initial experiment conditions, including particle layer thickness and temperature.
The particle type (and so particle porosity and thermal diffusivity) is also important; for example, all experiments with pumice particles produced negligible steam.
Temperature data from the experiments also provided insights into the movement of melt through the particle layer.In several experiments, a stepped reduction in particle temperature followed by temperature stabilization was recorded by the thermocouples.This was interpreted to be an upward-moving meltwater front.This meltwater front was cooler than the dry particles, which caused a step in the temperature time-series as the thermocouple came into contact with the melt, and a stabilization as the pore space surrounding the thermocouple became occupied by melt.These stepped reductions in temperature were recorded by sequentially higher thermocouples as the meltwater front rose through the particle layer.These meltwater fronts were observed for all three particle types, and rates of movement ranged from 0.04 to 0.59 mm/s.An example of these stepped temperature time-series can be seen in Figure 4c.As our model does not include a fluid phase, it cannot capture these perturbations.However, our model performs well at reproducing the leading-order thermal decay measured in the experiments.
The principal mechanism proposed for this upward-moving wetting front is particle sinking, displacing ponded ice-surface melt which is generated from the downward wasting of the ice layer.Walder (2000b) similarly reported on the presence of an upward-moving meltwater front as particles passively melted into the underlying snow.This hypothesis was tested for ballotini based on measured packing fraction, where 65% of any given volume is spheres and the remaining 35% is external pore space.The ballotini experiment with a 45 mm particle layer thickness and an initial temperature of 200°C produced 72.1 g melt during the experiment, equivalent to 16.76 mm melt within our experiment apparatus if no particles were present.Taking into account the packing fraction of ballotini for the first 45 mm and conservation of mass, this calculation produces a meltwater height of 46.0 mm, which exceeds the particle layer thickness by 1.0 mm.Melt ponding was observed above the particle surface in this experiment.An additional mechanism for the upward transport of meltwater, relating to steam escape, is discussed in Section 5.2.2.

Steam
For simplicity, higher-order terms related to steam generation are not included in our numerical model.It remains important however to understand the role of steam in PDC-ice interactions, in terms of both flow mobility and character, as PDC emplacement temperatures can be higher than the temperatures required to generate steam, for example, Mount St. Helens (Banks & Hoblitt, 1996), and Soufrière Hills Volcano (Cole et al., 2002;Scott & Glasspool, 2005).Initial particle temperatures in some experiments were sufficiently high to generate steam, revealing a range of additional behaviors pertaining PDC-ice interactions.High temperature experiments revealed the existence of steam-driven melt incorporation, fluidization, and elutriation of fines.Evidence for the presence of steam in experiments was recorded by the thermocouple data in the form of (a) data noise and (b) stability around 100°C.Temperature data noise was interpreted as sporadic and localized steam generation caused by meltwater contacting particles above the minimum temperature required to produce steam.Stability in the temperature time-series around 100°C was interpreted as sustained boiling and generation of steam caused by convective overturning of the particle layer in contact with the ice, resetting the thermal boundary condition.Additional evidence for steam generation was provided through experiment footage (see Supporting Information S1, Vale et al., 2023).
Several features observed during the experiments and in the resulting data provided evidence of steam-driven melt incorporation into the particle layer.Within the first few seconds of some high temperature experiments, we observed vapor-supported melt lenses skittering on the particle surface.After a few seconds the vapor bubbles burst, leaving a saturated area on the particle surface.Thermocouple data indicates that this boiling is occurring at or close to the particle-ice interface as the data show sustained temperature stability around 100°C in the thermocouples closest to the ice.This steam-driven melt mixing observation is also consistent with observations by Walder (2000b), who reported that under some conditions steam generation can cause complete convective overturning of the particles, which drives thermal scouring of the substrate, mixing, and slurry formation.
The thermocouple data provided further evidence for steam-driven melt mixing.In some experiments where steam escaped through the particle layer the temperature profile steps were recorded at multiple successive thermocouples within the same 1 s period, that is, apparent simultaneous temperature reductions recorded by several thermocouples (e.g., Figure S2b in Supporting Information S1, Vale et al., 2023).This sporadic steam escape is caused by dry particles, above the vaporization temperature, coming into contact with cool meltwater.This escape of steam via the upper particle surface also brought melt to the surface with it.This was recorded in experiment footage as localized "flash" wetting of the particle surface (e.g., Supporting Information S1, Figure 3b, Vale et al., 2023).The thermocouple data, in combination with these observations highlight the important role of steam for efficient melt incorporation and mixing within the particle layer.
Granular layer mobility is known to increase when pore fluid pressures are high (e.g., Breard, Dufek, et al., 2020;Druitt et al., 2007;Roche et al., 2002).Measuring pore fluid pressures was not possible in our small-scale experiments, however, fluidization, which occurs as a result of high pore pressure was able to be directly observed for various initial temperatures and layer thicknesses.Fluidization occurs when the upward flux of gas is sufficient to support the weight of the particles above, reducing interparticle contacts, causing the layer to behave in a fluid-like manner (Sparks, 1976).Based on the grain size range of experiment particles, and assuming a density for steam of effectively zero, our experiment particles relate to groups B and D of the Geldart fluidization classification system, which behave as a continuum (Geldart, 1973).Fluidization of the particle layer will have an influence on flow mobility in dynamic particle-ice interaction settings.In experiments fluidization occurred instantaneously following emplacement, and endured for several seconds.Ballotini particles were most readily fluidized across the widest range of initial conditions.This is consistent with the inferred measurements of steam production being greatest for ballotini particles.
Where fluidization occurred in Ruapehu ash experiments, fines were elutriated from the particle layer and spattered up the beaker sides.Elutriation of fines is a phenomenon frequently reported in pyroclastic literature, for example, Wilson (1980), Fisher (1995), and Kelfoun and Gueugneau (2022).Evidence for fines elutriation is recorded in PDC deposits in the form of fines depletion (Brand et al., 2014) and elutriation pipes (Pacheco-Hoyos et al., 2020;Stinton et al., 2014).The presence of this phenomenon in the experiments highlights that the initial particle temperature range in experiments was sufficient to reproduce naturally occurring behaviors.
Steam was not produced in all experiments, nor was it accounted for in our model.Nevertheless, steam plays an important role in the mobility of and melt incorporation into the particle layer due to vapor-driven fluidization and steam escape, which drives the efficient mixing of melt through the layer.Extrapolating these observations to dynamic settings and geophysical scales, it is anticipated that the presence of steam may increase the mobility and potential runout distances of PDCs, and accelerate the transformation from PDC to ice-melt lahar.
Ice-melt lahars have historically represented some of the most hazardous volcanic flows (Brown et al., 2017).For example, during the 1985 eruption of Nevado del Ruiz, Colombia, PDCs were emplaced onto the summit area resulting in thermal and mechanical scour of the snow and ice substrate, removing 16% of the surface area and 9% of the total volume of the summit ice cap (Thouret, 1990).The incorporation of melt, snow, and ice into these PDCs transformed them into ice-melt lahars that flowed down valleys on the volcano flanks.The propagation of these lahars through populated areas resulted in ∼25,000 fatalities (Naranjo et al., 1986), highlighting the extreme hazards posed by PDC-ice interactions and ice-melt lahars, and the need for robust modeling of these events.
Our experiments and numerical model provide constraints on the magnitude and timescales of melting in PDC-ice interactions.This can be used to inform source conditions for the generation of ice-melt lahars.Physics-based simulations of lahars (and other debris flows) used in hazard assessment typically use volumetric flux [L 3 /T] hydrographs as source conditions, which provide time-series fluxes of water and entrained solids (e.g., Jenkins et al., 2023).This source flux is typically distributed over an area A [L 2 ], meaning that the prescribed flux has units [L/T].In cases of deposition of hot, static ash, our model can provide a time-series of melt generation.Given that s(t) represents the total melt [L] generated in a time t, the meltwater flux q w [L/T] is given by where ρ w is the water density.In Section 4.2, we established that the early-time motion of the ash-ice interface obeys s ∼ λ ̅ ̅ t √ .Substituting this into Equation 13 and taking ρ w /ρ I ≈ 1 yields the early-time form of the meltwater flux, which obeys q w (t) ∼ λ/ (2 In Figure 8 we plot redimensionalized hydrographs that correspond to the simulations shown in Figure 6.Here, the meltwater flux decay follows the q w ∝ 1/ ̅ ̅ t √ early-time scaling, before sharply dropping off to zero as melting terminates.Our model predicts that melting occurs over nearly 10 hr, and the magnitude and duration of q w exceeds well-established empirical thresholds for analogous rainfall-driven debris flows (e.g., Guzzetti et al., 2008).By combining hydrographs generated by our model with calibrated empirical thresholds for lahar/debris flow initiation, insights can be provided into the triggering conditions of melt-driven lahars.Furthermore, where melt-driven lahar genesis is expected, our model provides a physical basis for the general form of melt-driven source hydrographs.Note that a slight perturbation is required to avoid a singularity in q w at t = 0.However, morphodynamic flow solvers typically require a short "ramp-up" period to avoid instabilities in the source region.Above, homogeneous ash emplacement, distributed over an area A, is considered.However, in natural settings, the thicknesses and thermal properties of the ash and ice will vary laterally.Given that the lateral extent of ash greatly exceeds the ash thickness (L/ d ≫ 1), melting dynamics and overall heat loss is dominated by vertical thermal transport.Therefore, it would be a reasonable assumption to neglect terms related to lateral thermal transport, and instead divide the system into a series of 1-D subdomains.

Limitations and Recommendations for Further Work
In the static melting experiments, the interactions between hot granular media and ice are investigated.In nature, where PDCs are emplaced onto frozen substrates, the substrate surface will typically be comprised of varying thicknesses of less-dense snow, and underlain by ice.Ice is used in our experiments for two reasons: (a) snow is difficult to reproduce in a laboratory environment, and (b) ice is a less complex substance and efforts were made to simplify the research problem for modeling purposes.Further research into PDC interactions with frozen substrates should consider how heat transfer, melt and steam generation would differ if the substrate consisted of snow, rather than ice.The experimental works of Walder (2000aWalder ( , 2000b) investigated hot particle-snow interactions and provided useful insights into thermal scour of snow substrate by convective vapor bubbling, and passive melting at lower temperatures.However, this research did not quantify the rate of heat transfer, or the generated melt and steam, limiting comparison with our numerical model.
Similar to Walder (2000b), our experiments isolated the thermal interactions between particle and ice, and did not account for mechanical scouring, which can occur in natural PDC-ice interactions.The static experiments presented here represent an important precursory investigation to dynamic experiments from phenomenological and numerical perspectives because in nature some emplaced particles will form a static bed, and therefore some passive melting will occur.Furthermore, prior to these experiments the dynamics and timescales of melt and steam generation were unknown.In nature where mechanical scour processes also operate, mechanical abrasion of the substrate can generate ice fragments which can be incorporated into the hot granular layer.This may enhance the production of meltwater as heat transfer from the hot granular layer to fragmented ice will be more efficient, compared with unfragmented ice, due to the larger surface area to volume ratio of the ice fragments (Stroberg et al., 2010).To investigate the effects of coupled thermal and mechanical scouring, hot granular layer interactions with ice should be studied in dynamic flowing configurations.
Though higher-order terms related to steam generation were not included in our numerical model, it was nonetheless able to accurately reproduce the leading order dynamics measured in the experiments.This is not considered an insurmountable limitation as there are currently no means to quantify steam in natural PDC-ice interactions, and steam generation at geophysical scales using the ratio of melt to steam generation in the small-scale laboratory experiments under different initial conditions can be estimated.At geophysical scales, it is most critical that (a) melt volume generated from pyroclast-ice interactions can be quantified, and (b) the physical role of steam for efficient melt incorporation, and its effects on flow mobility are understood.The static melting experiments in combination with the numerical model work to resolve these requirements.

Conclusions
A series of static melting experiments were conducted, where hot particles were emplaced onto an ice substrate as an analog to PDC-ice interactions.These experiments isolated the thermal interactions, as in order to fully understand the thermal and mechanical coupling in PDC-ice interactions, a detailed physical knowledge of particle-ice interactions in the simplest configuration must first be generated.The experiments revealed that melt and steam systematically increase with increasing particle temperature and layer thickness.The experiments were capable of reproducing phenomena seen in natural-scale flows, including fluidization and fines elutriation, indicating that the initial temperature conditions were within a representative natural range.Experiments also provided insights into melt mixing mechanisms.Based on the thermocouple data, in combination with visual observations, steam was identified to play critical role in the rate of melt incorporation through the particle layer.This has implications for the rate of transformation from PDC to ice-melt lahar.
A 1-D numerical model of heat transfer, which was calibrated against Ruapehu ash experiments is also presented.This model accurately predicts melt generation to within ∼5%.We also provided analytical similarity solutions for our numerical model at early-times and at typical geophysical scales, along with an example meltwater flux hydrograph which can be used to inform source conditions for simulations of melt-driven lahars.The ability to predict melt generation at geophysical scales when a PDC is emplaced onto an ice substrate represents a significant advancement toward robust modeling of the ice-melt lahar hazard.

Appendix A: Bilinear Remapping of the Numerical Model
A common approach in obtaining numerical solutions to moving boundary problems is to rescale the governing equations onto a fixed domain.Although introducing additional mathematical complexity (in the form of advective transport terms), solving advection-diffusion problems on a fixed domain removes numerical complexities associated with solving coupled PDEs on evolving domains (e.g., time-dependent grids and resolution).Bilinear mapping is adopted.That is, a separate linear transform is applied to each solid phase.The location of key interfaces in these rescaled domains (η for ash, ν for ice) are summarized in the table below.These rescaled domains are solved separately and coupled together through a shared Direchlet boundary condition for temperature and the Stefan condition, which scale advective terms that account for the motion of the ash-ice interface.

A1. Remapping the Ash Subdomain z⇔η
The ash subdomain is rescaled using the linear mapping η = z s. (A1) Accordingly, the spatial and temporal derivatives from independent variables (z, t) to (η, t) must be transformed.The transformed first and second spatial derivatives are and Finally, the transformed temporal derivative is given by

A2. Remapping the Ice Subdomain z⇔ν
The same procedure is now performed in the ice region, which is rescaled using the linear mapping: whose first and second spatial derivatives are and The transformed temporal derivative is given by Applying these transformed spatial and temporal derivatives to each subdomain yields our system of rescaled, non-dimensional governing Equations 8a-8d.

Figure 1 .
Figure 1.Redoubt volcano viewed from the northwest following the 4 April 2009 eruption.Annotations show volcanic processes and deposits.Incisions in the glaciated surface indicate thermal and mechanical scour by pyroclastic density currents (PDCs).Lahar deposits are observed downstream of the PDC deposits.Photo source: USGS (2009).

Figure 2 .
Figure 2. Grain size distributions of experiment particle types, glass ballotini (left column), crushed pumice (central column), and Ruapehu ash/lapilli (right column), shown with Camsizer X2 imagery for shape analysis and Microscope/scanning electron microscope (SEM) imagery used for componentry and vesicularity analysis.SEM Voltage 20 kV, and working distance 22.4 mm.

Figure 3 .
Figure 3. Experimental configuration and model domains.(a) Particles are poured through a thermocouple array onto a horizontal layer of ice.Temperature is recorded using a ring-mounted thermocouple array supporting eight type-k thermocouples set at different heights through the particle layer.Thermocouple heights are denoted by star markers.(b) Schematic diagrams of our numerical domains.For numerical convenience, we non-dimensionalize and rescale our dimensionless coordinate system (left) into separate subdomains (right) for each phase (see Appendix A).

Figure 5 .
Figure 5. Panels (a-c) Ice Mass Loss as a function of the mass of heated particles of different initial temperatures.The ice and particle masses are normalized by their horizontal surface area in contact.Panel (d) Total Meltwater against Ice Mass Loss, where deviation from a 1:1 slope equates to steam.Particle types are denoted by symbol and correspond to the same symbol shapes as panels (a-c).

Figure 6 .
Figure 6.(a) Evolution of the ash-ice interface plotted for various values of R α .The interface initially obeys the early-time similarity solution (dashed blue curves) s = λ ̅ ̅ t √ before transitioning to an end state where melting terminates.(b) Contours of λ (Equation11) are plotted to demonstrate its dependence on R α and R k .In both panels T A T m ≈ 1.73 and T I T m ≈ 0.17.

Figure 7 .
Figure7.End state position of the ash-ice interface versus (a) R α , with several different thickness ratios plotted to demonstrate the invariance of total melting to H; (b) Difference in ash and melting temperature T A T m , also plotted for various values of Nu which highlights the weak dependence on the Nusselt number.The dashed black and blue curves represent the Nu = 0 limit and melting threshold (see Equation12) respectively.

Figure 8 .
Figure 8. Example hydrographs calculated by dimensionalizing the curves in Figure 6a.Note that early-time melting produces a meltwater flux that decays with q w ∼ 1/ ̅ ̅ t √ .

Table 1
Grain Characteristics of the Particle Types Used in Static MeltingExperiments (a)

Table 2
Calibrated Model Parameters Note.c p,A and k A are mean values from five laboratory measurements.

Table 3
Measured and Simulated Melt After 10 Min