Effective models for the evolution of magma permeability are key to understanding shallow magma ascent and eruption dynamics. Models are generally empirical constructs, commonly focused on monodisperse systems, and unable to cope with the foam limit at high porosity. Here, we confirm that bubble size distributions in high-porosity pyroclasts are highly polydisperse. We combine collated experimental data and numerical simulations to test and validate a theoretically grounded percolation model for isotropic magma permeability, which accounts for the effect of polydispersivity of bubble sizes. We find that the polydispersivity controls the percolation threshold. It also serves as essential input into the scaling of permeability that is required to achieve universality in the description of permeability. Our model performs well against collated published data for the permeability of high-porosity volcanic rocks. We then extend this model to predict the viscous and inertial contributions to fluid flow that are required to model magma outgassing in all regimes. Our scaling relationship holds across the full range of porosity, from the percolation threshold to the open-foam limit.


Permeability can exert a first-order control on the explosive potential of magmas rising through the crust (e.g., Mueller et al., 2008; Degruyter et al., 2012; Cassidy et al., 2018). Yet models for magma permeability remain empirical, semi-empirical, or limited to specific systems (Klug and Cashman, 1996; Saar and Manga, 1999; Mueller et al., 2008; Farquharson et al., 2015; Gonnermann et al., 2017). These models are typically calibrated against data sets collected for natural or experimental volcanic materials quenched from magmas (Klug and Cashman, 1996; Mueller et al., 2008; Colombier et al., 2017), which, taken together, provide a rich data resource for testing new, theoretical models. The scatter and range that can exist in permeability data for materials produced in nature and in the laboratory are prodigious due to the wide variety of microstructural geometries possible in magmas (e.g., vesicular, fractured, granular). Therefore, the most powerful and universal laws for how permeability varies with porosity have been built with a direct accounting for the microstructural origins (Martys et al., 1994; Wadsworth et al., 2016; Vasseur and Wadsworth, 2017; Giachetti et al., 2019).

Multiple discrete bubble nucleation events, continuous nucleation of bubbles, bubble coalescence, and differential bubble growth rates can all contribute to polydisperse bubble size distributions in magmas and pyroclasts (e.g., Blower et al., 2001). If the full size distribution of vesicles is constrained, then polydispersivity can be parameterized by a single metric (Torquato, 2013):
where graphic is the nth moment of the radius distribution with graphicgraphic and graphic corresponding to the mean, variance, and skewness respectively. Collations of measured vesicle size distributions in erupted pyroclasts of pumice or scoria, and conversion of them to S via Equation 1, demonstrate that for all natural volcanic products, the vesicle sizes are highly polydisperse (Fig. 1). Despite this fact, no theoretical treatment of permeability exists for polydisperse bubble sizes that is valid from the lower limit of the percolation threshold (i.e., the point at which a cluster of bubbles first spans the system edge to edge) to the upper limit of foam at high porosity. In this contribution, we use a combination of collated published experimental data and novel numerical simulation results to test and validate a scaling for permeability that fully accounts for polydisperse bubble sizes. We show, in particular, that for the low values of S typical of natural systems, our new model predicts values of permeability that are many orders of magnitude different from those yielded by existing modeling approaches. We focus on the prograde path of bubble growth, in which the porosity is an increasing function of time.


One highly effective approach to validating a model for magma permeability is to approximate a high-porosity bubbly magma as a system of overlapping spheres of the same porosity (Blower, 2001; Vasseur and Wadsworth, 2017). We adopt this geometric approach in the design of our numerical simulations. We place spheres of radius distribution p(R) randomly in a periodic cube domain of edge length L (cf. Rintoul and Torquato, 1997; Vasseur and Wadsworth, 2017). We allow the spheres to overlap freely as we add them one by one. We define the spheres as the “bubble” phase, with porosity ϕ, and for this two-phase system, the inter-sphere phase represents the “groundmass” phase, with volume fraction 1 – ϕ. In contrast to previous models that have used monodisperse spheres, we here use a power-law size distribution such that the probability that a sphere has a radius between R and R + dR is p(R) = αR−(α+1), where α > 3 for three-dimensional (3-D) domains (α > 2 for the two-dimensional [2-D] equivalent). We vary graphic across a wide range and then convert the resultant distributions of sizes into the polydispersivity S. For the probability density function p(R) defined here, graphic = α/(α – n) for n = 1, 2, and 3. To illustrate our method, in Figure 2 we present rendered 3-D volumes in which we show only the fluid (non-solid) phase in a sub-volume.

For each new sphere we add to the domain, we use a Monte Carlo union-find algorithm to check if there is a cluster of connected spheres that span the domain from face to face (Newman and Ziff, 2001). Each new sphere increases the total number density by 1/L3 to N/L3, where N is the total number of spheres up to that point. As we increase N, there is a point beyond which the system has a connected phase across its length, which equates to a percolation porosity threshold ϕc. By repeating this process 10,000 times for each S and N, we can find the probability that ϕ ≥ ϕc, which we term Π. The variation of Π with ϕ is sigmoidal, to which we fit the function: Π = {1 + tanh[[ϕ – ϕ′c(L)]/Δ(L)]}/2. The fit parameters Δ(L) and ϕ′c(L) are the width of the transition and the effective percolation threshold, respectively. In the GSA Data Repository1, we plot these two parameters against one another, and following Sasidevan (2013), we can find the limiting ϕc as Δ → 0 (equivalent to an infinite domain size L3 → ∞). This represents constraint of the 3-D percolation threshold for overlapping spheres for any S.

We calibrate and check our method for constraining ϕc by performing the same analysis but in 2-D and comparing directly with the results of Sasidevan (2013), as well as by comparing our 3-D results in the monodisperse limit (for S = 1) with values reported previously (Rintoul and Torquato, 1997; Lorenz and Ziff, 2001; Vasseur and Wadsworth, 2017). In 2-D, S = graphic

For each domain we generate for ϕ > ϕc, we use the LBflow code (Llewellin, 2010) to simulate fluid flow across the system through the sphere phase. We use input conditions for which the Reynolds (Re) and Mach (Ma) numbers are ≪1, and for which Darcy’s law applies at steady state. We impose a driving pressure gradient of ∇P = 0.01 Pa·m−1, and find 10−10 < Re < 10−6 and 10−14 < Ma < 10−9. We use a steady-state criterion such that the average speed across the entire lattice must not vary by more than a relative factor of 10−5 over 50 time steps, two times consecutively (Llewellin, 2010). The output is a simulation value of permeability k for each domain of a given ϕ and S. Example steady-state snapshots of the fluid-flow vector field are shown in the insets to Figure 2. Taken together, these methods provide us with constraint of permeability and the percolation threshold as a function of bubble polydispersivity. Finally, we measure the specific surface area for each domain using a Lewiner marching cubes algorithm (Lewiner et al., 2003; Vasseur and Wadsworth, 2017).


We use our method to predict the percolation threshold, and we find that it is dependent on the polydispersivity of bubble sizes. Our results show that in the monodisperse limit (S = 1), ϕc ≈ 0.28934 ± 0.00033 (or equivalently ϕc ≈ 0.67630 ± 0.00001 in the 2-D case). These values are within 0.08% of predictions given in previous work (Rintoul and Torquato, 1997; Lorenz and Ziff, 2001; Vasseur and Wadsworth, 2017). This result allows us to confidently explore polydisperse systems (S < 1). For S < 1, the percolation threshold diverges nonlinearly from the monodisperse limit toward very high values (ϕc ≈ 1) as the polydispersivity rises (S → 0), which is shown in Figure 3A. This implies that polydisperse systems such as those found in natural pumice and scoria would be characterized by a higher percolation threshold than equivalent monodisperse systems, which verifies hypotheses made previously (Colombier et al., 2017). In the Data Repository, we collate these numerical results.

Our results also reveal that the permeability-porosity trend k(ϕ) is highly dependent on how polydisperse the growing bubble population is. To show this, we look for a simple, universal scaling law that can capture all of our results. The most successful universal scaling relationships for k(ϕ) for random heterogeneous media rely on the specific surface area of the pore network, s, and have the form k ∝ 1/s2 (Martys et al., 1994; Carman, 1997; Wadsworth et al., 2016). Torquato (2013) showed that for overlapping polydisperse sphere systems:
where graphic in the monodisperse limit, and graphic is the normalized specific surface area. This relationship agrees well with our simulation result across all S (Fig. 3B). A scaling for k(ϕ) that relies on s has the advantage that the bubble size is not a direct input, which for polydisperse systems can be challenging to measure or constrain (Blower et al., 2001), and relies instead on bulk properties only. We use a simple form
where b is the percolation exponent. For bubbly geometry relevant here, b = 2.4 has been predicted theoretically (Feng et al., 1987; where those authors use the term “inverted swiss cheese” to describe this geometry). We find excellent agreement between our measured k and the prediction of Equation 3 (Fig. 4). In Figure 4, we show the dimensionless result for a normalized permeability graphic as well as the dimensional results from our simulations.

Previous models for the permeability of systems with polydisperse bubble sizes do not include an account of ϕc (Costa, 2006). Here we demonstrate that ϕc is critical at high polydispersivity (Fig. 4A) and is clearly important in experimental data (Colombier et al., 2017). Much of the experimental data for magma or volcanic materials from which k and ϕ have been constrained are neither associated with a measured s, nor with the information required to predict s from the size distribution of pores, vesicles, or bubbles (exceptions include Farquharson et al. [2015] and Baker et al. [2019]). Therefore, in order to further test our model against experimental data, we have chosen to scrutinize the family of open-celled foams in ceramic materials. Collating data for the permeability of ceramic solid open foams for which either s or p(R) were measured (Philipse and Schram, 1991; Acosta et al., 1995; Innocentini et al., 1998, 1999; Moreira et al., 2004), we can apply our model directly. We find that our model accurately predicts experimental results in the foam limit without requiring empirical adjustment (Fig. 4A) to within approximately an order of magnitude at the highest porosities (with increasing accuracy for lower porosity). The ceramic foams used here have textures that are similar to those of, and are relevant for, high-porosity pyroclasts up to and including the reticulite open-foam limit (see Innocentini et al., 1998, their figure 2).

We also compare our model with a large compiled data set for the permeability of volcanic rocks (Colombier et al., 2017, and references therein). In Figure 4B, we focus on the permeability of rock products of explosive eruptions, and show that when our model is solved for the global mean S = 0.28 (Fig. 1), we predict that ϕc = 0.44 (Fig. 3) and that the trends of k(ϕ) for reasonable mean vesicle sizes are consistent with the data. Therefore, we propose that our model captures the principal features of prograde degassing and the onset and evolution of the propensity for outgassing, as it is recorded in eruptive products here.


The parameterization of permeability k given here is useful only for predicting outgassing rates in the low-Reynolds-number range where Darcy’s law applies. However, outgassing through magma can have high velocity at shallow levels, such that the Reynolds number is unlikely to be low. Therefore, an additional constraint of the inertial contribution to flow through porous magma is required for our model to be of widest utility.

The average steady-state fluid velocity graphic that arises from a given pressure gradient ∇P can be predicted for any Reynolds number by the Forchheimer equation (Whitaker, 1996):
where μ and ρ are the fluid properties viscosity and density, respectively, and kI is the inertial permeability. kI can be predicted from k using a simple scaling kI = Akc where A ≈ 1010 and c = 3/2 (Zhou et al., 2019). This simple scaling coupled with Equations 14 and constitutive laws for the fluid properties provides a full quantitative description of fluid flow rates through porous magmas.

In our approach, we have made the simplifying assumption that a polydisperse bubbly magma can be approximated as a system of overlapping spheres. While this model approach is used widely (e.g., Blower, 2001; Vasseur and Wadsworth, 2017), the geometry imposed neglects the effects of bubble-bubble flattening and deformation prior to coalescence (Gonnermann et al., 2017) and the development of tube pumice (Dingwell et al., 2015), which have been observed and modeled. We note that textures that seem to record the capillary resistance of magma bubbles to coalescence may be one reason for the high percolation threshold invoked to explain magma bubbly permeability (Colombier et al., 2017; Giachetti et al., 2019). Similarly, it should be noted that the model neglects the effect of high crystal volume fractions, which can influence the permeability in a complex way (Lindoo et al., 2017). Nevertheless, we show quantitatively here that bubble polydispersivity also has a first-order effect on the percolation threshold (Fig. 3A) that has not been considered previously.

Coupling our model with Equation 4 provides a framework in which magma outgassing can be predicted. If this framework were embedded in a numerical model for magma expansion and compaction (e.g., Gonnermann et al., 2017), then the nonlinear feedback between gas pressure–driven magma porosity change and the local gas flow rates would be within our ability to predict.


Polydispersivity of bubble sizes in magma foams is crucial to predicting percolation threshold, pore surface areas, and permeability. Using numerical simulations, measurements on quenched volcanic materials, and constraints from analogous ceramic foams as a validation step, we provide a simple-to-use model for magma permeability. Using the Forchheimer equation, our model can be used to predict the local outgassing rates from open-celled magma foams in the shallow subsurface, of wide utility to gas monitoring efforts and to models that seek to predict the explosive potential of magma in various scenarios worldwide.


We acknowledge support of the European Research Council (ERC) Advanced Grant on Experimental Access to Volcanic Eruptions: Driving Observational Potential (EAVESDROP, grant 834225). We thank the Leibniz Supercomputing Centre (LRZ) of the Bavarian Academy of Sciences and Humanities (Garching bei München, Germany) for computational support, as well as Jamie Farquharson, Thomas Shea, and one anonymous reviewer for constructive comments.

1GSA Data Repository item 2020152, rigorous constraint of the percolation transition in 2-D and 3-D and the full reference list pertaining to Figure 1, is available online at http://www.geosociety.org/datarepository/2020/, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY license.