Proc R Soc A

volcanology, mathematical modelling, asymptotic methods

Emma Greenbank

# A theoretical model of Surtseyan bomb fragmentation

Emma Greenbank1 Mark J. McGuinness1 C. Ian Schipper2 1School of Mathematics
and Statistics, Victoria University of Wellington, New Zealand.
2School of Geography, Environment and Earth Sciences, Victoria University of
Wellington, New Zealand<EMAIL_ADDRESS>
###### Abstract

Surtseyan eruptions are an important class of mostly basaltic volcanic
eruptions first identified in the 1960’s, where erupting magma at an air-water
interface interacts with large quantities of slurry, a mixture of previously
ejected tephra that re-enters the crater together with water. During a
Surtseyan eruption, hot magma bombs are ejected that initially contain pockets
of slurry. Despite the formation of steam and anticipated subsequent high
pressures inside these bombs, many survive to land without exploding. We seek
to explain this by building and solving a simplified spherical mathematical
model that describes the coupled evolution of pressure and temperature due to
the flashing of liquid to vapour within a Surtseyan bomb while it is in
flight. Analysis of the model provides a criterion for fragmentation of the
bomb due to steam pressure buildup, and predicts that if diffusive steam flow
through the porous bomb is sufficiently rapid, the bomb will survive the
flight intact. This criterion explicitly relates fragmentation to bomb
properties, and describes how a Surtseyan bomb can survive in flight despite
containing flashing liquid water, contributing to an ongoing discussion in
volcanology about the origins of the inclusions found inside bombs.

###### keywords:

Surtseyan eruption, magma bombs, magma fragmentation, flow in porous media

††journal: rspa

##  1 Introduction

Surtseyan eruptions are characterised by significant bulk interaction of water
with ascending magma [1, 2, 3, 4, 5], as exemplified by the volcano off the
coast of Iceland that rose above the sea to become the island of Surtsey in
the 1960s [1, 6].

Like Hunga Ha’apai in Tonga in 2009, and Copelinhos in the Azores from
1957–1958, volcanoes that have risen under the sea and have a crater rim that
is at sea-level may exhibit this style of eruption. Mt Ruapehu in New Zealand
with its crater lake sometimes erupts in this way.

Magma-water interaction plays an important role in determining volcanic
eruption styles [3]. The style of interaction can vary greatly, from deep
submarine scenarios where quenching of magma causes relatively passive thermal
granulation, to terrestrial scenarios where rapid heat exchange between magma
and ground or surface water can enhance magma fragmentation and explosivity
[7]. Surtseyan eruptions, initially thought of as shallow submarine events,
are now understood to involve interaction with any shallow standing water
body, including lakes, rivers and marine waters [5].

Whatever the conditions under which magma and water meet, their interaction is
limited in large part by the extreme difficulty of mixing fluids with such
strongly contrasting viscosities and other physical properties [8]. Surtseyan
eruptions, illustrated in the cartoon in Fig. (1), are unique, in that the
water-magma interaction occurs in near-surface, periodically-flooded vents,
where there is re-entry of in-vent water-saturated slurry, a mixture of water
and previously erupted ejecta that has fallen back into the crater [2, 4].
Fresh magma that rises to encounter this slurry can more readily mix with it
than with pure water, due to the slurry’s higher bulk viscosity and lower heat
capacity [9]. The almost silent jets of tephra that are uniquely
characteristic of Surtseyan eruptions have shapes that have been compared to
cypress trees and cocks’ tails, and are a direct result of magma-slurry
interaction [1].

The component of the ejecta that we model here is the Surtseyan bomb, a lump
of magma which trails a black tail of tephra as it shoots out of an erupting
tephra mass, hence appearing at the tips of the cypress tree tephra shapes
illustrated in Fig. (1). This black tail turns white within seconds, due to
cooling and condensation of superheated water vapour within it [1, 4].
Examination of bombs collected from Surtseyan deposits [4] reveals they are
vesicular (bubbly) with porosities ranging from 0.35 to 0.8, and they are
hosts for entrained or engulfed smaller groups of clasts. The vesicularity is
due to the usual process of water coming out of solution in the magma due to
pressures dropping as the magma rises in the throat of the volcano. Samples of
bombs are shown in Fig. (2), and there is a video in the repository referenced
in [4] showing multiple inclusions in virtual x-ray sections through a bomb.
Each enclosed clast is usually found to be loosely held in a void space in its
host bomb. The first few millimetres of the material surrounding a void shows
evidence of compression and densification, suggestive of high historical
pressures [4], and indicative of the magma having been above the glass
transition temperature for some time after slurry entrainment [10]. The bombs
are permeable to fluid flow, due to their high vesicularity [11].

Figure 1: Cartoon of a Surtseyan eruption. Numbers refer to the water vapour
cloud (1), ash in cypress tree shapes with bombs at their tips (2), crater
(3), abundant water at crater level (4), layers of lava and ash (5), stratum
(6), magma conduit (7), magma chamber (8), dike (9). Copyright: Creative
Commons license, https://commons.wikimedia.org/wiki/File:Surtseyan_Eruption-
numbers.svg, attributed to ©Sèmhur / Wikimedia Commons / CC-BY-SA-3.0 (or Free
Art License) .

Figure 2: Composite clasts from Surtsey (Iceland) Island, showing entrained
smaller clasts, near the outer surface or after breaking off a piece of bomb.

If a volume of liquid is flashed instantaneously to steam while constrained to
that same volume, the resulting pressure [12] far exceeds the tensile strength
of the surrounding vesicular magma, and the bomb is then expected to fragment.
However, observations indicate that many bombs survive intact, suggesting that
while steam is created by heat transfer from magma, raising the pressure
inside a bomb, it can also escape through the surrounding vesicular bomb,
relieving the pressure increase. So there is a race between heating that
creates steam to raise the pressure and steam flow that relieves the pressure
inside a bomb. Resolving the winner of this race provides the motivation for
the development of a model for the transient behaviour of pressure and
temperature inside Surtseyan bombs while they are in free fall after ejection,
to find the conditions under which a bomb is expected to survive the
development of high pressures around a flashing slurry inclusion.

Other possible fragmentation mechanisms that we do not address here include
impact with the ground, and vesiculation of the interior together with
quenching of the outer surface.

The most relevant previous work modelling the flashing of liquid inside
vesicular magma is a relatively simple model [12] that takes a cavalier
approach in ignoring temperature transients, despite their importance in
driving up the pressure. By ignoring the possibly very large initial
temperature gradients that drive the flashing of liquid to steam, that work is
likely to lead to a serious under-estimate of the maximum pressures that
arise. The main purpose of the present work is to properly include the driving
mechanism of transient temperature gradients at the boiling front surrounding
a slurry inclusion, and hence provide more reliable estimates for the maximum
pressures that are developed in the inherently transient interaction between
heating and steam escape.

There is no other directly comparable previous work on modelling pressure
increase inside Surtseyan bombs that we are aware of. The style of modelling
here is inspired by that used to explain shock tube experiments [13] in which
volcanic rock samples are fragmented by pressure differences that propagate
into the samples as high pressure gas escapes from them. However, in that
modelling and experiments the thermal behaviour is adiabatic, in contrast to
the transient energy conservation principles needed here.

A similar approach can be seen in recent modelling of the flashing of liquid
when frying potato snacks [14]. The flashing of liquid water to steam is
central to that modelling, but the development of high pressures is not of
interest due to the relatively large permeabilities of the cooked snack, so
that they are able to use a steady-state diffusion equation for vapour flow
outwards from the flashing front.

As in the model presented in [12], we make no attempt here to model the
development of vesicularity in rising magma, or to consider viscous flow
effects or changes in porosity due to compression. Instead, vesicular molten
magma is treated as a competent porous rock matrix of constant porosity, due
to the high viscosity of the magma at temperatures at and below 1275K [10, 15,
16]. The Deborah number

${\rm De}=\frac{t_{\rm vis}}{t_{\rm elas}}$

may be used to separate the time scales on which viscous (flowing) and elastic
(brittle) responses occur in magma under stress. Here $t_{\rm vis}$ is the
time required for viscous relaxation of magma, and $t_{\rm elas}$ is the time
needed to deform magma as a competent solid. Brittle responses occur for De>1
[16], and viscous flow is more important for relieving stress if De<0.01.
Viscous relaxation times are strongly dependent on temperature, and at 1275K
$t_{\rm vis}\sim 1$ s. Then the assumption of brittle behaviour is good for

${\rm De}\sim\frac{1}{t_{\rm elas}}>1\;,$

that is, for times $t_{\rm elas}<1$ s. It will be seen in the numerical
simulations that maximum pressure is reached in a millisecond or less,
consistent with this assumption.

In the remainder of this paper, we present new measurements of porosity and
permeability in intact Surtseyan bombs, and a regression relationship between
them. We then build a new fully transient model from scratch by starting with
physically accurate and properly coupled conservation equations in the vapour
and liquid regions, with a moving boiling front between them. Once we have
obtained a consistent set of coupled nonlinear partial differential equations,
we nondimensionalise, choosing appropriate scalings so that the essential
transport mechanisms are captured in a reduced set of equations. A particular
focus in this paper is to then explore the theoretical consequences for the
maximum pressure developed in the model, if the initial temperature profile is
allowed to have a step change at the flashing front. We combine numerical
solutions with asymptotic arguments that suggest the step change initial
temperature case is mathematically ill-posed. An initial temperature profile
that is ramped from cold to hot over a distance the size of a typical pore
provides a criterion for fragmentation in terms of magma permeability or
porosity. The criterion indicates that the permeabilities of the intact bombs
are all high enough that the maximum pressures reached are not expected to
fragment the bombs.

##  2 Field Work and Data

Porosity and permeability were determined on small (<64 mm diameter) ejecta
collected from the Island of Surtsey, Iceland [17], and is available at
https://doi.org/10.5281/zenodo.5214799. Deposits formed during typical
Surtseyan tephra jetting were examined in the field in 2015. The deposits
consisted of poorly sorted discontinuous beds, 10–20 cm in thickness. They
were dominated by the fine ash (<2 mm) that typifies the products of
phreatomagmatic eruptions [19], but also included abundant lapilli (2–64 mm)
and strikingly outsized bombs (>64 mm).

Bombs were universally observed to be composites, containing identifiable
remnants of slurry incorporation (Fig. (2)). The slurry inclusions are usually
of a relatively large grainsize, but there is no reason to believe that the
overall particle size distribution of the slurry itself should be any
different to that in the preserved deposits. Our porosity and permeability
data is from lapilli within this deposit, as they are the most representative
of portions of vesicular magma that evaded fragmentation, and they retain the
textural characteristics of the magma as it erupted. Ash particles are too
small to retain useful information about pore networks, and bombs have greater
potential for post-eruptive expansion that can overprint the syn-eruptive
magmatic textures that are of interest here [20].

The lapilli were examined using X-Ray computed Microtomography ($\mu$-cT).
Each sample was scanned with a Phoenix Nanotom 180 X-ray $\mu$-cT at
l’Université d’Orléans. Totals of 2000–2300 scans of each sample were
collected during 360∘ rotation, using a tungsten filament and molybdenum
target. Operating voltages were in the range 80 to 100 keV, with currents of
50–90 nA. Voxel edge lengths ranged from $\sim$2.5 to 6.5 $\mu$m. Raw scans
were reconstructed into stacks of greyscale images with an offline PC
microcluster running Phoenix reconstruction software. Quantification of
porosity and modelling of permeability were performed on representative
elementary volumes of 500–750 px3 isolated within each $\mu$-cT scan. Multiple
sub-volumes capturing the range of vesicle heterogeneity in natural samples
were isolated and analyzed separately, and all isolated sub-volumes were
chosen to avoid any entrapped slurry or other inclusions that were not
representative of the magmatic foams themselves. Vesicles were measured using
the 3D object counter plugin for ImageJ [21], after greyscale thresholding to
isolate void space from glass and phenocrysts.

We modelled permeability using a gas flow simulation program [22]. This
parallel computing program measures Darcian permeability by simulating single-
phase gas flow with lattice Boltzmann simulations on the Palabos computational
fluid dynamics platform. For each sub-volume, three calculations were made in
three orthogonal directions. Simulations were conducted using low inlet
pressures, ensuring low Reynolds number flow conditions, thus neglecting
inertial contributions to total permeability [23]. The investigated samples
are marginally smaller than the bombs described in this work and shown in Fig.
(2), but have similar matrix porosities, and are considered to be
unfragmented, intact bodies that are representative of Surtseyan magma at the
time of eruption.

These calculations of permeability and porosity of the vesicular network
within intact, unfragmented samples of Surtseyan ejecta are summarised in Fig.
(3), together with the fitted straight line

$\log_{10}k=6.4\phi-14.1\;.$ (1)

This relationship between permeability and porosity is similar to that
observed in other studies of vesicular basaltic magma [11], and is used later
when analysing model behaviour. Even though some vesicles adjacent to slurry
domains in composite bombs appear modestly compressed, we consider this to not
devalue the use of the porosity-permeability relationship shown in Figure 2,
because the compression occurs in a narrow halo around slurry domains, and
because the hysteresis effect on permeability in vesicular magmas ensures that
modestly compressed bubble networks will retain significant permeabilities
developed during their previously uncompressed states [23].

Figure 3: The log of permeability plotted against porosity, as measured in
samples of Surtseyan ejecta, together with a best fit straight line.

##  3 Mathematical Model

We model the boiling of the water at the surface of a single slurry inclusion,
due to the surrounding relatively hot vesicular magma, and the subsequent
transport of vapour out of the magma. We focus on the pressure changes in
connected pores consequent on boiling, in an approach which is simpler than
that in [13], positing that an excessive pore pressure will fragment the
surrounding magma, and making no attempt to model stress and strain in the
magma itself.

We represent both slurry and magma as nested equivalent spheres of solid
porous material. We place the inclusion at the centre of the magma and set
this centre to be origin, giving a frame of reference that is travelling with
the bomb. The entire Surtseyan bomb is ejected from the volcanic vent and is
travelling through the air in free fall for the first several seconds of its
existence, so gravity effects are ignored. We neglect any effects from
possible rotation of the bomb, or air friction on the outside of the bomb. The
slurry inclusion is modelled as a porous medium filled with liquid water. The
slurry temperature is assumed to be initially near boiling point at
atmospheric pressure. Our simulations indicate that initial slurry temperature
has little effect on pressure development compared to the magma temperature.
The enveloping hot magma is treated as a porous medium with its connected
porosity containing only steam. A two-phase region will develop at the
interface between liquid and steam phases. We assume, supported later by the
large size of the Stefan number, that the two-phase region is relatively thin.
The entire bomb is assumed to start at atmospheric pressure.

The pointwise conservation of enthalpy equation for a moving fluid (which here
may be the liquid or vapour phase of water) in the absence of sources or
sinks, takes the form [24]

$\frac{\partial}{\partial t}(\varrho h)+\nabla\cdot(\varrho h{\bf
v})=-\nabla\cdot{\bf q}+\frac{Dp}{Dt}+\tau:\nabla{\bf v}$ (2)

where $\varrho$ is fluid density, $h$ is specific enthalpy, $p$ is fluid
pressure, $\bf v$ is the local fluid velocity vector, $\bf q$ is the heat
flux, and $D/Dt=\partial/\partial t+{\bf v}\cdot\nabla$ is the total
derivative operator. The symbol $\tau$ is the deviatoric stress tensor in the
viscous dissipation term.

As noted in [25] and in [26, Appendix D], the pointwise conservation of
enthalpy equation can be written in terms of pressure and temperature $T$ as

$\varrho c_{p}\frac{DT}{Dt}-\beta T\frac{Dp}{Dt}=\nabla\cdot(\kappa_{T}\nabla
T)+\tau:\nabla{\bf v}\;,$ (3)

where the specific heat at constant pressure is

$c_{p}=T\left(\frac{\partial S}{\partial T}\right)_{p}\;,$

the specific enthalpy is $h=U+p/\varrho$, $S$ is the specific entropy, and $U$
is specific internal energy. We have also used $dh=TdS+dp/\varrho$. and
$TdS=c_{p}dT-\frac{\beta T}{\varrho}dp$ to obtain equation (3). The
coefficient of isothermal compressibility is

$\beta=\varrho\left(\frac{\partial(1/\varrho)}{\partial p}\right)_{T}\;.$

We have also used Fourier’s law for heat conduction ${\bf q}=-\kappa_{T}\nabla
T$, where $\kappa_{T}$ is the thermal conductivity of the fluid, and the mass
conservation equation for fluid,

$\frac{\partial\varrho}{\partial t}+\nabla\cdot(\varrho{\bf v})=0\;.$ (4)

When considering the rock component of the porous medium that we take to
comprise both the slurry inclusion and the surrounding hot magma, we ignore
small displacements, velocities, and compressibility, and we consider it to be
a competent solid material. We assume our representative elementary volume is
small enough that there are no appreciable local temperature differences
between rock and steam, so that after averaging as in [24, 26] over a
representative elementary volume of porosity $\phi$, we have conservation of
energy in the rock component as

$\varrho_{m}c_{pm}(1-\phi)\frac{\partial T}{\partial
t}=(1-\phi)\nabla\cdot(\kappa_{m}\nabla T)\;,$ (5)

where $\varrho_{m}$ is rock density in the absence of any porosity, $c_{pm}$
is the specific heat capacity of the rock, and $(1-\phi)\kappa_{m}$ is its
effective thermal conductivity.

We now consider separately the liquid-filled slurry region $r<s(t)$ inside the
flashing front, and the vapour-filled hot magma region $R_{2}>r>s(t)$ of hot
magma outside the flashing front. The only distinction we make between slurry
and hot magma is that cool liquid occupies the porosity in the slurry region
and hot vapour occupies the porosity in the magma region. For computational
simplicity here, but also consistent with our observations in section 2 that
inclusions have similar properties to surrounding rock, we assume the same
(constant) porosity $\phi$ in both regions.

###  3.1 Slurry Region

Conservation of mass for liquid water in the slurry region $0<r<s(t)$ is given
by

$\phi\frac{\partial\varrho_{l}}{\partial t}+\nabla\cdot(\varrho_{l}{\bf
u}_{l})=0$ (6)

and the momentum conservation equation for liquid is given by Darcy’s Law, as
liquid velocities arising from the compressibility of liquid water are found
to be small:

${\bf u}_{l}=-\frac{k}{\mu_{l}}\nabla p\;,$ (7)

where $p$ is the pressure in the fluid, and we do not model pressure in the
rock matrix. An equation of state for liquid density is provided by the
expression

$\varrho_{l}=\varrho_{0l}+\beta_{\ell}(p-p_{a})-\alpha(T-T^{w}_{0})\;,$ (8)

which fits within a 10% accuracy the specific volume data present in [27], for
pressures in the range 1-100 bars, and temperatures between saturation and 573
K. The terms $\varrho_{0l}$, $p_{a}$, and $T^{w}_{0}$ are reference values of
liquid density, pressure and temperature, listed with $\alpha,\beta_{\ell}$ in
Table 1.

Averaging the point-wise energy equation (3) over a representative elementary
volume as in [24, 26], the energy equations for liquid water and solid rock
combine to give

$\varrho^{\prime}c^{\prime}\frac{\partial T}{\partial t}+\varrho_{l}c_{pl}{\bf
u}_{l}\cdot\nabla T-\phi\beta T\frac{\partial p}{\partial t}-\beta T{\bf
u}_{l}\cdot\nabla p=\kappa_{el}\nabla^{2}T+\overline{\tau:\nabla{\bf u}_{l}}$
(9)

where the temperature $T$ is assumed to be the same in fluid and adjacent
rock, $\varrho_{l}$ is the density of liquid water, ${\bf u}_{l}=\phi{\bf
v}_{l}$ is the Darcy velocity of the liquid water,
$\kappa_{el}=(1-\phi)\kappa_{m}+\phi\kappa_{l}$ is the average thermal
conductivity of the liquid-filled region, $c_{pl}$ is the specific heat of
liquid water, and $\kappa_{l}$ is the thermal conductivity of liquid water.
Typical values for these parameters are given in Table 1. The symbol $\tau$ is
the deviatoric stress tensor in the averaged viscous dissipation term,
accounting for heat generated in liquid water due to visous shear flow within
pores.

Table 1: Physical Constants Constant | Name | Typical Value | Units
---|---|---|---
$c_{pl}$ | specific heat of liquid water | 4200 | J.kg-1.K-1
$c_{m}$ | specific heat of magma | 840 | J.kg-1.K-1
$c_{ps}$ | specific heat of steam | 2000 | J.kg-1.K-1
$h_{sl}$ | specific heat of vapourisation | $2.3\times 10^{6}$ | J.kg-1
$k$ | permeability | $10^{-12}$ | m2
$\kappa_{e}$ | thermal conductivity [28] | 2 | W.m-1.K-1
$\kappa_{el}$ | thermal conductivity [28] | 3 | W.m-1.K-1
$M$ | molar mass of water | $18\times 10^{-3}$ | kg.mol-1
$p_{a}$ | atmospheric pressure | $10^{5}$ | Pa
$R$ | universal gas constant | 8.314 | J.K-1.mol-1
$R_{1}$ | initial inclusion radius | 0.01 | m
$R_{2}$ | magma (bomb) radius | 0.1 | m
$T_{i}$ | initial inclusion temperature | 373 | K
$T_{m}$ | initial magma temperature | 1275 | K
$T_{c}$ | critical temperature of water | 647 | K
$T_{0}^{w}$ | reference temp, $T_{i}$ |  |
$T_{0}^{e}$ | reference temp, $T_{i}$ |  |
$\alpha$ | thermal exp coeff for liquid water | 0.5 | kg.m-3.K-1
$\beta_{\ell}$ | isothermal comp. of liquid water | $4.6\times 10^{-10}$ | kg.m-3.Pa-1
$D_{m}$ | thermal diffusivity $\kappa_{e}/(\varrho c)$ in magma | $1.4\times 10^{-6}$ | m2.s-1
$\mu_{v}$ | dynamic viscosity of water vapour | $3\times 10^{-5}$ | Pa.s
$\phi$ | porosity | 0.3 |
$\varrho_{m}$ | density of basalt | 2750 | kg.m-3
$\varrho c$ | $(1-\phi)\varrho_{m}c_{m}+\phi\varrho_{s}c_{ps}$ | 1.6 $\times$106 | J.m-3.K-1
$\varrho^{\prime}c^{\prime}$ | $(1-\phi)\varrho_{m}c_{m}+\phi\varrho_{l}c_{pl}$ | $3\times 10^{6}$ | J.m-3.K-1
$\varrho_{0l}$ | reference liquid density | 1000 | kg.m-3

###  3.2 Magma Region

In the region $R_{2}>r>s(t)$ there is only steam in the pores, mostly above
the critical temperature. Initially this steam would be magmatic vapour from
the vesiculation process that created the pores, which we are not modelling.
It will be displaced by steam generated by flashing the liquid in the slurry,
the process we are modelling here. Mass conservation of steam is given by

$\phi\frac{\partial\varrho_{s}}{\partial t}+\nabla\cdot(\varrho_{s}{\bf
u})=0\;.$ (10)

The momentum conservation equation for fluid flow is typically given by
Darcy’s Law for laminar flow in a porous medium. As it is possible that steam
flow might be turbulent, we use the Forcheimer equation for steam pressure
$p$, which combines Darcy flow with the turbulent Ergun equation [29],

$\nabla p=-\frac{\mu_{v}}{k}{\bf u}-\frac{\varrho_{s}c_{F}}{\sqrt{k}}{\bf
u}|{\bf u}|\;,$ (11)

where $c_{F}$ is an order one coefficient, and $\varrho_{s}$ is the density of
steam.

For an ideal gas, $\beta T=1$. Averaging the point-wise energy equation (3)
for steam and the rock energy equation (5) over a representative elementary
volume of the porous magma gives the averaged energy equation for steam and
rock,

$\varrho c\frac{\partial T}{\partial t}+\varrho_{s}c_{ps}{\bf u}\cdot\nabla
T-\phi\frac{\partial p}{\partial t}-{\bf u}\cdot\nabla
p=\kappa_{e}\nabla^{2}T+\overline{\tau:\nabla{\bf u}}$ (12)

where the temperature $T$ in steam is assumed to be equal to that in adjacent
rock, ${\bf u}=\phi{\bf v}$ is the Darcy velocity of steam,
$\kappa_{e}=(1-\phi)\kappa_{m}+\phi\kappa_{s}$ is the effective thermal
conductivity of magma with steam in the pores, and $\kappa_{s}$ is the thermal
conductivity of steam.

The last (drag) term is the viscous dissipation in the steam flow, averaged
over a representative elementary volume. This is the heat generated by
viscosity in shearing flow for steam as it moves through porous magma. This
term is often neglected [26, E.4]. We estimate its size here using [30, eqn
10.7.24] (see also [31]) and Darcy’s Law,

$\overline{\tau:\nabla{\bf u}}=-\frac{\mu}{k}{\bf
u}^{2}=-\frac{k}{\mu}\left(\nabla p\right)^{2}\;.$

Pressure gradient is estimated using a maximum pressure difference that equals
a tensile strength of 1MPa for magma over a lengthscale of 0.1m to give a
magnitude for averaged viscous dissipation of the order of 106. This is a
factor of $10^{3}$ smaller than the sensible heating rate $\varrho
c\frac{\partial T}{\partial t}\sim 10^{9}$, using a timescale of one second,
and a factor of $10^{2}$ smaller (on a lengthscale of one pore) than the
diffusion term, $\kappa_{e}\nabla^{2}T\sim 10^{8}$. So we neglect viscous
dissipative heating due to steam drag in the magma region, compared to heating
due to diffusion from nearby hot magma.

The variation of $\varrho c$ due solely to changes in steam temperature of
1200∘C and changes in steam pressure of 20 bars (typical tensile strength of
volcanic rock) is less than 0.5%. It is dominated by the thermal capacity of
the magma, which typically varies from values near 800 J.kg-1.K-1 at 400K,
rising rapidly to values near 1300 at glass temperatures near 900-1000K [32].
We have assumed a constant value for $\varrho c$, neglecting variations of
$\pm 30$%.

The ideal gas law is used in the vapour region, so that for $r>s(t)$,

$\varrho_{s}=\frac{pM}{RT}\;.$ (13)

###  3.3 Flashing Front

We assume that boiling occurs in a thin moving region located at $r=s(t)$ that
separates the liquid and vapour regions. We acknowledge the spherical symmetry
of our model by taking all variables to depend only on $r$ and $t$. Then
vapour and liquid velocities have only a radial component $u=\phi v$,
$u_{l}=\phi v_{l}$ respectively. We write the energy equations in enthalpy
form, and integrate them and the mass conservation equations with respect to
volume across the moving flashing front, to obtain

$\phi\varrho_{s}h_{sl}({v}-\dot{s})=\phi\varrho_{l}h_{sl}({v}_{l}-\dot{s})=\left[\kappa\nabla
T\right]_{-}^{+}+\phi({v}-{v}_{l})p\;,$ (14)

where $h_{sl}=h_{s}-h_{l}$ is the specific heat of vapourisation, and

$\left[\kappa\nabla T\right]_{-}^{+}=\kappa_{e}\nabla
T(0^{+})-\kappa_{el}\nabla T(0^{-})\;.$

At the flashing front, pressure and temperature are related by the Clausius-
Clapeyron equation,

$p=p^{e}_{0}e^{\frac{Mh_{sl}}{RT^{e}_{0}}\left[\frac{T-T^{e}_{0}}{T}\right]}$
(15)

where $T^{e}_{0}$ and $p^{e}_{0}$ are reference temperature and pressure
values for liquid and vapour phases of water at equilibrium.

These equations (6)–(15) form our dimensional model equations. Boundary
conditions are

$T(R_{2})=300K\;,\quad p(R_{2})=p_{a}\;,\quad\frac{\partial T}{\partial
r}=\frac{\partial p}{\partial r}=0\;\;\;{\rm at}\;r=0\;.$

Temperature and pressure are assumed to be continuous across the flashing
front. Initial conditions are that the temperature of the magma is $T_{m}$ and
the temperature of the inclusion is at boiling point for atmospheric pressure,
$T_{i}$. Initial pressures are taken to be $p_{a}$ everywhere.

###  3.4 Fragmentation Criterion

We use the criterion for fragmentation that steam pressure exceeds the
critical value

$p_{c}=(1-\phi)\sigma_{Y}\;,$ (16)

where $\sigma_{Y}=2$MPa is a typical value for tensile strength of magmatic
rock. This criterion arises out of the use of a soil mechanics approach in
shock tube fragmentation research and presented in [13], which is also based
in part on Biot’s work on wave propagation in fluid filled porous media [33].
The factor $(1-\phi)$ accounts for the fraction of the magma that is load-
bearing. A similar criterion is recommended in other studies of shock tube
fragmentation of volcanic rock, including [34, eqn 11] and [35].

##  4 Model Reduction

We rescale (see also Table 1) using

$r=\tilde{r}R_{2}\;,\quad\;s=\tilde{s}R_{2}\;,\quad\;T=\tilde{T}T_{m}\;,\quad\;p=\tilde{p}p_{a}\;,\quad
t=\tilde{t}t_{0}\;,$ $v_{l}=v_{0l}\tilde{v}_{l}\;,\quad
v=v_{0s}\tilde{v}\;,\quad\varrho_{l}=\varrho_{0l}\tilde{\varrho}_{l}\;,\quad\varrho_{s}=\varrho_{0s}\tilde{\varrho}_{s}\;.$

Noting that flashing of liquid to vapour is the driving force for pressure
change inside the bomb, the time scale $t_{0}$ is chosen to be the time
required to flash all of the liquid in the slurry inclusion to steam. A
balance between the energy required to do this, and the heat provided by
conduction across the slurry surface at $r=R_{1}$ gives the estimate
$t_{0}=\frac{\phi\varrho_{0l}h_{sl}R_{1}^{2}}{3\kappa_{e}(T_{m}-T_{i})}\approx
17$.

We use the ideal gas law to provide $\varrho_{0s}=p_{a}M/(RT_{m})$, and since
the source of steam is the flashing front we balance $\dot{s}\varrho_{l}$ with
$v_{s}\varrho_{s}$ giving

$v_{0s}=\frac{R_{2}\varrho_{0l}}{t_{0}\varrho_{0s}}\;.$ (17)

In the following equations, we drop the tilde on dimensionless variables for
simplicity of notation. The critical pressure for fragmentation is rescaled to
give the nondimensional value

$p_{c}=20(1-\phi)\;,$ (18)

Table 2: Parameters Parameter | Definition | Typical Value | Units
---|---|---|---
$\delta_{1}$ | $\beta_{\ell}\,p_{a}/\varrho_{0l}$ | $4.6\times 10^{-8}$ |
$\delta_{2}$ | $\dfrac{t_{0}v_{0s}\phi\varrho_{0s}c_{ps}}{R_{2}\varrho c}$ | 0.4 |
$\delta_{3}$ | $\dfrac{\phi p_{a}}{\varrho cT_{m}}$ | $1.5\times 10^{-5}$ |
$\delta_{4}$ | $\dfrac{\phi v_{0s}p_{a}t_{0}}{R_{2}\varrho cT_{m}}$ | 0.09 |
$\delta_{5}$ | $\dfrac{\kappa_{e}t_{0}}{\varrho cR_{2}^{2}}$ | $2\times 10^{-3}$ |
$\epsilon_{\varrho}$ | $\varrho_{0s}/\varrho_{0l}$ | $1.7\times 10^{-4}$ |
$\epsilon_{1}$ | $v_{0l}/v_{0s}$ | $2\times 10^{-7}$ |
$\epsilon_{2}$ | $\phi\beta_{\ell}\,p_{a}/(\varrho^{\prime}c^{\prime})$ | $5\times 10^{-12}$ |
$\epsilon_{3}$ | $\kappa_{el}t_{0}/(\varrho^{\prime}c^{\prime}R_{2}^{2})$ | $10^{-3}$ |
$\epsilon_{4}$ | $\dfrac{kp_{a}}{\mu_{v}R_{2}\phi v_{0s}}$ | $2\times 10^{-3}$ |
$\epsilon_{5}$ | $\dfrac{kp_{a}t_{0}}{\phi\mu_{v}R^{2}_{2}}$ | 14 |
$\epsilon_{6}$ | $\dfrac{\sqrt{k\;}\varrho_{0s}c_{F}\phi v_{0s}}{\mu_{v}}$ | $0.08$ |
$\lambda_{1}$ | $\alpha T_{m}/\varrho_{0l}$ | 0.6 |
$\lambda_{2}$ | $p_{a}/(\varrho_{0s}h_{sl})$ | 0.26 |
$\varrho_{0s}$ | $\dfrac{p_{a}M}{RT_{m}}$ | 0.17 | kg.m-3
$\varrho_{0l}$ | reference value of $\varrho_{l}$ | 1000 | kg.m-3
$H$ | $Mh_{sl}/(RT_{0}^{e})$ | 13 |
St | $h_{sl}\phi\varrho_{0s}R_{2}v_{0s}/(T_{m}\kappa_{e})$ | 212 |
$t_{0}$ | $\dfrac{\phi\varrho_{0l}h_{sl}R_{1}^{2}}{3\kappa_{e}(T_{m}-T_{i})}$ | 13 | s
$T_{0}$ | $T_{i}/T_{m}$ | 0.29 |
$v_{0s}$ | $\dfrac{R_{2}\varrho_{0l}}{t_{0}\varrho_{0s}}$ | 46 | m.s-1
$v_{0l}$ | $\dfrac{\alpha\kappa_{el}(T_{c}-T^{w}_{0})}{3R_{1}\varrho c(\varrho_{0l}-\alpha(T_{c}-T^{w}_{0}))}$ | $10^{-5}$ | m.s-1
$v_{dl}$ | $v_{0l}t_{0}/R_{2}$ | $10^{-3}$ |

###  4.1 In the liquid region

In the slurry, the equation of state (8) nondimensionalises to

$\varrho_{l}=1+\delta_{1}(p-1)-\lambda_{1}(T-T_{0})\;.$ (19)

Due to the symmetry of the problem the only mechanisms that can cause liquid
to move in the slurry are liquid density changes. The relatively small size of
$\delta_{1}$ (Table 2) means that thermal expansion is the main factor. We
over-estimate the velocity scale $v_{0l}$ by calculating the change in
dimensional radius due to thermal expansion as the slurry is heated to near
critical temperature $T_{\rm crit}$, and the time to heat the ball to this
temperature by conduction. Ignoring losses of liquid from the slurry, the
slurry mass is the volume times the density, so that if the slurry starts at
radius $R_{1}$, density $\varrho_{0l}$, and temperature $T_{0}^{w}$, and
expands to a new radius $R_{1}+\Delta R_{1}$, then $(R_{1}+\Delta
R_{1})^{3}\left[\varrho_{0l}-\alpha(T_{\rm
crit}-T^{w}_{0})\right]=R_{1}^{3}\varrho_{0l}$. For small $\Delta R_{1}$ this
gives $\Delta R_{1}\approx R_{1}\alpha(T_{\rm
crit}-T^{w}_{0})(3(\varrho_{0l}-\alpha(T_{\rm crit}-T^{w}_{0})))$. Estimating
the time $t_{s}$ required to heat the slurry by conduction gives $t_{s}\approx
R_{1}^{2}\varrho c/\kappa_{el}$, so that the liquid velocity scale is
approximated by

$v_{0l}\approx\frac{\Delta R_{1}}{t_{s}}=\frac{\alpha\kappa_{el}(T_{\rm
crit}-T^{w}_{0})}{3(\varrho_{0l}-\alpha(T_{\rm crit}-T^{w}_{0}))R_{1}\varrho
c}\approx 10^{-5}{\rm m}.{\rm s}^{-1}\;.$

The conservation of mass equation in the slurry becomes, after rescaling and
dropping the tildes,

$\phi\frac{\partial\varrho_{l}}{\partial t}+v_{dl}\frac{\partial}{\partial
r}(\varrho_{l}{u}_{l})=0\;.$ (20)

Because $v_{dl}\sim 10^{-3}\ll 1$ and $u_{l}=0$ at origin, equation (20)
implies that to leading order in $v_{dl}$, $\varrho_{l}\sim 1$, and $u_{l}=0$.
Then pressures in the slurry do not vary appreciably with radius, and can be
taken to be the same as the time-varying pressure value at the flash front. It
follows also that $\varrho^{\prime}c^{\prime}$, which depends mainly on
$\varrho_{l}$, can now be treated as constant.

The conservation of energy equation (9) in the slurry rescales to give, after
dropping the tildes and neglecting terms involving liquid velocity,

$\frac{\partial T}{\partial t}-\epsilon_{2}T\frac{\partial p}{\partial
t}=\frac{\epsilon_{3}}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial T}{\partial r}\right)\;,$ (21)

where the parameters $\epsilon_{2}$, $\epsilon_{3}$ are given in Table 2. The
parameter $\epsilon_{2}\sim 10^{-11}$ is eight orders of magnitude smaller
than $\epsilon_{3}\sim 10^{-3}$, indicating that this pressure-work term may
be neglected. However, we retain for now the term with the small parameter
$\epsilon_{3}$, anticipating that the net temperature gradient in the thermal
boundary layer at the flashing front is what drives changes in pressure. Hence
in the slurry, we reduce to the energy equation

$\frac{\partial T}{\partial
t}=\frac{\epsilon_{3}}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial T}{\partial r}\right)\;,\quad r<s(t)\;.$ (22)

###  4.2 In the vapour region

In the hot magma vapour transport region $r>s(t)$, Forcheimer’s equation (11)
nondimensionalises to $\epsilon_{4}\nabla p=-v-\epsilon_{6}\varrho_{s}v|v|$.
Since $\epsilon_{6}\approx 0.008\ll 1$ (see Table 2), the Ergun term is not
required111Note that the pore Reynolds number
${\text{R}e}_{p}=\rho_{f}\sqrt{k}\phi|u|/\mu_{v}$, with simulations giving
$v\approx 10$m/s, leading to ${\text{R}e}_{p}\approx 0.01$. Hence the flow is
laminar., and Forcheimer’s equation reduces to Darcy’s law,

$v=-\epsilon_{4}\nabla p\;.$ (23)

The mass equation (10) becomes after nondimensionalising and substituting for
$v$,

$\frac{\partial\varrho_{s}}{\partial
t}=\epsilon_{5}\nabla\cdot\left(\varrho_{s}\nabla p\right)\;.$ (24)

The energy equation (12) takes the nondimensional form

$\frac{\partial T}{\partial t}+\delta_{2}\varrho_{s}v\frac{\partial
T}{\partial r}-\delta_{3}\frac{\partial p}{\partial
t}-\delta_{4}v\frac{\partial p}{\partial
r}=\frac{\delta_{5}}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial T}{\partial r}\right)\;,$ (25)

where the meanings and values of the parameters are given in Table 2.

Using Darcy’s law (23) to replace $v$ gives

$\frac{\partial T}{\partial t}-\epsilon_{4}\delta_{2}\varrho_{s}\frac{\partial
p}{\partial r}\frac{\partial T}{\partial r}-\delta_{3}\frac{\partial
p}{\partial t}+\delta_{4}\epsilon_{4}\left(\frac{\partial p}{\partial
r}\right)^{2}=\frac{\delta_{5}}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial T}{\partial r}\right)\;.$ (26)

The most significant terms in the above equation give a diffusion equation for
temperature,

$\frac{\partial T}{\partial
t}=\frac{\delta_{5}}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial T}{\partial r}\right)\;,\quad r>s(t)\;.$ (27)

We retain the diffusion term containing the small parameter $\delta_{5}\sim
1.7\times 10^{-3}$, as we need the temperature gradients in the thermal
boundary layers at the flashing front to calculate the speed of the flashing
front. We have neglected the heat advection term involving
$\epsilon_{4}\delta_{2}\sim 7.4\times 10^{-4}$, which is nearly half the size
of the diffusion term. We have also neglected the pressure-work terms
involving $\delta_{3}$ and $\delta_{4}\epsilon_{4}$, which are 10-2 times
$\delta_{5}$. We have checked a posteriori that they remain relatively small
despite pressure increases of order ten simulated later in this paper.

The ideal gas law takes the nondimensional form

$p=\varrho_{s}T\;.$ (28)

###  4.3 At the flashing front

Conservation of mass and energy across the flashing front $r=s(t)$ as
expressed in equations (14) give

$\dot{s}(\varrho_{l}-\epsilon_{\varrho}\varrho_{s})=-\varrho_{s}v+v_{dl}\varrho_{l}v_{l}$
(29)

and

$\varrho_{s}(v-\epsilon_{\varrho}\dot{s})=\frac{1}{{\rm St}}\left[\nabla
T\right]_{-}^{+}+\lambda_{2}(v-\epsilon_{1}v_{l})p\;,$ (30)

where constants are defined in Table 2. We have neglected the small relative
difference between $\kappa_{e}=2$ and $\kappa_{el}=3$ to obtain equation (30).
The Stefan number St$\sim$200 is relatively large, consistent with the
assumption of a narrow two-phase region at the flashing front.

Noting the small size of $\epsilon_{\varrho}\sim 10^{-4}$ and
$\epsilon_{1}\sim 10^{-7}$, we drop the advective mass and heat transport
terms in equations (29) and (30), and use Darcy’s law (23) to obtain

$\dot{s}=\epsilon_{4}\varrho_{s}\frac{\partial p}{\partial r}=-\frac{1}{{\rm
St}}\left[\frac{\partial T}{\partial
r}\right]_{-}^{+}+\lambda_{2}\epsilon_{4}p\frac{\partial p}{\partial r}\;.$

The parameter combination $\lambda_{2}\epsilon_{4}\sim 5\times 10^{-4}$ is ten
times smaller than 1/St$\sim 5\times 10^{-3}$, so we discard the pressure-work
term involving $\lambda_{2}\epsilon_{4}$ to give the reduced jump conditions

$\dot{s}=\epsilon_{4}\varrho_{s}\frac{\partial p}{\partial r}=-\frac{1}{{\rm
St}}\left[\frac{\partial T}{\partial r}\right]_{-}^{+}\;.$ (31)

The Clausius-Clapeyron condition becomes, in dimensionless terms, if we choose
the reference values to be $p_{0}^{e}=p_{a}$ and $T_{0}^{e}=T_{i}=373K$,

$p=\exp\left[H\left(\frac{T-T_{0}}{T}\right)\right]\;.$ (32)

This condition, giving the pressure dependence of the vapourization
temperature of water, is only valid for the two-phase conditions that apply at
the flashing front $r=s(t)$, and require that $T\in[T_{0},0.5]$ approximately.

Our reduced system then consists of equation (22) for temperature diffusion in
the slurry, where pressures and densities are spatially constant, equation
(24) for nonlinear pressure diffusion (coupled with temperature) in the
surrounding vapour region, equation (27) for temperature diffusion in the
vapour region, the ideal gas law (28) relating density, pressure and
temperature in the vapour region, with boundary conditions (31) and (32) at
the moving flashing front $r=s(t)$ providing the temperature, pressure and
speed $\dot{s}$ there, and boundary conditions at origin and at the surface of
the bomb

$p=1\;,\;\;r=1\;;\quad\frac{\partial T}{\partial r}=0\;,\;\;r=0\;;\quad
T=T_{0}\;,\;\;r=1\;.$

The temperature gradients at the flashing front provide a vapour flux boundary
condition, dominated by the gradient on the hot magma side that drives
pressure up at the flashing front and forces vapour outwards into the bomb.
Typical initial conditions would be

$T=T_{0}\;,\;\;r<s(0)\;;\quad T=1\;,\;\;r>s(0)\;;\quad p=1\;;\quad
s(0)=R_{1}/R_{2}\;,$

with a step change in temperature at the flashing front.

###  4.4 Thermal Boundary Layers

In the magma and in the slurry, thermal diffusivity $\sim 0.002$ is small
compared to the nonlinear diffusivity of pressure given by $\epsilon_{5}p\sim
70$. So we note that temperature changes are expected to propagate more slowly
than pressure changes, and we seek approximations to the temperatures in the
thermal boundary layers near the flashing front, which drive steam production
via equation (31), which specifies the flux at the moving flash boundary for
equation (24).

In the magma, we consider an inner region described by a radial coordinate
$\sigma$ given by $r=\epsilon+\sqrt{\delta_{5}\;}\;\sigma$, that is close to
the flashing front that starts at $\epsilon=R_{2}/R_{1}$, at times that are
early enough to ignore movement of the flashing front. Then the temperature
equation in the magma in this inner region becomes

$T_{t}=\frac{1}{(\epsilon+\sqrt{\delta_{5}}\sigma)^{2}}\frac{\partial}{\partial\sigma}\left[(\epsilon+\sqrt{\delta_{5}}\sigma)^{2}\frac{\partial
T}{\partial\sigma}\right]\;.$

Considering the limit $\delta_{5}\to 0$, and taking the inner solution valid
for $\sqrt{\delta_{5}\;}\sigma\ll\epsilon$, we have the boundary layer
equation (with subscripts $t$ and $\sigma$ indicating partial derivatives)

$T_{t}=T_{\sigma\sigma}\;,$ (33)

with boundary conditions $T(0,t)=T_{f}$ and $T\to 1$ as $\sigma\to\infty$. The
temperature at the flashing front $T_{f}\sim 0.4$ varies with time, but only
by about 5%, being limited by the critical temperature. The outer solution is
$T=1$ in the rest of the magma, away from the boundary layer. We take
advantage of the relatively small variation in $T_{f}$, by taking it to be a
constant value.

The boundary layer equation (33) then admits a similarity solution.
Considering the similarity variable $\eta=\sigma^{2}/t$, the partial
differential equation (33) becomes an ordinary differential equation
$4T_{\eta\eta}+\frac{2}{\eta}T_{\eta}+T_{\eta}=0$. This is first order in
$T_{\eta}\equiv dT/d\eta$, and may be integrated twice to obtain the solution

$T=(1-T_{f})\;\;{\rm erf}\left(\frac{\sigma}{2\sqrt{t}}\right)+T_{f}\;,$ (34)

where ${\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-u^{2}}du$, and the
constants have been chosen to match the flash temperature $T_{f}$ when
$\sigma=0$ and the outer solution $T=1$ as $t\to 0$. This inner solution
provides the value of the temperature gradient at the flashing front, on the
magma side, as

$\frac{\partial T}{\partial
r}=\frac{1-T_{f}}{\sqrt{\pi\delta_{5}t\;}}\approx\frac{0.6}{\sqrt{\pi\delta_{5}t\;}}\;.$
(35)

The same approach with $\sigma_{2}=\frac{\epsilon-r}{\sqrt{\epsilon_{3}}}$
provides the temperature in the boundary layer in the slurry,

$T=T_{f}+(T_{0}-T_{f})\;\;{\rm
erf}\left(\frac{\sigma_{2}}{2\sqrt{t}}\right)\;.$ (36)

Then the temperature gradient on the slurry side of the flashing front is
estimated as

$\frac{\partial T}{\partial
r}=\frac{T_{f}-T_{0}}{\sqrt{\pi\epsilon_{3}t\;}}\approx\frac{0.04}{\sqrt{\pi\epsilon_{3}t\;}}\;.$
(37)

Noting that $\epsilon_{3}\approx\delta_{5}$, we see that the contribution of
the temperature gradient in the slurry is about one-fifteenth that of the
temperature gradient in the magma, due to the smaller temperature differences
between flash temperature and slurry temperature, and to the similar thermal
diffusivities.

##  5 Numerical Solutions

(a) Pressures in the hot magma versus time and $\xi=(r-s)/(1-s)$, where $r$ is
the radius.

(b) The location of the flashing front $s(t)$ versus time.

(c) Simulated temperatures (symbols) and the similarity solution (solid line)
in the hot magma vs the similarity variable $\eta$.

(d) Simulated temperatures (symbols) and the similarity solution (solid line)
in the slurry vs the similarity variable $\eta$.

(e) Pressure at the flashing front versus time

(f) Maximum simulated pressures versus the base ten log of permeability
(symbols).

Figure 4: Simulation of flashing to steam inside a Surtseyan bomb. All
variables are nondimensional. Physical constant values are as listed in Table
(1), except in figure (f), where permeability values are taken to depend on
porosity as does the fit (1) to measured values for Surtseyan bombs. Eleven
porosity values were chosen evenly spaced in the range [0.1,0.95], providing
the indicated values of permeability $k$. The subhorizontal dashed line in (f)
is the critical pressure (18) above which fragmentation occurs. The vertical
dashed line indicates the smallest permeability measured in our data from
intact bombs.

Our reduced system is solved using the method of lines, and coded in Matlab.
The moving boundary at the flashing front is fixed in place using Landau
transformations, which lead to advective terms depending on the speed of the
boundary, in the partial differential equations describing the reduced model.
The thermal boundary layers are resolved by transforming the spatial variable
adjacent to the flashing front, which is now fixed in place, in both the
slurry and hot magma regions, so as to obtain greater resolution at the
flashing front. Spatial derivatives are then replaced by equispaced
differences in the new spatial variable. The resulting system of coupled
ordinary differential equations is stiff. We found consistent results using
1200 spatial mesh points and the stiff solver ode15s in Matlab. Examples of
the code are available at
https://doi.org/10.5281/zenodo.5214799.

Typical results are plotted in Fig. (4), from simulations run using the
parameter values listed in Table (1). Pressures rise rapidly to a global
maximum at the flashing front, diffusing out into the hot magma region, then
slowly decaying as the liquid in the slurry boils away to nothing.

The smallest mesh size used in these simulations corresponds to ten microns
adjacent to the flashing front, a typical value for pore size in Surtseyan
bombs and pumice. Initial temperatures ramp from boiling point in the slurry
to magma temperatures over one pore size. Plot (f) in Fig. (4) shows the
maximum pressure reached at the flashing front versus permeability. The
initial temperature profile used for this plot matches the similarity
solution, ramping from flash temperature to close to magma values over a
distance of ten microns, with a smallest computational mesh size corresponding
to one micron.

In the time taken for pressures to reach their maximum value, the flashing
front has moved less than $10^{-4}$ times its original position. Nevertheless,
turning off the advective terms in the conservation equations solved
numerically, gives a much closer match than that shown in Fig. (4) to the
similarity solutions for temperatures in slurry and in magma, which were
derived by ignoring movement of the flashing front. The multiple temperatures
plotted at $\eta=0$ in Fig. (4), especially noticeable in the slurry
temperature plot, correspond to the changing boiling temperature at the
flashing front at early times, due to pressure changes there.

(a)

(b)

Figure 5: Maximum dimensionless pressure reached (a), and the log of the
dimensionless time taken to reach it (b), versus the log of the dimensionless
size of the smallest mesh used in computer simulations of the reduced model.
Logs are to the base ten. Initial temperatures are a step function,
effectively a ramp over the smallest mesh size. Parameter values used are as
listed in Table (1).

Pressure rises rapidly at the flashing front, due to the influx of steam
there, driven mainly by the temperature gradient on the hot magma side. The
initial temperature profile is at present theoretically a step function, with
a gradient that is in theory infinite. In the computer code this step change
is approximated by a ramp with a steepness that depends on the smallest mesh
size. The results plotted in Fig. (5) show how the maximum pressure at the
flashing front and the time taken to reach it, vary with the steepness of the
initial temperature profile. The largest of these sizes used in Fig. (5)
corresponds to a dimensional value of a few microns. Ten microns would be a
typical pore size, and might be taken to be the smallest mesh size in terms of
the present application. It is clear that the maximum pressure increases as
the ramp tends towards a step, and that the time taken decreases roughly
proportionately to the distance over which the ramp operates.

If the initial temperature profile ramps over a (dimensional) distance that is
fixed at a representative pore size of 10$\mu$, simulated maximum pessures
converge rapidly as the smallest mesh reduces through 1$\mu$, so that the
pressures we present in Fig. (4) are accurate to within 0.1 bar. We note that
the size of one pore is less than the size of a representative elementary
volume, so that we are technically exploring beyond the limits of the
continuum approach used here.

The dimensionless time range within which maximum pressure is realised is 10-8
to 10-6, corresponding to a dimensional range of about 10-7 to 10-5s,
indicating consistency with the assumption that viscous flow effects can be
neglected, since viscous flow effects take more than a second to be
significant.

These results raise mathematical questions, irrespective of any geophysical
(textural) indications of a practical smallest mesh size or representative
elementary volume, about what determines the time taken to reach the maximum
pressure at the flashing front, and whether the maximum pressure value seen in
Fig. (5) is approaching some limit as the ramp approaches a step (as smallest
mesh size approaches zero). We now derive an upper bound on the maximum
pressure. We find that our bound diverges as the initial temperature gradient
diverges, suggesting that the pressure maximum may be theoretically unbounded
in this limit.

##  6 Maximum Pressure

We seek an approximation to the maximum pressure, in the form of a formula
that relates it to the material properties of the magma and the enclosed
slurry. Such a formula was obtained previously [12] by using the steady-state
pressure behaviour to obtain an upper bound, and was possible only because the
temperature field was approximated as constant. The present model makes it
clear that temperature gradient changes are significant, and drive boiling, so
it is important to properly account for the effect of early high temperature
gradients at the flashing front on the maximum pressure developed by heating.

###  6.1 Upper Bound

We begin by deriving an upper bound on the maximum pressure developed at early
times at the flashing front. We consider the reduced equation (24) for steam
pressure or density, driven by the gradient of the temperature on the magma
side of flash. We ignore the relatively small contribution from the
temperature gradient on the slurry side, and movement of the flashing front is
ignored.

Initial temperature in the magma is taken to be a similarity solution that
starts at an earlier time $t=-t_{e}<0$, so that eqn (37) is modified to read

$\frac{\partial T}{\partial
r}=\frac{T_{f}-T_{0}}{\sqrt{\pi\epsilon_{3}(t+t_{e})\;}}\;,\;\;r=\epsilon\;.$
(38)

This is the flux from an initial temperature profile that is ramped, and a
step function initial profile is recovered in the limit $t_{e}\to 0$ . Noting
that the error function is close to one in value when its argument is two, the
distance $\Delta r_{e}$ over which the initial temperature ramps from the
flash value to the value one is given by $\Delta r_{e}\approx
4\sqrt{\delta_{5}t_{e}\;}$.

Early pressure behaviour is governed by equation (24), which close to the
flashing front can be written in the form

$\frac{\partial}{\partial
t}\left(\frac{p}{T}\right)=\epsilon_{5}\frac{\partial}{\partial
r}\left(\frac{p}{T}\frac{\partial p}{\partial r}\right)\;.$ (39)

At early times, pressure changes due to the influx of vapour at the flashing
front propagate a distance $\Delta r\approx\sqrt{\epsilon_{5}t\;}$ into the
magma, which together with the spatially constant initial pressures in the
magma suggests the approximation at the flashing front

$\frac{\partial}{\partial r}\left(\frac{p}{T}\frac{\partial p}{\partial
r}\right)\approx-\frac{p}{T\Delta r}\frac{\partial p}{\partial r}\;.$

We put these together and take advantage of the relatively slow time rate of
change of temperature to write equation (39) as

$\frac{\partial p}{\partial
t}\approx-\sqrt{\frac{\epsilon_{5}}{t}\;}\;p\frac{\partial p}{\partial r}\;.$
(40)

At the flashing front, we can then approximate the early time pressure
behaviour by using the reduced jump conditions (31) and ignoring heat flow
into the slurry to set

$p\frac{\partial p}{\partial r}\approx-\frac{T(1-T)}{\epsilon_{4}{\rm
St}\sqrt{\pi\delta_{5}(t+t_{e})\;}}\;.$ (41)

Substituting this into equation (40) gives the following approximation for
early pressure changes at the flash point:

$\frac{\partial p}{\partial t}\approx\frac{B_{1}}{\sqrt{t^{2}+tt_{e}\;}}\;,$
(42)

where $B_{1}=\frac{T(1-T)}{\epsilon_{4}{\rm
St}}\sqrt{\frac{\epsilon_{5}}{\pi\delta_{5}}\;}$ and $T\approx 0.4$. The
solution is

$p^{e}\approx
p_{0}+B_{1}\left[\ln\left(2\sqrt{t^{2}+t_{e}t\;}+2t+t_{e}\right)-\ln(t_{e})\right]\;.$
(43)

The timescale $t^{*}$ for reaching a maximum flash pressure is estimated by
calculating when $p^{e}$ crosses the pressure null surface at flash, where
$\partial p/\partial t=0$. We estimate this by considering the density null
surface, where $\partial\varrho/\partial t=0$, since temperatures remain of
order one. Then equation (24) becomes the quasi-steady state equation

$\frac{\partial}{\partial r}\left(r^{2}\frac{p}{T}\frac{\partial p}{\partial
r}\right)=0\;,\quad r>s(t)\;,$ (44)

with solution

$p_{q}^{2}=2C_{1}\int\frac{T}{r^{2}}dr\;.$

Since $T\leqslant 1$, $p_{q}$ has an upper bound given by
$\left(p_{q}^{u}\right)^{2}=2C_{1}\int r^{-2}dr=-2C_{1}/r+C_{2}$. The boundary
conditions $p(1)=1$ and the flux condition (41) give $C_{2}=1+2C_{1}$,
$C_{1}=-\epsilon^{2}(1-T)/(\epsilon_{4}{\rm
St}\sqrt{\pi\delta_{5}(t+t_{e})\;})$. At $r=\epsilon$, this upper bound is
closely approximated by

$p_{q}^{u}=\frac{B_{2}}{\left(t+t_{e}\right)^{\frac{1}{4}}}\;,\;\;B_{2}^{2}=\frac{2\epsilon(1-\epsilon)(1-T)}{\epsilon_{4}{\rm
St}\sqrt{\pi\delta_{5}\;}}$ (45)

Dropping the term $p_{0}$ in equation (43) and equating just the leading terms
in equations (43) and (45) as $t\to 0$, the time $t^{*}$ at which flash
pressure reaches a maximum is estimated as the solution to
$2B_{1}\sqrt{t^{*}/t_{e}}=B_{2}/((t_{e})^{\frac{1}{4}})$, which is
$t^{*}=B_{2}^{2}\sqrt{t_{e}\;}/(4B_{1}^{2})$. The maximum pressure at the
flashing front has an upper bound estimate of
$p_{\text{max}}=B_{2}t_{e}^{-1/4}$, and replacing $t_{e}$ with its equivalent
lengthscale $\Delta r_{e}$ gives
$p_{\text{max}}^{2}=8\epsilon(1-\epsilon)(1-T_{f})/(\epsilon_{4}{\rm
St}\sqrt{\pi}\Delta r_{e})$. where $T_{f}\approx 0.4$ is the nondimensional
flash temperature.

This can be rewritten as

$p_{\text{max}}^{2}=\frac{8\epsilon(1-\epsilon)(1-T_{f})T_{m}\mu_{v}\kappa_{e}}{\sqrt{\pi\;}\Delta
r_{e}\,k\;p_{a}h_{sl}\rho_{0s}}\;,$ (46)

providing an estimate of an upper bound on the maximum pressure that develops
at the flashing front, and indicating its dependence on key properties of the
magma bomb.

This estimate is compared to simulated maximum pressures in Fig. (6), and we
see the estimated pressures are about three times the simulated values. So the
estimate provides a weak upper bound for simulated maximum pressures.

(a) Simulated maximum pressure (red line with symbols), theoretical upper
bound on pressure (46) (solid black line), and a representative value
20(1-$\phi$) (dashed black line) for the tensile strength of magmatic rock
when cooled, versus permeability.

(b) The log of the theoretical upper bound versus the log of the simulated
maximum pressure (red symbols), together with a line indicating where values
are equal.

Figure 6:  Log-log plots of simulated maximum pressures and theoretical upper
bounds. Pressures are dimensionless, but are effectively in bars. Permeability
values (m2) depend on porosity as in the fit (1) to measured values for
Surtseyan bombs, with eleven porosity values chosen evenly spaced in the range
[0.2,0.8], providing the indicated values of permeability k. Units of
permeability are m2. All logarithms are to the base ten.

###  6.2 Fragmentation Criterion

While the upper bound illustrated in Fig. (6) is correct, and it is in the
form of a formula, it is too weak to provide a useful fragmentation criterion
for Surtseyan bombs. The log-log graph of the upper bound in Fig. (6)
indicates that it would be improved as an estimate of the numerical solutions
by simply shifting it in log-log space. This suggests trying a simple linear
approximation over the relevant range of porosities and permeabilities, to
relate the upper bound more closely to the more accurate numerical results.
The average of numerical maximum pressure values is compared to the average of
upper bound values over the range illustrated, and gives the following
approximation $P^{\tt est}_{\text{max}}$ to simulated pressures,

$P^{\tt est}_{\text{max}}=0.29p_{\text{max}}$ (47)

where $p_{\text{max}}$ is given by eqn 46.

Figure 7: Log-log plot of simulated maximum nondimensional pressure (line with
red symbols), and the simple fitted theoretical estimate $P^{\tt
est}_{\text{max}}$ (Eqn. (47) (black line), for various permeabilities with
porosities matching our data. The dashed subhorizontal black line is a
representative value 20(1-$\phi$) bars for the tensile strength of the
magmatic rock when cooled. The blue vertical dashed line indicates the
smallest permeability measured for intact samples of Surtseyan ejecta. Units
of permeability are m2. Logs are to the base ten. Other parameter values are
as listed in Table (1).

This approximation $P^{\tt est}_{\text{max}}$ is compared to numerical values
in Fig. (7), and to the smallest permeability $k=0.3$ darcys observed in data
from intact Surtseyan bombs. Fragmentation of a bomb corresponds to a maximum
pressure that exceeds the tensile strength approximated by the dashed line in
Fig. (7). The numerical simulations (red circles) predict fragmentation when
$k<0.6$ darcys, while the approximate formula $P^{\tt est}_{\text{max}}$
(solid black line) predicts fragmentation when $k<3$ darcys. The data in Fig.
(3) indicate variability in measured permeability, with the smallest measured
value for permeability of intact bombs in the confidence range 0.1–1.0 darcys.
This contains the value $k=0.6$ darcys, that is, the numerically predicted
smallest permeability value is inside the confidence range for smallest
measured permeability of intact bombs. This provides some support from field
measurements for the theoretical modelling described here. There is
significant uncertainty in the assumed value of tensile strength, which also
affects the minimum permeability predictions for intact bombs. Note that no
model fitting to data has been conducted here.

##  7 Discussion

The estimate $P^{\tt est}_{\text{max}}$ is a simple formula that approximates
the numerical values of maximum pressure in a Surtseyan bomb. This provides us
with a fragmentation criterion, that a Surtseyan bomb should fragment if the
nondimensional pressure exceeds the nondimensional tensile strength
20(1-$\phi$), that is, combining eqns (47) and (46), if

$P^{\tt
est}_{\text{max}}=\sqrt{\frac{0.7\epsilon(1-\epsilon)(1-T_{f})T_{m}\mu_{v}\kappa_{e}}{\sqrt{\pi\;}\Delta
r_{e}\,k\;p_{a}h_{sl}\rho_{0s}}}>20(1-\phi)$ (48)

We can compare our fragmentation criterion (48) with the dimensional estimate
for maximum pressure $p^{\tt
old}_{\text{max}}=\sqrt{\frac{7(1-\epsilon)(T_{m}-T_{0})\mu_{v}\kappa_{e}}{\,k\;h_{sl}}\left(\frac{RT_{m}}{M}\right)}$
obtained in earlier work [12] by writing our result in equation (48) in terms
of the dimensional value for pressure

$p^{\tt
est}_{\text{max}}=\sqrt{\frac{0.7(1-\epsilon)(1-T_{f})T_{m}\mu_{v}\kappa_{e}}{\sqrt{\pi\;}\,k\;h_{sl}}\left(\frac{RT_{m}}{M}\right)\left(\frac{\epsilon}{\Delta
r}\right)\;}\;\;\;{\rm Pa}\approx\sqrt{\frac{0.1\epsilon}{\Delta
r}\;}\;\;p^{\tt old}_{\text{max}}\;.$ (49)

The form of $p^{\tt old}_{\text{max}}$ differs from $p^{\tt est}_{\text{max}}$
most significantly in that our new result has an extra factor
$\approx\sqrt{0.1\epsilon/\Delta r\;}$, involving the ratio of the inclusion
size to the ramping distance for the assumed initial temperature profile. This
extra factor means that our present estimate of minimum permeability of 0.6
darcys for a bomb to survive in flight is significantly higher than the
estimate of 0.02 darcys in [12]. This previous estimate is well below the
smallest permeability (0.1–1.0 darcys) measured in intact bombs, while our new
result lies inside this range. In this sense, our present model is more
consistent with the data from intact bombs, since the previous result says
that some intact bombs would be expected to have permeabilities below the
measured smallest values 0.1–1.0 darcys.

The term $\Delta r$ is a measure of the size of the region over which initial
temperature ramps from boiling values at the flashing front up to magma
values. In the limit as it goes to zero, and the initial temperature profile
approaches a step function, our upper bound tends to infinity. Our approximate
early time solution (43) also grows without bound in this limit, suggesting
our model may be mathematically ill-posed if it is assumed that the initial
temperature profile is a step function. Such an assumption is common for the
initial temperature profiles in heat conduction problems, where thanks to the
stable behaviour of solutions to diffusion equations, it does not cause any
problems. The ill-posedness might be regularised by re-scaling pressure to
include previously neglected small parameters like the pressure-work term in
the $\dot{s}$ equation just before equation (31), or by recognising the small
distance over which initial flash pressures change and hence including viscous
dissipation terms.

However, a step function initial temperature profile is not physically
realistic, and a ramped profile already provides us with a mathematically
better behaved model. In a volcanological context, the length-scale $\Delta r$
might be considered to be the size of a representative elementary volume,
which would be several times the mean pore size.

The measurements made on intact bombs shown in Fig (3) suggest minimum
permeabilities in the range 0.1–1.0 darcys. This is broadly consistent with
our model results. We note that this is a one-sided view of all ejected bombs,
applying only to intact bombs. While it is consistent with the hypothesis that
bombs with smaller permeability will fragment, there are no measurements
available for fragmented bombs to confirm this. Even if there were such
measurements, we would have no assurance of the cause of fragmentation.
Fragmentation due to impact with the ground is not modelled here — we consider
only the mechanism of steam generation and the associated steam pressures.
Given the highly vesicular nature of the ejected magma, it is possible that
all ejected bombs have higher permeabilities than the critical value for
fragmentation. If so, our model then provides a physical explanation for why
typical Surtseyan ejecta are not expected to fragment due to steam pressure
buildup from flashing of slurry inclusions after ejection, in that the vapour
typically escapes before pressures reach fragmentation values.

##  8 Conclusion

We have developed and simplified a fully transient sphere-within-a-sphere
model for the pressure increase expected to occur inside Surtseyan bombs after
ejection, due to flashing to steam of liquid in slurry inclusions. We reduced
our results to a single formula (49) for the maximum pressure developed at the
flashing front, revealing how that pressure depends on permeability and
relative size of the inclusion. This formula, when compared to the tensile
strength of the bomb, provides a new fragmentation criterion for Surstseyan
bombs.

Our model neglects the effect of rotation, which we expect to assist steam
transport and reduce the maximum pressure if included in the model. The
thermal capacity of magma has been approximated as a constant value, ignoring
its variation of $\pm$30% over the temperature range modelled. The geometry of
bombs is not spherical, and multiple inclusions are found inside each one at
varying distances from the surface of the bomb. The deformation of the porous
matrix near slurry inclusions that is observed in samples suggests that a more
sophisticated approach like that in [13] might be useful, explicitly
considering stress and strain in the rock matrix. Our work already shows that
pressure changes due to flashing will propagate more rapidly than temperature
changes, relevant to questions about compression effects near inclusions. A
region of reduced permeability in the matrix next to an inclusion may also be
a useful future modification of our model.

The mechanism for insertion of slurry into the hot magma immediately prior to
ejection is important for initial temperature profiles, and is a more
complicated flow and heat transport problem that may lead to estimates of
$\Delta r$ that better reflect the volcanology of Surtseyan eruptions.

Our results highlight the importance of the initial temperature gradient in
the hot magma adjacent to the flashing front, driving boiling and hence
pressure increase there. Numerical simulations illustrated in Fig. (4) using a
ramping distance for the initial temperature that is given by a typical pore
size indicate that fragmentation is expected to occur for any bombs with
permeabilities less that about one darcy. For these relatively low
permeability bombs, steam pressure build-up due to heating by surrounding hot
magma is not adequately relieved by steam escape through the porous magma.
This, together with the observation that all measured permeabilities of intact
ejecta are greater than one darcy, provides an explanation for why most
Surtseyan bombs survive steam pressures developed inside due to flashing of
slurry inclusions.

EG and MM contributed equally to the mathematical modelling, asymptotics, and
numerical simulations; IS provided all data, the Field Work and Data section,
and most of the volcanological expertise, together with the fundamental
fragmentation question that motivated this work.

We are grateful to Andrew Fowler (University of Limerick and University of
Oxford) for fruitful discussions during the early development of our model.

## References

  * [1] Thorarinsson S. 1967 Surtsey: The New Island in the North Atlantic. Viking, New York.
  * [2] Kokelaar B. 1983 The mechanism of Surtseyan volcanism. Jour. Geological Soc. 140, 939–944.
  * [3] Kokelaar P. 1986 Magma-water interactions in subaqueous and emergent basaltic volcanism. Bull Volcanol 48, 275–289.
  * [4] Schipper CI, White JD. 2016 Magma-slurry interaction in Surtseyan eruptions. Geology 44, 1–7.
  * [5] Németh K, Kósik S. 2020 Review of Explosive Hydrovolcanism. Geosciences 44, 1–27.
  * [6] Schipper CI, White JD, Jakobsson S, Palin J, Bush-Marcinowski T. 2015 The Surtsey Magma Series. Scientific Reports 5, 1–17.
  * [7] Büttner R, Zimanowski B, Mohrholz CO, Kummel R. 2005 Analysis of thermohydraulic explosion energetics. Journal of Applied Physics 98.
  * [8] Zimanowski B, Büttner R, Lorenz V. 1997 Premixing of magma and water in MFCI experiments. Bulletin of Volcanology 58, 491–495.
  * [9] Schipper CI, Whilte JD, Zimanowski B, Büttner R, Sonder I, Schmid A. 2011 Experimental interaction of magma and "dirty" coolants. Earth Planet Science Letters 303, 232–336.
  * [10] Dingwell DB. 1996 Volcanic Dilemma: Flow or Blow?. Science 273, 1054–1055.
  * [11] Saar MO, Manga M. 1999 Permeability-porosity relationship in vesicular basalts. Geophysical Research Letters 26, 111–114.
  * [12] McGuinness MJ, Greenbank E, Schipper CI. 2016 Modelling vapour transport in Surtseyan bombs. Journal of Volcanology and Geothermal Research 318, 103–113.
  * [13] Fowler AC, Scheu B, Lee WT, McGuinness MJ. 2010 A theoretical model of the explosive fragmentation of vesicular magma. Proceedings of the Royal Society London A 466, 731–752.
  * [14] Babb T, Benham G, Bows J, Gonzalez-Farina R, Kiradjiev K, Lee W, Tibos S. 2021 Predicting lift-off time when deep-frying potato dough snacks. SIAM Journal on Applied Mathematics 81, 574–590.
  * [15] Alidibirov M, Dingwell DB. 2000 Three fragmentation mechanisms for highly viscous magma under rapid decompression. Journal of Volcanology and Geothermal Research 100, 413–421.
  * [16] Wadsworth FB, Witcher T, Vasseur J, Dingwell DB, Scheu B. 2017 When does magma break?. Advances in Volcanology DOI: 10.1007/11157_2017_23.
  * [17] SchipperI CI, Voyer ML, Moussallam Y, Whilte JD, Thordarson T, J-I JIK, Chang Q. 2016 Degassing and magma mixing during the eruption of Surtsey Volcano (Iceland, 1963-1967): the signatures of a dynamic and discrete rift propagation event. Bulletin of Volcanology 78, 1–19.
  * [18] McGuinness M, Greenbank E, Schipper CI. 2021 "A theoretical model of Surtseyan bomb fragmentation" code and data. .
  * [19] Walker G, Croasdale R. 1972 Characteristics of some basaltic pyroclasts. Bulletin of Volcanology 35, 303–317.
  * [20] Thomas R, Sparks R. 1992 Cooling of tephra during fallout from eruption columns. Bulletin of Volcanology 54, 542–553.
  * [21] Schneider CA, Rasband WS, Eliceiri KW. 2012 NIH Image to ImageJ: 25 years of image analysis. Nature Methods 9, 671–675.
  * [22] Degruyter W, Burgisser A, Bachmann O, Malaspinas O. 2010 Synchotron X-ray microtomography and lattice Boltzmann simulations of gas flow through volcanic pumices. Geosphere 6, 470–481.
  * [23] Rust AC, Cashman KV. 2004 Permeability of vesicular silicic magma: inertial and hysteresis effects. Earth Planet Science Letters 228, 93–107.
  * [24] Ishii M, Hibiki T. 2006 Thermo-fluid dynamics of two-phase flow. NY: Springer.
  * [25] Batchelor G. 1967 An Introduction to Fluid Dynamics. Cambridge University Press.
  * [26] Fowler AC. 2011 Mathematical Geoscience. Springer.
  * [27] Sato H, Uematsu M, Watanabe K, Saul A, Wagner W. 1988 New international skeleton tables for the thermodynamic properties of ordinary water substance. Journal of Physical Chemistry Reference Data 17, 1439–1540.
  * [28] Robertson EC. 1988 Thermal properties of rocks. Open-file report 88-441 US Geological Survey Reston, Virginia, USA.
  * [29] Nield D, Bejan A. 2006 Convection in porous media. Springer New York.
  * [30] Bear J. 1988 Dynamics of Fluids in Porous Media. Dover, NY.
  * [31] Pilotti M, Succi S, Menduni G. 2002 Energy dissipation and permeability in porous media. Europhysics Letters 60, 72–78.
  * [32] Nicholls A, Potuzak M, Dingwell D. 2009 Cooling rates of basaltic hyaloclastites and pillow lava glasses from the HSDP2 drill core. Geochimica et Cosmochimica Acta 73, 1052–1066.
  * [33] Biot MM. 1956 Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range. The Journal of the Acoustical Society of America 28, 168–178.
  * [34] Koyaguchi T, Mitani N. 2005 A theoretical model for fragmentation of viscous bubbly magmas in shock tubes. Journal of Geophysical Research 110, 1–21.
  * [35] Koyaguchi T, Scheu B, Mitani N, Melnik O. 2008 A fragmentation criterion for highly viscous bubbly magmas estimated from shock tube experiments. J. Volc. & Geoth. Res. 178, 58–71.


Suppose $L=K\times K$, and let $\varpi_{1}=(\varpi_{K},1)$ and
$\varpi_{2}=(1,\varpi_{K})$ be uniformizers of the two maximal ideals of
$\mathcal{O}_{L}$. Then

$\sum_{n_{1},n_{2}\geq
0}W^{\mathrm{new}}\left(\left(\begin{smallmatrix}\varpi_{1}^{n_{1}}\varpi_{2}^{n_{2}}&0\\\
0&1\end{smallmatrix}\right)\right)q_{K}^{(1/2-s)(n_{1}+n_{2})}=L^{\mathrm{std}}(\pi_{1},s)L^{\mathrm{std}}(\pi_{2},s).$

###  8.2. The Asai zeta integral

####  8.2.1. Definitions

For $W\in\mathcal{W}(\pi)$ and $\Phi\in\mathcal{S}(K^{2})$, we define

$Z(W,\Phi;s)\coloneqq\int_{Z_{H}N_{H}\backslash
H}W(h)f^{\Phi}(h;\chi,s)\,\mathrm{d}h$

where $\chi=\omega_{\pi}|_{K^{\times}}$. This is well-defined for $\Re(s)\gg
0$ and has analytic continuation as a meromorphic function of $q_{K}^{s}$; as
usual, we define $L^{\operatorname{As}}(\pi,s)$ as the unique $L$-factor which
generates the fractional ideal of $\mathbb{C}[q_{K}^{\pm s}]$ given by the
values of this integral as $W$ and $\Phi$ vary.

####  8.2.2. Torus integral

We define for $W\in\mathcal{W}(\pi)$ a “partial zeta integral”

$Z(W;s)=\int_{K^{\times}}W\left(\left(\begin{smallmatrix}x&0\\\
0&1\end{smallmatrix}\right)\right)|x|_{K}^{s-1}\,\mathrm{d}^{\times}x.$

Cf. [Gro21, Definition 3.3.1]. Again, the integral is convergent for
$\Re(s)\gg 0$ and has analytic continuation as a meromorphic function of
$q_{K}^{s}$.

######  Proposition 8.2.1.

The function $Z(hW;\chi,s)$ transforms as an element of the induced
representation $I(|\cdot|^{1/2-s},|\cdot|^{s-1/2}\chi)$, i.e. the dual of the
representation in which the Siegel section lies; and we have

$\int_{B_{H}\backslash H}Z(hW;s)f^{\Phi}(h;\chi,s)\,\mathrm{d}h=Z(W,\Phi,s).$

###### Proof.

This follows by splitting up the integral over $Z_{H}N_{H}\backslash H$ into
an integral over $Z_{H}N_{H}\backslash B_{H}$ followed by an integral over
$B_{H}\backslash H$. ∎

Note that it suffices to take the integral over
$K_{H}=\operatorname{GL}_{2}(\mathcal{O}_{K})$, since $B_{H}K_{H}=H$.

####  8.2.3. A special Schwartz function

Let $r\geq 1$, and consider the Schwartz function

$\Phi_{r}=(q_{K}^{2}-1)q_{K}^{2r-2}\cdot\operatorname{ch}(\mathfrak{p}_{K}^{r},1+\mathfrak{p}_{K}^{r}).$

From [LSZ22, Lemma 3.2.5], if $\chi$ is trivial on $1+\mathfrak{p}^{r}$, then
$f^{\Phi_{r}}(-;\chi,s)$ vanishes outside the coset $B_{H}\backslash
B_{H}K_{0}(\mathfrak{p}_{K}^{r})$ (the “small cell” modulo
$\mathfrak{p}^{r}$), and we have

$f^{\Phi_{r}}(1;\chi,s)=(q_{K}+1)q_{K}^{r-1},$

which is also the volume of the small cell for the unramified Haar measure on
$B_{H}\backslash H$. (On the other hand, if $\chi$ is nontrivial on
$1+\mathfrak{p}^{r}$ then $f^{\Phi_{r}}(-;\chi,s)$ is identically zero.)

Since the small cell modulo $\mathfrak{p}^{r}$ is the orbit of the identity
under the subgroup $\left(\begin{smallmatrix}1&0\\\
\mathfrak{p}^{r}&1+\mathfrak{p}^{r}\end{smallmatrix}\right)$, we have the
following: if $W\in\mathcal{W}(\pi)$ is invariant under this subgroup, the
integrand $h\mapsto Z(hW;s)f^{\Phi}(h;\chi,s)$ is equal to
$f^{\Phi_{r}}(1)Z(W,s)$ on the small cell in $(B_{H}\cap K_{H})\backslash
K_{H}$, and vanishes everywhere else on $K_{H}$. So for any $W$ invariant
under this subgroup we have

$Z(W,\Phi_{r},s)=Z(W;s).$

###  8.3. Imprimitive Asai $L$-factors

######  Definition 8.3.1.

We define the _imprimitive Asai $L$-factor_ by

$L^{\mathrm{imp}}_{\operatorname{As}}(\pi,s)=\begin{cases}Z(W^{\mathrm{new}},s)&\text{if
$\pi$ is ramified}\\\ L(\chi,2s)\cdot Z(W^{\mathrm{new}},s)&\text{if $\pi$ is
unramified}.\end{cases}$

(Note that “unramified” here means “has invariants under
$\operatorname{GL}_{2}(\mathcal{O}_{L})$”, even when $L$ is a ramified
extension.)

######  Proposition 8.3.2.

Let $r_{0}$ be such that $\varpi_{K}^{r_{0}}=\mathfrak{c}\cap\mathcal{O}_{K}$,
where $\mathfrak{c}\trianglelefteqslant\mathcal{O}_{L}$ is the conductor of
$\pi$. Then we have

$L^{\mathrm{imp}}_{\operatorname{As}}(\pi,s)=Z(W^{\mathrm{new}},\Phi_{r_{0}},s).$

If $r_{0}\geq 1$ then we have in fact
$L^{\mathrm{imp}}_{\operatorname{As}}(\pi,s)=Z(W^{\mathrm{new}},\Phi_{r},s)$
for any $r\geq r_{0}$.

####  8.3.1. Unramified case

Suppose $\pi$ has trivial conductor (but we allow $L$ to be a ramified
extension). Then the new vector is simply the normalised spherical Whittaker
function $W^{\operatorname{sph}}$; and the function $h\mapsto
Z(hW^{\operatorname{sph}},s)$ is invariant under $K_{H}$, so we want to take
$r=0$. In this case $f^{\Phi_{r}}(1;\chi,s)=L(\chi,s)$, so by a well-known
explicit computation we find

$L_{\mathrm{As}}^{\mathrm{imp}}(\pi,s)Z(W^{\operatorname{sph}},\Phi_{0},s)=L(\chi,s)Z(W^{\operatorname{sph}},s)=L_{\mathrm{As}}(\pi,s).$

######  Remark 8.3.3.

This computation is well-known if $L=K\times K$ or $L$ is an unramified
quadratic extension, in which case $L_{\mathrm{As}}(\pi,s)$ is a degree 4
$L$-factor. However the computation is also valid for $L$ a ramified field
extension, in which case $L_{\mathrm{As}}$ is a degree 3 $L$-factor (and
coincidentally agrees with the symmetric square $L$-factor of $\pi$).

####  8.3.2. Explicit formulae

We give explicit formulae for the imprimitive Asai $L$-factor defined above.

######  Proposition 8.3.4.

Suppose $L$ is a field extension (ramified or unramified), and $\pi$ has non-
trivial conductor. Let $C$ be the multiset of nonzero complex numbers such
that

$L^{\mathrm{std}}(\pi,s)=\prod_{\alpha\in C}(1-\alpha q_{L}^{-s})^{-1}.$

Then

$L^{\mathrm{imp}}_{\mathrm{As}}(\pi,s)=\prod_{\alpha\in
C}(1-\alpha^{e_{L/K}}q_{K}^{-s})^{-1}$

where $e_{L/K}$ is the ramification index (1 or 2).∎

######  Proposition 8.3.5.

Suppose $L=K\times K$ (so again $q_{K}=q_{L}\eqqcolon q$), and
$\pi=\pi_{1}\times\pi_{2}$ is arbitrary (ramified or unramified). For $i=1,2$
let $C_{i}$ be the multiset of complex numbers such that
$L^{\mathrm{std}}(\pi_{i},s)=\prod_{\alpha_{i}\in
C_{i}}(1-\alpha_{i}q^{-s})^{-1}$. Then we have

$L^{\mathrm{imp}}_{\mathrm{As}}(\pi,s)=\prod_{\alpha_{1}\in
C_{1}}\prod_{\alpha_{2}\in C_{2}}(1-\alpha_{1}\alpha_{2}q^{-s})^{-1}.$

####  8.3.3. Relation to the primitive Asai $L$-factor

In all cases there exists an $L$-factor $P(s)$ such that

$L_{\mathrm{As}}(\pi,s)=P(s)\cdot L^{\mathrm{imp}}_{\mathrm{As}}(\pi,s).$

We have seen that if $\pi$ is unramified, then $P(s)=1$. One checks similarly
that if $L=K\times K$, then we have $P(s)=1$ in any of the following cases:

  * •

at least one of $\pi_{1}$, $\pi_{2}$ is unramified;

  * •

one of the $\pi_{i}$ is supercuspidal, but the other is not.

###  8.4. Zeta integrals with $U^{t}$-eigenvectors

We now consider the local zeta integral at the split finite place $v=p$. Since
we want to consider twists by ramified characters, we no longer assume the
character $\nu$ is trivial.

We can identify
$\mathcal{W}(\Pi_{1},\psi_{p})\otimes\mathcal{W}(\Pi_{2},\psi_{p})$ with
$\mathcal{W}(\Pi_{p},\psi_{F,p})$ via mapping $W_{1}\otimes W_{2}$ to the
function $g\mapsto W_{1}(\left(\begin{smallmatrix}1/\sqrt{D}&\\\
&1\end{smallmatrix}\right)g_{1})W_{2}(\left(\begin{smallmatrix}-1/\sqrt{D}&\\\
&1\end{smallmatrix}\right)g_{2})$ (so normalised new vectors go to normalised
new vectors). A straightforward substitution gives

$Z(W,\Phi;\nu,s)=\hat{\nu}_{p}^{-1}(-\sqrt{D})\cdot\int_{\mathbb{Q}_{p}^{\times}}W_{1}(\left(\begin{smallmatrix}-1&\\\
&1\end{smallmatrix}\right)g)W_{2}(g)\hat{\nu}_{p}^{-1}(\det
g)f^{\Phi}(g;\chi\nu^{-2},s)\,\mathrm{d}g,$

and the integral on the right-hand side is a Rankin–Selberg zeta integral. We
are interested in the following situation:

  * •

$\pi_{1}$ and $\pi_{2}$ are unramified;

  * •

$\nu$ is a Dirichlet character of $p$-power conductor (possibly trivial), so
$\hat{\nu}_{p}^{-1}(-\sqrt{D})=\nu(-\sqrt{D})$;

  * •

$\Phi=\Phi_{p,\nu^{-1},\nu}$ as defined above;

  * •

$W_{2}$ is the normalised spherical vector;

  * •

$W_{1}$ is the $U_{\mathfrak{p}_{1}}^{t}$-eigenvector given by
$(1-\tfrac{\beta_{1}}{U_{\mathfrak{p}_{1}}^{t}})W_{1}^{\operatorname{sph}}$,
where $W_{1}^{\operatorname{sph}}$ is the normalised spherical vector.

######  Proposition 8.4.1.

With the above choices, we have

$Z(W,\Phi;\nu,s)=\tfrac{p}{p+1}\nu(-\sqrt{D})\mathcal{E}_{p}(\operatorname{As}\Pi,\nu^{-1},s)L_{\operatorname{As}}(\Pi_{p}\times\hat{\nu}_{p}^{-1},s).$

This is a rather standard computation (compare e.g. [CH20]) so we shall not
give the proof in detail here.

###  8.5. Archimedean case

It remains to compute the Archimedean zeta integral, for which we use the
formula of 6.1.2 for $W_{\infty}$. Since $\hat{\nu}_{\infty}$ is trivial on
the support of $W_{\infty}$, we need to compute

$Z(W_{\infty},{\Phi_{\infty}}^{(k_{1}-k_{2})},s)\coloneqq\int_{Z(\mathbb{R})N(\mathbb{R})\backslash\operatorname{GL}_{2}(\mathbb{R})}W_{\infty}(h)f^{\Phi_{\infty}}(h;1,s)\,\mathrm{d}h.$

As shown in [Jac72, p.85], since the integrand is invariant under
$\mathbb{R}^{\times}\times{\rm SO}_{2}(\mathbb{R})$, we can write this as

$f^{\Phi_{\infty}^{(k_{1}-k_{2})}}(1;1,s)\cdot\int_{\mathbb{R}^{\times}_{+}}y^{(k_{1}+k_{2})/2}e^{-4\pi
y/\sqrt{D}}y^{s-1}\,\mathrm{d}^{\times}y.$

We have

$\int_{\mathbb{R}^{\times}_{+}}y^{(k_{1}+k_{2})/2}e^{-4\pi
y/\sqrt{D}}y^{s-1}\,\mathrm{d}^{\times}y=\tfrac{1}{2}\cdot\left(\tfrac{\sqrt{D}}{2}\right)^{s-1+\tfrac{k_{1}+k_{2}}{2}}\cdot\Gamma_{\mathbb{C}}\left(s-1+\tfrac{k_{1}+k_{2}}{2}\right),$

where $\Gamma_{\mathbb{C}}(s)\coloneqq 2(2\pi)^{-s}\Gamma(s)$, and similarly

$f^{\Phi_{\infty}^{(k_{1}-k_{2})}}(1;1,s)=i^{k_{1}-k_{2}}2^{s-\tfrac{k_{1}-k_{2}}{2}}\Gamma_{\mathbb{C}}\left(s+\tfrac{k_{1}-k_{2}}{2}\right),$

so we conclude that

$Z(W_{\infty},{\Phi_{\infty}}^{(k_{1}-k_{2})},s)=i^{k_{1}-k_{2}}2^{-k_{1}}(\sqrt{D})^{s-1+(k_{1}+k_{2})/2}\cdot\Gamma_{\mathbb{C}}(s-1+\tfrac{k_{1}+k_{2}}{2})\Gamma_{\mathbb{C}}(s+\tfrac{k_{1}-k_{2}}{2}).$

###  8.6. The global zeta integral

We can finally put together the zeta integrals of the above sections to prove
our main result.

######  Theorem 8.6.1.

Let $\Pi$ be an automorphic representation as above, with central character
$\chi_{\Pi}$. Let $\mathcal{L}_{\Pi}$ the one-variable measure defined in
Definition 7.3.2. Then for any locally-algebraic $s+\nu$ of
$\mathbb{Z}_{p}^{\times}$, with $s$ an integer such that
$\tfrac{1-k_{1}+k_{2}}{2}\leq s\leq\tfrac{k_{1}-k_{2}}{2}$ and $\nu$ of finite
order, we have

$\mathcal{L}_{p,\operatorname{As}}(\Pi)(s+\nu)=\mathcal{E}_{p}(\operatorname{As}(\Pi),\nu^{-1},s)\cdot\frac{\Gamma(s+\tfrac{k_{1}+k_{2}-2}{2})\Gamma(s+\tfrac{k_{1}-k_{2}}{2})}{2^{(k_{1}-2)}i^{(1-k_{2})}(-2\pi
i)^{(2s+k_{1}-1)}}\cdot\frac{L_{\operatorname{As}}\left(\Pi,\nu^{-1},s\right)}{\Omega(\Pi)}.$

The same holds with $\mathcal{L}_{p,\operatorname{As}}^{\mathrm{imp}}(\Pi)$
and $L_{\operatorname{As}}^{\mathrm{imp}}(\Pi,\nu^{-1},s)$ in place of the
primitive $L$-functions.

###### Proof.

We first consider the primitive $L$-function. By definition and the discussion
in section 7.3, we have

$\mathcal{L}_{p,\operatorname{As}}(\Pi)(s+\nu)=(\star)\cdot\prod_{v}\int_{Z_{H}(\mathbb{Q}_{v})N_{H}(\mathbb{Q}_{v})\backslash
H(\mathbb{Q}_{v})}\nu^{-1}(\det(h_{v}))W_{v}(h_{v})f^{\Phi_{v}}(h_{v};\chi_{v}\nu^{-2},s)\,\mathrm{d}h_{v}.$

where $(\star)$ is the term from (7.5). We have chosen our test data so that
the terms in this product for all places $v\notin\\{p,\infty\\}$ are the local
$L$-factors $L_{\operatorname{As}}(\Pi_{v}\times\nu_{v}^{-1},s)$. From 8.4.1,
the local term at $p$ is given by

$\tfrac{p}{p+1}\nu(-\sqrt{D})\cdot
L_{\operatorname{As}}(\Pi_{p}\times\nu_{p}^{-1},s)\cdot\mathcal{E}_{p}(\operatorname{As}\Pi,\nu^{-1},s).$

The local integral at $v=\infty$ was computed in the previous section, and
combining these gives the stated result.

In the case of the imprimitive $p$-adic $L$-function, we have seen above that
our choice of local data at the ramified primes (the “imprimitive test data”
defined in section 7.3) gives the imprimitive Asai $L$-factor, and the
computation at $p$ and $\infty$ is the same as before. ∎

### Acknowledgements

This paper was written while all three authors were in residence at the
Mathematical Sciences Research Institute in Berkeley, California, for the
program “Algebraic Cycles, $L$-values and Euler Systems” during the spring of
2023 (supported by NSF grant no. DMS-1928930). It is a pleasure to thank MSRI
and its staff for providing such a splendid working environment.

## References

  * [BP22] G. Boxer and V. Pilloni, _Higher Hida and Coleman theories on the modular curve_, Épijournal Géom. Algéb. 6 (2022), epiga:6112. MR 4482376.
  * [CH20] S.-Y. Chen and M.-L. Hsieh, _On primitive $p$-adic Rankin-Selberg $L$-functions_, Development of Iwasawa theory—the centennial of K. Iwasawa’s birth, Adv. Stud. Pure Math., vol. 86, Math. Soc. Japan, Tokyo, 2020, pp. 195–242. MR 4385083.
  * [CPR89] J. Coates and B. Perrin-Riou, _On $p$-adic $L$-functions attached to motives over ${\bf Q}$_, Algebraic number theory (in honour of K. Iwasawa), Adv. Stud. Pure Math., vol. 17, Academic Press, Boston, MA, 1989, pp. 23–54. MR 1097608.
  * [Far11] L. Fargues, _La filtration canonique des points de torsion des groupes $p$-divisibles_, Ann. Sci. Éc. Norm. Supér. (4) 44 (2011), no. 6, 905–961, with collaboration of Yichao Tian. MR 2919687.
  * [Gro20] G. Grossi, _On norm relations for Asai–Flach classes_, Int. J. Number Theory 16 (2020), no. 10, 2311–2377.
  * [Gro21] by same author, _Higher Hida theory for Hilbert modular varieties in the totally split case_, preprint, 2021, arXiv:2106.05666.
  * [GLZ23] G. Grossi, D. Loeffler, and S. Zerbes, _Asai–Flach classes and $p$-adic L-functions_, in preparation, 2023.
  * [Har90] M. Harris, _Period invariants of Hilbert modular forms. I: Trilinear differential operators and L-functions_, Cohomology of arithmetic groups and automorphic forms (Luminy, 1989), Lect. Notes Math., vol. 1447, Springer, 1990, pp. 155–202. MR 1082965.
  * [HT01] M. Harris and J. Tilouine, _p-adic measures and square roots of special values of triple product L-functions_, Math. Ann. 320 (2001), no. 1, 127–147.
  * [HK92] M. Harris and S. Kudla, _Arithmetic automorphic forms for the nonholomorphic discrete series of $\operatorname{GSp}(2)$_, Duke Math. J. 66 (1992), no. 1, 59–121. MR 1159432.
  * [Hid85] H. Hida, _A $p$-adic measure attached to the zeta functions associated with two elliptic modular forms. I_, Invent. Math. 79 (1985), no. 1, 159–195. MR 774534.
  * [Jac72] H. Jacquet, _Automorphic forms on ${\rm GL}(2)$. Part II_, Lecture Notes in Mathematics, vol. 278, Springer, 1972. MR 0562503.
  * [Kat78] N. Katz, _$p$ -adic $L$-functions for CM fields_, Invent. Math. 49 (1978), 199–297.
  * [KM76] F. Knudsen and D. Mumford, _The projectivity of the moduli space of stable curves. I: Preliminaries on ”det” and ”Div”_, Math. Scand. 39 (1976), 19–55.
  * [LLZ18] A. Lei, D. Loeffler, and S. L. Zerbes, _Euler systems for Hilbert modular surfaces_, Forum Math. Sigma 6 (2018), e23, 67. MR 3876763.
  * [LPSZ21] D. Loeffler, V. Pilloni, C. Skinner, and S. Zerbes, _Higher Hida theory and $p$-adic $L$-functions for ${GSp}_{4}$_, Duke Math. J. 170 (2021), no. 18, 4033–4121. MR 4348233.
  * [LSZ22] D. Loeffler, C. Skinner, and S. L. Zerbes, _Euler systems for $\operatorname{GSp}(4)$_, J. Eur. Math. Soc. 24 (2022), no. 2, 669–733. MR 4382481.
  * [LW20] D. Loeffler and C. Williams, _$p$ -adic Asai $L$-functions of Bianchi modular forms_, Algebra Number Theory 14 (2020), no. 7, 1669–1710. MR 4150247.
  * [Oh23] G. Oh, _Higher Hida theory and the p-adic L-functions for $U(2,1)\times U(1,1)$_, preprint, 2023.
  * [Pan82] A. Panchishkin, _Le prolongement $p$-adique analytique des fonctions $L$ de Rankin. I_, C. R. Acad. Sci. Paris Sér. I Math. 295 (1982), no. 2, 51–53. MR 676360.
  * [Pan94] by same author, _Motives over totally real fields and $p$-adic $L$-functions_, Ann. Inst. Fourier (Grenoble) 44 (1994), no. 4, 989–1023. MR 1306547.
  * [Pil20] V. Pilloni, _Higher coherent cohomology and $p$-adic modular forms of singular weights_, Duke Math. J. 169 (2020), no. 9, 1647–1807. MR 4105535.
  * [Shi81] G. Shimura, _On certain zeta functions attached to two Hilbert modular forms. II. The case of automorphic forms on a quaternion algebra_, Ann. of Math. (2) 114 (1981), no. 3, 569–607. MR 634429.
  * [Urb14] E. Urban, _Nearly overconvergent modular forms_, Iwasawa theory 2012: state of the art and recent advances, Springer, 2014, pp. 401–441.