## Abstract

Recent experimental work has shown that, when a vertical column of rock under large pressure is suddenly depressurized, the column can ‘explode’ in a structured and repeatable way. The observations show that a sequence of horizontal fractures forms from the top down, and the resulting blocks are lifted off and ejected. The blocks can suffer secondary internal fractures. This experiment provides a framework for understanding the way in which catastrophic explosion can occur, and is motivated by the corresponding phenomenon of magmatic explosion during Vulcanian eruptions. We build a theoretical model to describe these results, and show that it is capable of describing both the primary sequence of fracturing and the secondary intrablock fracturing. The model allows us to suggest a practical criterion for when such explosions occur: firstly, the initial confining pressure must exceed the yield stress of the rock, and, secondly, the diffusion of the gas by porous flow must be sufficiently slow that a large excess pore pressure is built up. This will be the case if the rock permeability is small enough.

## 1. Introduction

Depending on magma composition and evolution, explosive volcanic eruptions manifest themselves in a broad variety of activity ranging from mild Strombolian eruptions and Hawaiian fire fountaining to vigorous Vulcanian and Plinian eruptions. The range of different types of explosive eruption is mirrored in the variety of their trigger mechanisms. Vulcanian-style eruptions are typically caused by the sudden unplugging of a sealed volcanic vent. The consequent depressurization of the underlying magma generates stresses in the vesicular magma and may initiate further volatile exsolution; if the magma is sufficiently viscous, the gases are not able to escape from the pressurized vesicles, thus adding further stress to the magma which can cause brittle rupture of the magma itself, resulting in an explosive eruption.

It is this fragmentation of magma which distinguishes explosive volcanic eruptions from the more quiet effusive form of volcanic activity. Hence, the conditions leading to magma fragmentation as well as the fragmentation process itself are the key points for the understanding of the dynamics of volcanic eruptions. So far, various models have been proposed for fragmentation, based on the investigations of the deposits of explosive eruptions, theoretical models and laboratory experiments. Early models mentioned magma disruption by bubble coalescence as the main mechanism (Verhoogen 1951). However, McBirney & Murase (1970) and Sparks (1978) demonstrated that coalescence is unlikely to be the leading fragmentation mechanism for highly viscous magmas, although possibly applicable for eruptions of low-viscosity magmas. Basically, two main groups of processes leading to fragmentation can be discriminated.

### (a) Fragmentation owing to rapidly accelerating two-phase flow

In this case, the driving force for expansion and acceleration is derived from vesiculation and bubble growth. High strain rates within the magma cause the fragmentation, owing to either instabilities in the fluid (ductile) or brittle fracture. Low-viscosity magma is likely to fragment by fluid instabilities resulting in Hawaiian fire fountains or Strombolian bubble bursts. Fragmentation by brittle fracture owing to high strain rates may be important for Plinian eruptions of more viscous magma (Dingwell & Webb 1989; Gilbert & Sparks 1998; Papale et al. 1998; Papale 1999; Cashman et al. 2000).

### (b) Fragmentation of vesicular magma owing to rapid decompression

Here, gas overpressure builds up in the vesicles of the magma, mainly depending on melt properties such as viscosity, volatile content and diffusivity as well as the overburden pressure. Rapid decompression of the vesicular magma can be triggered by the sudden removal of an edifice (e.g. dome collapse or landslide) or by the ejection of a plug from the conduit. A pressure gradient is built up at the top of the magma body between the pressurized gas in the vesicles and the low-pressure area behind the unloading wave. At a certain pressure differential, the tensile strength of the magma is overcome and brittle disruption of the upper layer within the vesicular magma will occur. Now, a newly formed free surface is exposed to low pressures and again a pressure gradient is built up, leading to the fragmentation of the next layer. Thus, a layer-by-layer fragmentation process (a fragmentation wave) moves downwards into the magma (Alidibirov & Dingwell 1996, 2000; Scheu *et al.*2008*a*,*b*). Calculations demonstrate that the fragmentation wave velocity should be far slower than the speed of sound. This was confirmed subsequently experimentally, with values ranging from 2 to 110 m s^{−1} (Kennedy et al. 2005; Scheu *et al.* 2006). This fragmentation process plausibly accounts for the eruption style of most silicic events, e.g. for lateral blasts and Vulcanian to sub-Plinian eruptions. It may further contribute to Plinian eruptions of silicic magma, most probably in combination with the other mechanism mentioned above. Layer-by-layer fragmentation is nowadays widely accepted as being the predominant fracturing process caused by rapid decompression (Cashman *et al.* 2000; Melnik 2000; Ichihara *et al.* 2002; Namiki & Manga 2005; Scheu *et al.*2006, 2008*a*,*b*).

Several physical models have been proposed accounting for the stress distribution caused by pressurized gas within the vesicles of a magma body. McBirney & Murase (1970) used elasticity theory and the Griffith theory of fracture. More recent models are based on the stress distribution of thin-walled spheres (Alidibirov 1994) or thick-walled spheres (Zhang 1999). Based on laboratory experiments, Spieler *et al.* (2004) linked overpressure and porosity in an empirical fragmentation criterion. Comparative studies of these models, including comparisons with experimental datasets, showed that none of the approaches fit the entire range of data accurately enough, leading to the suggestion that the process is intrinsically more complex (Spieler *et al.* 2004; Koyaguchi *et al.* 2008). Koyaguchi *et al.* (2008) added the Griffith theory of crack propagation to the models of Zhang (1999) and Koyaguchi & Mitani (2005), and thus were able to reproduce satisfactorily not only the threshold for magma fragmentation but also its propagation speed (Scheu *et al.* 2006).

Our aim in this paper is to describe a mathematical model based on the first principles of mass and momentum conservation which can explain the experimental results. In principle, this will provide a theoretical framework for the subsequent investigation of the more realistic situation which occurs in an actual volcanic vent. The experimental set-up is described in the following section.

## 2. Experimental investigations

### (a) Experimental set-up

To investigate the fragmentation of porous magma induced by rapid decompression, Alidibirov (1994) designed a vertical shock tube apparatus. This was improved and modified several times to suit different scientific problems (Spieler *et al.* 2004; Scheu *et al.*2006, 2008*a*,*b*), yet the basic principle remains the same: a cylindrical sample (in this study 12 mm radius and 60 mm length) of a porous volcanic rock is glued into a sample holder and placed in an autoclave. A set of diaphragms separates the autoclave from a large chamber at atmospheric pressure. The autoclave with the sample inside is slowly pressurized with nitrogen gas. After an equilibration time, rapid decompression is triggered by a systematic failure of the diaphragm system. When the rarefaction wave reaches the sample, a pressure gradient starts to build up in the sample leading to disruption of the sample and ejection of the fragments (above the fragmentation threshold) or to permeable degassing of the entire sample (below the fragmentation threshold).

Figure 1 depicts the specific set-up used for this study, which was optimized for visual observation (Scheu *et al.* 2008*a*). The autoclave as well as the sample holder was manufactured out of Plexiglas, and all experiments were performed at ambient temperature. Special care was taken to use a transparent glue neither obscuring the view nor changing the appearance of the sample. Three pressure sensors monitor the gas pressure evolution within the autoclave (above the sample), at the side of the sample as well as directly below it. A high-speed video camera (Photron Fastcam) records the experiment at a rate of 10 000 frames per second.

In these experiments, magma is represented by a sample of porous volcanic rock. The use of a rock to represent magma can be justified by the rapidity of the decompression. On the millisecond time scale of these experiments, magma behaves like a brittle solid rather than a viscous liquid (Dingwell 1996, 2006). When the autoclave is decompressed, the sample remains pressurized for a certain time owing to the gas in its open pore space, causing it to fragment.

### (b) Experimental results

In a series of experiments, different glues were tested to ensure optimal experimental conditions. The chosen glue, a special superglue, does not penetrate into the sample and is applied to the entire inner surface of the sample holder. The glue is brittle and only just strong enough to hold the sample in place against the pressure difference developed between its top and bottom when the autoclave is decompressed. If no glue is used, the entire sample is propelled upwards by compressed gas ejected from its base but shows a similar pattern of fractures. If the glue holds the sample in place too strongly, the sample fractures internally leaving a shell of rock attached to the sleeve. In both these cases, it is more difficult to observe the fracture pattern developing.

Figure 2 shows the snapshots of a fragmenting sample recorded by a high-speed camera at 10 000 frames per second. These indicate that the fragmentation starts at the top of the sample and continues downwards. Two types of fractures are observed. The first type is the fracture of the glued column, and usually occurs within the upper third of the sample. The fractures are parallel to the sample surface and dissect the entire sample into discs (layers). Crystals, large pores and pre-existing small cracks may cause a slight deviation from this fracture pattern. The second type is the internal fracturing of an expelled fragment. After a disc is ejected, further fracturing may be observed during ejection. This consecutive fracturing often leads to an entire disaggregation of the sample into ash and small blocky fragments.

A more complicated picture of the second type of fracture is revealed by figure 3, which gives the locations of the upper and the lower surfaces of cracks as a function of time. Slices of rock are ejected from the upper part of the remaining sample, but these slices can break into two or more discs, exhibiting layered fractures while in flight before they shatter into small pieces. This fracturing behaviour was systematically observed in the higher pressure experiments, and cannot be explained by the existing models of layer-by-layer fragmentation previously mentioned. Thus, a comprehensive mathematical model of this process must be able to account for both these types of fracture.

## 3. A mathematical model of fragmentation

We consider the column of rock to be an elastic porous medium, and as such the basic theory for the deformation of the medium is that given by Biot (1956, 1962). The gas can flow through the porous rock, and we suppose the rock fractures if the gas pore pressure exceeds the stress in the solid matrix by a yield stress *σ*_{Y}. When the pressure is lowered at the surface of the sample, the solid stress is lowered very rapidly, while a decompression wave propagates downwards through the pore space. It is the resulting excess pore pressure which can cause fracturing of the sample. In our theoretical discussion, we will describe the way in which sequential primary fracturing can occur, and also how secondary fracture within separated blocks can occur. A distinction that we will make is between experiments where the sample is glued to the container and those where it is not. In the latter, we generally find (as is observed) that the column lifts off the base initially, before subsequently fracturing. We interpret this as being due to a primary fracture at the base (where the yield stress is zero), followed by a subsequent secondary fracture. If the sample is glued, then we suppose that a wall friction acts on the sample owing to the glue, and this allows the possibility for primary internal fracturing. In the model which follows, we describe the elastic displacement of the rock and the turbulent porous flow of the (adiabatically compressible) gas through the rock, allowing for the poroelastic effect of the gas pressure on the rock deformation.

In Biot’s (1956) theory, the solid displacement is denoted **u**, while the displacement of the pore fluid is denoted **U**. The solid strain tensor is denoted *e*_{ij}, so that
3.1
and the dilatations of solid and fluid are defined by
3.2
Biot’s equations for the stresses are then
3.31
and
3.32
where *N*, *A*, *Q* and *R* are the four elastic constants; *A* and *N* are equivalent to Lamé constants for the solid matrix, and *Q* and *R* are associated with the deformability of the pore space and the pore fluid. Note that, by eliminating *ε*, we can write equation (3.3)_{1} in the form
3.4
where
3.5
(Biot writes the stress in the solid in the form .)

For small displacements, the velocity of the pore fluid is defined by
3.6
Biot then poses the conservation of momentum equations for the pore fluid and the matrix, which can be written in the form
3.71
and
3.72
where **A** and **D** are the added mass and the interfacial drag terms, taken by Biot in the form
3.8
The subscripts *t* denote time derivatives, here and subsequently. *ρ*_{a} is an added mass coefficient having units of density, and *b* is an interfacial drag coefficient. We discuss suitable forms for these coefficients below.

The model is completed by the addition of the conservation of mass equations
3.9
for the pore fluid, and
3.10
for the solid, together with an equation of state for the pore fluid, which in the present case we take to be the perfect gas law
3.11
and an energy equation for the gas, which we take to be the adiabatic equation
3.12
where *c*_{p} is the specific heat at constant pressure.

The provenance of the conservation of the mass of solid is a curiosity, because in the elastic description of the solid matrix, mass conservation is not generally considered. In fact, the elastic description implicitly assumes small strains and essentially constant porosity, and we can interpret equation (3.10) as simply stating that *ϕ* is approximately constant. We shall thus not consider this equation further.

We now consider the specification of the coefficients *ρ*_{a} and *b*. The added mass coefficient commonly takes the form
3.13
in dispersed two-phase flows (Robinson *et al.* 2008), where *ϕ* is the volume fraction of the dispersed phase, *ρ*_{c} is the density of the continuous phase, and *C*_{VM} is an *O*(1) number. An example would be the bubbly flow of gas bubbles in a liquid phase, where *ρ*_{c} is the density of the liquid. The added mass force arises through the existence of relative motion between the phases. Because the liquid has to get round the gas bubbles, its corresponding deformation implies that *ρ*_{a}∼*ρ*_{c}. In the present case, the relative motion of solid and pore gas is accommodated by the deformation of the gas phase (both phases are continuous, i.e. connected), and therefore we suppose that *ρ*_{a} is given by
3.14
with *C*_{VM}=*O*(1).

The drag coefficient *b* is most simply determined through analogy with Darcy’s Law. If we ignore the (usually small) acceleration terms in equation (3.7)_{2}, then (bearing in mind that the Darcy flux is *ϕ*(**v**−**u**_{t})) we see that Darcy’s Law is obtained by choosing
3.15
(Biot 1956), where *k* is the permeability. In the present situation, the pore gas flow is sufficiently rapid that the Ergun equation may be appropriate. One form of this which caters for both Darcy (laminar) and Ergun (turbulent) forms is the Forchheimer equation
3.16
(Nield & Bejan 2006), where in the present context **V**=*ϕ*(**v**−**u**_{t}), and *c*_{F} is an *O*(1) coefficient. The second term is dominant if the pore Reynolds number is large, i.e.
3.17
and we anticipate that this is the case. It is then appropriate to define, equivalently to equation (3.16),
3.18

### (a) One-dimensional model

As the experimental results are in essence one-dimensional, we suppose that the variables depend only on the vertical coordinate *z*, where the initial sample lies in 0<*z*<*l*. The vertical solid displacement is denoted *w*, the vertical component of the solid stress is denoted *σ* and the vertical velocity of the pore gas is denoted *v*. A complication arises depending on whether the sample is glued to, or otherwise pressed against, the sample holder. In this case, it is appropriate to prescribe a wall friction *F*. Strictly, this renders a one-dimensional model inexact, although, for long thin samples, cross-sectional averaging will still lead to a one-dimensional model. Here, we simply allow for side-wall friction by the addition of a friction term to the solid momentum equation. The equations describing the flow then take the form
3.19
where
3.20
and *ρ*_{0} and *p*_{0} are the reference values of gas density and pressure. As before, subscripts *t* denote time derivatives, while subscripts *z* denote space derivatives. We suppose the wall friction term is given by
3.21
where *μ*_{g} is the shear modulus of the glue, *d*_{g} is the glue thickness, *R* is the sample radius, and *w*_{0}(*z*) is the initial displacement of the sample (this may be non-zero because we take as the reference state the unstressed sample before it is pressurized).

The model is augmented by the initial conditions 3.22 Suitable boundary conditions are to prescribe solid and fluid stress at the surface, and zero displacement at the base, thus 3.23

The pressure in the chamber at the surface of the sample, *p*_{c}, is a function of time which decreases from its initial value *p*_{0}. The dynamics of this decrease depend on the rate of escape of gas from the roof of the chamber. In fitting the experimental data, we use an exponential of the form
3.24
More general forms could be adopted, but equation (3.24) points out the importance of the chamber relaxation time scale *t*_{c}, which has an important part to play in the dynamics. Just as with a champagne cork, if the pressure decreases sufficiently slowly, the gas will be able to escape without fracturing the rock.

To describe the mechanism of fracture, we follow precepts first put forward in the context of soil mechanics. The total stress on the rock is the weighted sum of the stresses on rock grains and pore space, thus 3.25 and the effective stress is 3.26

We hypothesize that fracturing occurs if the effective stress exceeds a yield stress (1−*ϕ*)*σ*_{Y}, weighted by the solid volume fraction (1−*ϕ*) to indicate the dependence of yield on the number density of intergranular bonds. Fracturing of the sample thus occurs if
3.27
and in this case the boundary conditions in the fractures will be modified, as described later.

### (b) Non-dimensionalization

We scale the variables as follows:
3.28
so that in dimensionless form (using the same letters to denote the dimensionless variables^{1}), we have the dimensionless model
3.29
where the dimensionless parameters are
3.30
*α* was defined in equation (3.5), and is presumed to be *O*(1).

To estimate the values of these parameters, we use the values in table 1. From these, we compute the scales of the model, shown in table 2, and hence we compute the typical values of the dimensionless parameters, as indicated in table 3.

The parameters *ε*, *ν* and *δ* are all small, and the Reynolds number *R**e*_{p} is moderately large. We therefore neglect the terms of order , *ν* and *δ*, but retain temporarily the singular term *ε**w*_{tt}, which is instrumental in describing the elastic waves which might propagate through the sample. The resultant model reduces to
3.31
where
3.32
and *β*=*O*(1). By eliminating *v*, we find that the pore gas pressure satisfies the nonlinear diffusion equation
3.33
This can also be written in terms of the scaled gas density *ρ*=*p*^{1/γ}, which provides the simpler equation for numerical purposes
3.34

The dimensionless pressure in the chamber above the sample is also denoted as *p*_{c}(*t*), and our choice of exponential in equation (3.24) yields the dimensionless form
3.35
where
3.36
Initially, the sample is at rest, the chamber pressure *p*_{c}=1 and this is also the initial pore pressure and solid stress through the sample. Up until the time when the sample lifts off the base, we take the solid displacement and gas flux to be zero at the base. At the surface of the sample, both solid stress and gas pressure equal the chamber pressure. Thus, suitable initial and boundary conditions for the sample displacement and gas pressure, up until fracture or lift off occurs, are
3.37
The fracture criterion (3.27) can be written in the dimensionless form
3.38
where
3.39

The quantity (1−*β*) represents the initial compressive displacement (at the surface) of the sample owing to the imposition of the initial chamber pressure, with *β*<1 representing a compression, and *β*>1 representing inflation. Mainly for simplicity, we assume that *β*=1, and in this case we have
3.40
The relaxation and fracture parameters *a* and *p** are given in table 3; typically, *a* is large and *p** is small. The consequences of these estimates are that diffusive relaxation of the pore gas pressure is relatively slow, and the fractures occur relatively close to each other in the sample. In the following sections, we describe the nature of the solutions to these equations.

### (c) Primary fracture

The basic picture we have of what happens following the release of the chamber pressure is as follows. The surface pressure at *z*=1 decreases (over a time scale of *O*(1/*a*)), and as it does so the pore pressure profile *p* decreases also, propagating a wave down into the sample. At the same time, the solid stress relaxes quasi-statically, since the acceleration term is small. To get an indication of the behaviour of the solid displacement *w*, we can reason as follows. Neglecting *ε* and taking *β*=1, equation (3.3)_{1} implies
3.41
representing force balance in the rock. At early times, the pore pressure has only changed near the surface, and its gradient with *z* is large; thus, *w*_{zz} is larger than *λ**w*, so that approximately , where
3.42

Consulting equation (3.42), we see that is a monotonically decreasing function of *z*, and thus the fracture criterion (3.38) is first satisfied at the base *z*=0 when *p*−*p*_{c}>*p**. In particular, without glue (when *λ*=0), the model predicts initial lift-off at the sample base.

An approximate correction to equation (3.42) at early times, or when *λ* is small, is given by
3.43
representing the solid strain in the rock. Since is a monotonically increasing function of *z*, we see that *w*_{z} increases with *z* near the base (where *p*_{z}≈0) and then decreases near the surface (where *p*_{z} is large). Thus, *w*_{z} reaches *p** first at an interior point below the surface, as observed.

### (d) Attenuation

In the above description, we have assumed a quasi-static elastic response of the sample based on the small value of *ε*. At closer inspection, this seems risky, because the general solution of equation (3.31), given *p*_{z}, will be the sum of the particular quasi-static solution, together with a pair of high-speed elastic waves which propagate up and down the sample. The consequent rapid fluctuations in *w* might then be expected to lead to the rapid disintegration of the sample without the ordered top-down sequence which is observed. In fact, so long as the surface pressure decreases over a time scale greater than , thus explicitly if
3.44
we can see that the quasi-static displacement profile satisfies not only the boundary conditions but also the initial conditions. A simple analogous problem is the forced oscillator equation
3.45
where it is known (Fowler & Kember 1996) that the slowly varying solution *w*∼*R*(*t*) is accompanied by the free oscillations of exponentially small amplitude. On this basis, we expect the accompanying elastic waves here to have negligible amplitude. This conclusion would not be true if the surface pressure were ramped down instantaneously (or very rapidly). Attenuation is present in the model through the terms of *O*(*ν*) and *O*(*δ*), but is not operative over the short time scales of the experiment.

### e) Explosion

Once an initial fracture has formed, either in the interior or at the base of the sample, the separated block can lift off. In this second phase of the motion, we suppose that the glue is also fractured, so that there is no longer any resistive traction at the walls. After separation, the gas flow within the separated blocks is still described by the previous model equations, but we must now consider the evolution of the crack widths. In addition, we address the issue of secondary fragmentation within the separated blocks. Suppose during the experiment we have a sequence of fractures formed at *z*_{1},*z*_{2},…, where the fracture at *z*_{i} opens into a crack whose upper and lower surfaces are denoted , . The (dimensional) gas flux into the crack is *ϕ*(*v*−*w*_{t}), and therefore the conservation of gas mass in the *i*-th crack yields the dimensional equation
3.46
where *ρ*_{i} is the gas density in the *i*-th crack. This equation determines the gas density, and thus the pressure, in the crack.

To render the equation dimensionless, we scale the variables as before, and in addition we scale the crack boundaries as 3.47 this leads to the dimensionless form of equation (3.46), 3.48 where .

In the *i*-th block , it is appropriate to put
3.49
(i.e. the displacement is relative to the base of the block), and then the crack width equation (3.48) becomes
3.50
bearing in mind that we may take . In addition, the momentum equation (3.31)_{1} can be written
3.51
Ignoring the displacement acceleration term (but not the block acceleration term), and bearing in mind that the total stress *σ*_{tot}=(1−*ϕ*)*σ*−*ϕ**p*=−*p*_{i} at the *i*-th crack, we see that the integration of equation (3.51) over the block leads to the approximate block acceleration equation
3.52
where is the (constant) length of the *i*-th block. If the *i*-th fracture forms at *z*=*z*_{i} at time *t*=*t*_{i}, then the initial conditions for equation (3.52) are
3.53

In summary, the gas pore pressure in the blocks continues to satisfy equation (3.33), with the boundary conditions that 3.54 The crack pressures are determined by solving equation (3.50), which takes the approximate form 3.55 and the block motion is determined by solving equation (3.52).

When the fracture criterion (3.27) is written dimensionlessly, it becomes 3.56 this is identical to equation (3.38), but the form above is more useful within the blocks. As in the blocks we have, from equation (3.51), 3.57 the criterion for secondary fracture within the blocks is simply 3.58 where 3.59 is the linear profile joining the pressures at the top and bottom of the block.

When an internal fracture is generated at *z*_{i}, then *p*_{i}>*p*_{i+1}, and, as *ε*/*δ*∼10^{−3}, the block accelerates rapidly upwards. As it does so, the crack pressure decreases according to equation (3.55), and we surmise that *p*_{i} rapidly approaches *p*_{i+1}, and the block speed reaches a constant, as is observed in the experiments.

The physical mechanism for this relaxation is akin to that of the removal of a piston from a tube. The necessity to fill the expanding void in the crack causes an effective suction which acts as a drag on the ascending block. Mathematically, we can provide some insight into this by the consideration of equation (3.55). Considering for simplicity the case of the first fracture, where *i*=1, is constant, and *p*_{i+1}=*p*_{c} is prescribed, we may write equation (3.55) in the form
3.60
where *f* is the porous gas flux from the adjacent blocks into the expanding crack. A simple idea of how *f* depends on the crack pressure *p*_{1} follows from elementary consideration of the way in which diffusion operates in the rock. As *p*_{1} reduces, the flux to the crack will increase, thus (∂*f*)/(∂*p*_{1})<0. (Of course, it is more complicated than this: *f* will also decrease with time, as the diffusive boundary layers grow into the blocks.) Assuming this, then equation (3.52) suggests . Blithely ignoring exponents and other detail, the structure of equation (3.60) is thus
3.61
whence, loosely,
3.62
and substituting back into equation (3.52) yields
3.63
which would indicate a damped, possibly oscillatory, approach to the equilibrium where *p*_{1}=*p*_{c}. This damping is found in our numerical results, as well as being implicit in the experimental results.

After (rapid) equilibration, *p*_{i}≈*p*_{i+1}, and the crack pressures will all decrease with the surface pressure *p*_{c}. Further internal fractures in the blocks may then occur, according to equations (3.58) and (3.59), if *p*−*p*_{c}>*p**, while the fracture in the lowest (attached) block mimics the initial primary fracture.

## 4. Numerical results

Equations (3.33) and (3.41) together with the associated initial and boundary conditions have been solved numerically using Matlab to determine the solid displacement *w* and gas pore pressure *p*. The gas diffusion equation (written as equation (3.34) for the density) was solved using small time steps and the pdepe solver for a glued-down sample of rock, and the boundary value problem (3.41) for *w* in the glued-in piece was solved for by using a Green function approach. Then, the fragmentation criterion was checked at each time step.

After fragmentation had occurred, equations (3.52) and (3.55) for the evolution of gap size and pressure were solved using the Runge–Kutta ode45 solver, again with small time steps, and the gas fluxes from adjacent rock segments were calculated as required by ode45. The process was continued, checking for further fracture in every piece of rock at every time step, and increasing the number of coupled differential equations to be solved each time a new break was detected.

Parameter values are essentially as given in table 1, but with slightly different values of *ϕ* and *σ*_{Y}, as noted in the caption of figure 4. The value of *ϕ* is that appropriate for the experimental data in figure 3.

Figure 4 shows the evolution with time using snapshots of the gas pore pressure *p* and the solid compressive stress −*σ*. As the surface pressure at 60 mm lowers, a diffusive boundary layer profile (solid line) of the gas pressure migrates into the column. The compressive stress in the solid responds quasi-statically, and according to equations (3.29) and (3.43) is given (approximately if *λ* is small) by
4.1
as *ϕ*/(1−*ϕ*)≈1, this can be seen in figure 4*a*. The slight minimum in the solid stress is due to the non-zero value of *λ*. Consequently, fracture first occurs near this minimum, as shown in figure 4*b*, where, immediately following fracture, the solid stress in the basal part of the column relaxes rapidly towards the gas pore pressure, which is still approximately *p*_{c}. Continued diffusion of the gas pressure draws down the pressure in the ejected block, while the pressure and stress profiles in the basal column repeat their evolution as in figure 4*a*. This process repeats itself, as indicated in figure 4*d*,*e*, causing a sequence of fractures and consequent ejection of the blocks. Figure 5 shows the consequent pattern of fractures, plotted similarly to figure 3, with which it can be (favourably) compared.

## 5. Discussion

In our model, fragmentation is largely controlled by two quantities. The first of these is the yield stress *σ*_{Y}. As indicated by equation (3.39), if this is less than the initial confining pressure, then fragmentation is probable, providing the contained gas cannot escape sufficiently rapidly. In the experiment, the time scale for gas escape is *t*_{0}, and this is determined through the gas pore velocity scale *v*_{0}, defined in equation (3.28). Thus, for low permeability, *v*_{0} is low, *t*_{0} is high and in this case fragmentation is likely to occur.

While our model is designed to simulate the fragmentation experiments, which are themselves designed to provide insight into explosive volcanic eruptions, it is less easy to draw theoretical conclusions which can be applied to volcanic processes. In the experiment, the rock is already vesicular and porous when depressurization commences.

However, this might not always be the case in volcanic eruptions. At an eruption, the sudden depressurization by the removal of an overburden or volcanic plug can cause volatile exsolution to a very different extent, depending on the amount of depressurization and on magma properties such as viscosity, diffusivity and volatile disequilibrium. For instance, at sustained Plinian-type eruptions, an equilibrium can be assumed between magma ascent and depressurization-driven volatile exsolution. This results in a highly vesicular magma; consequently, magma disintegration (fragmentation) can occur, if the stresses generated in the melt phase by the vesiculation-driven expansion exceed the yield strength of the magma.

A completely different picture can be drawn for the typically short-lived Vulcanian-style eruptions. Here, the erupting magma is commonly highly viscous and, given the usually short time scale between plug removal and the consequent initiation of sudden depressurization and explosive ejection of the magma, minor additional volatile exsolution driven by the rapid depressurization will occur. The laboratory experiments described and modelled here bear most resemblance to the latter eruption style. Still, the details of the fracturing process in the volcanic situation are assumed to be somewhat different from that in the experiment, and the details of the model will thus also be necessarily different.

Even if we restrict our attention to the present experiment, there are issues in our model which are at least questionable. The precise form of the interphasic drag terms used in equation (3.16) is based on experimental results obtained under steady equations. In the present case, it is probable that inertial effects may be important (Chojnicki *et al.* 2006), and if, for example, we introduced an inertial term proportional to ∂**V**/∂*t* into equation (3.16), the resulting model (3.33) for *p* would become hyperbolic rather than parabolic. It seems unlikely that this would substantially alter the results, however.

The wall friction term (3.21) was introduced into the model because, without it, the model always predicted the initial lift-off at the base of the column. And, indeed, this is what is seen, because if the rock is not glued in place, then basal lift-off occurs first. Initial efforts to model the effect of the wall glue used a viscous formulation, but it is thought that the elastic form better suits the resistance prior to rupture. That is to say, once the rock ruptures, the glue, having less strength, ruptures also, and offers no further resistance. We have chosen a value of the glue shear modulus indicating its weaker strength, and comparable to the typical values for polymer glues.

An issue of interest concerns the robustness of the results we have obtained. Although more detailed explorations will be conducted in future work, we can state that the model fits the experimental results both qualitatively and quantitatively without the need for special choices of parameters. The results obtained depend in their quantitative predictions on the choice of the parameters *β*, *λ*, *p** and *a*. While the values of *β* and *λ* are open to doubt, the choices presented here were in fact our original choices. The same applies to the choice of yield stress *σ*_{Y}. Variation of these parameters does vary the results, but not to any marked degree.

## 6. Conclusions

We have constructed a mathematical model of the fragmentation process in the experiments described by Spieler *et al.* (2004). In our model, the rock is elastic and subject to a fracture criterion which depends on the local effective tensile stress in the rock. Following a fracture, the ejected fragment is expelled by the acceleration applied by the higher gas pressure below the fragment. However, the expulsion of the fragment and consequent widening of the fracture beneath requires that the pore gas be sucked from the fragment to fill the expanding fracture, and this suction causes a deceleration of the ejected fragment. Because the gas depressurization wave propagates relatively slowly, the effective tensile stress reaches a maximum at a finite depth in the sample, and thus fracturing occurs first at finite depth. As each fragment is ejected, the crack pressure rapidly decreases towards the surface pressure, and the process is repeated. In addition, secondary fractures can occur within the ejected samples.

All of this behaviour is seen in the laboratory experiment, and our numerical solutions of the model indicate good qualitative and quantitative agreement with the experiments. In future work, we will explore in more detail how the numerical results compare with the range of experimental results, and in particular whether some simple rules of thumb, for example concerning fracture threshold porosity, can be deduced from the model, and we will also explore whether yet simpler representations of the model solutions can be found in certain parameter regimes.

## Acknowledgements

We acknowledge the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland mathematics initiative grant 06/MI/005. The experimental work was funded by the Japan Society for the Promotion of Science (no. P05–318). We are thankful to Mie Ichihara for her help in building the laboratory as well as for useful discussions. Further, we thank Kei Kurita for sharing the high-speed camera with us. We also thank Peter Howell for his customary illuminating insights.

## Footnotes

↵1 The notation

*z*∼*l*, for example, is an abbreviation for the process of writing*z*=*l**z**, substituting into the model, and then in the resulting dimensionless equations, omitting the asterisks on the starred variables.- Received July 21, 2009.
- Accepted October 9, 2009.

- © 2009 The Royal Society