To illustrate the advances made in permeability calculations combining X-ray microtomography and lattice Boltzmann method simulations, a sample suite of different types of pumices was investigated. Large three-dimensional images at high spatial resolution were collected at three different synchrotron facilities (Elettra, SLS, and ESRF). Single phase gas flow simulations were done on computer clusters with a highly parallelized lattice Boltzmann code, named Palabos. Permeability measurements obtained by gas flow simulation were compared to lab measurements of pumices produced by the Kos Plateau Tuff eruption to validate the method. New permeability data for pumices from other silicic volcanic deposits is presented, and an empirical model for permeability is tested using geometrical and topological data, i.e., tortuosity, specific surface area, and total and connected porosity.


Macroscopic transport of fluids in porous media has numerous applications in various fields of geosciences (e.g., Turcotte and Schubert, 2002). In the case of laminar flow of a single phase fluid, this transport is described by Darcy's law, 
with x the direction of flow, −dP/dx the pressure gradient, μ the fluid viscosity, and U the average velocity per unit area. The permeability k is determined by the characteristics of the porous medium (geometry and topology) at the microscopic scale. Empirical models (e.g., Kozeny-Carman equation) relate geometrical properties such as specific surface area, and channel shape together with topological quantities such as tortuosity and connected porosity quantitatively to permeability (Bear, 1972). However, different types of rock require different types of empirical equations, and therefore a detailed examination of rocks combined with model testing is needed.

Both X-ray tomographic imaging (Flannery et al., 1987) and lattice Boltzmann numerical techniques (Wolf-Gladrow, 2000; Succi, 2001; Chopard et al., 2002; Sukop and Thorne, 2006) emerged at the end of the 1980s, and it was soon realized that these techniques were very compatible for permeability studies (Ferreol and Rothman, 1995; Heijs and Lowe, 1995; Auzerais et al., 1996). The lattice Boltzmann method (LBM) is excellent for modeling flow in complex geometries, and X-ray tomography provides the three-dimensional (3D) geometries needed as boundary conditions for the flow. Compared to laboratory measurements, this numerical technique has the great advantage that the velocity distribution can immediately be linked to the microscale geometry of the porous medium. The main drawback of the method, simulating a large enough volume (representative elementary volume [REV]; Bear, 1972) at a high enough spatial resolution, becomes less prominent as X-ray tomography and computer power (clusters) advance and become more available. Today, X-ray tomography facilities at synchrotron beamlines easily produce >10243 voxels of data, in combination with submicron spatial resolutions (Stampanoni et al., 2006; Weitkamp et al., 2010), while lattice Boltzmann codes are highly parallelized to be able to handle such large volumes (Harting et al., 2005; Fredrich et al., 2006).

The development of permeability in the volcanic conduit during an eruption is believed to play an important role in the shift between eruptive styles (Eichelberger et al., 1986; Woods and Koyaguchi, 1994; Gonnermann and Manga, 2007; Mueller et al., 2008). Pumices (highly vesicular volcanic rock) record the state of permeability in the volcanic conduit just prior to fragmentation, provided they are quenched fast enough and are highly viscous (Thomas et al., 1994). Hence, they give us insight into the degassing mechanisms during a volcanic eruption.

Here, we combine X-ray tomography with lattice Boltzmann numerical simulations of fluid transport to obtain geometrical properties and permeability estimates of pumices. In the line of previous studies (Wright et al., 2006; Polacci et al., 2009) images were obtained at very high spatial resolutions and flows were simulated on large enough grids to capture the REV. We focus on pumices produced by the Kos Plateau Tuff eruption and compare the results to published laboratory results (Bouvet de Maisonneuve et al., 2009). We also present the first permeability results for silicic pumices of five other explosive eruptions (Okmok Volcano, USA; Nevado de Toluca, Mexico; Mont Dore, France; Cappadocia, Turkey; Mont Pelée, Martinique, Lesser Antilles). These analyses are accompanied by the quantification of geometrical and topological characteristics of the pore space, i.e., total and connected porosity, tortuosity, and specific surface area to test an empirical model suggested by Degruyter et al. (2010).


Samples from six different silicic volcanic deposits were selected for analysis. The general characteristics of each of these deposits are summarized in Table 1. The main set of samples is composed of pumices from the 160 k.y. B.P. Kos Plateau Tuff (KPT) eruption. This eruption produced various types of pumices that were well characterized in previous studies (Allen, 2001; Bouvet de Maisonneuve et al., 2009; Degruyter et al., 2010). Tube, frothy, and microvesicular (μVes) pumices were used to validate the method (Fig. 1). The samples were taken from the same clasts as the study of Bouvet de Maisonneuve et al. (2009). The rest of the samples were taken from Okmok volcano, Alaska (AOK; Burgisser, 2005; Larsen et al., 2007); Nevado de Toluca volcano, Mexico (ATO; Burgisser and Gardner, 2006); Cappadocia, Turkey (CAP; Druitt et al., 1995); Mont Dore, Massif Central, France (GNP; Mossand, 1983); and Mont Pelée, Martinique, Lesser Antilles (MT; Martel and Poussineau, 2007). From each clast, one or several 5 mm cylindrical cores were drilled as this is the most optimal shape to scan a maximum volume of the sample (Ketcham and Carlson, 2001). In total 12 samples were prepared from 9 clasts (Table 2).


Modern synchrotron radiation facilities produce a high flux source of monochromatic X-rays, ideally suited for micrometer and submicrometer 3D imaging. Together with large cameras (≥ 1024×1024 pixels), this allows the quantification of large enough volumes of material. For this study, X-ray microtomographic imaging was performed at three different synchrotron facilities: the SYRMEP beamline at Elettra (Trieste, Italy; Polacci et al., 2009), the TOMCAT beamline at SLS (Villigen, Switzerland; Stampanoni et al., 2006), and the ID-19 beamline at ESRF (Grenoble, France; Weitkamp et al., 2010). At all these facilities the beam energy was varied between 20 and 25 keV. The samples were rotated 180° and between 1000 and 3000 projections were taken per scan. At the SYRMEP a field of view (FOV) of 7.9 × 3.7 mm2 and an effective pixel size of 7.74 × 7.74 μm2 were used. At both the TOMCAT and the ID-19 a FOV between 0.75 × 0.75 mm2 and 1.5 × 1.5 mm2 was used with submicron effective pixel sizes between 0.37 × 0.37 μm2 and 0.74 × 0.74 μm2 (see Table 2). For the submicron scans we used the “local tomography” technique to avoid reducing the size of the fragile samples (Faridani et al., 1992; Faridani et al., 1997). With this technique the camera zooms into the sample, leaving it partly outside the field of view. Noise is suppressed by rotating the sample continuously. Some of the samples were long enough in the direction perpendicular to the scan plane to capture several subvolumes. In total 27 subvolumes were scanned and analyzed (Table 2). The reconstruction of the images was done using the in-house software developed at the respective synchrotron beamlines. All images were converted to a stack of 8-bit grayscale images, and were then cropped into cubes for further analysis (Fig. 2A). Tube pumice images, for which the deformation was slightly tilted with respect to the z direction, were realigned and recropped with ImageJ.

Image segmentation, the process by which a grayscale image is converted into black and white voxels, was done using the Indicator Kriging method of 3dma-rock (Oh and Lindquist, 1999) and/or the Seed threshold method of Blob3D (Ketcham, 2005) giving a stack of binary images (Fig. 2B). Segmentation can introduce some errors. Hence, a postsegmentation clean up is performed by removing all objects with an equivalent sphere diameter of < 3 voxels (total of 27 voxels). Therefore the accuracy of the measurements will be three times the effective voxel size. Volume renderings of the segmented images were done with Paraview (Fig. 1). Total porosity was measured by counting the number of white voxels and dividing by the total volume. The connected porosity was calculated in each perpendicular direction by labeling the separated connected clusters of white voxels with ImageJ. The medial axis calculated in 3dma-rock can be used to search for the shortest paths between two end faces (Lindquist and Venkatarangan, 1999). The average of this value divided by the length between the two faces gives us the tortuosity. The surface area of the porous medium was calculated on the segmented volume with a Marching Cubes algorithm (Lorensen and Cline, 1987), which is incorporated in the 3dma-rock code. This algorithm creates triangulated surfaces from the original voxel data and therefore avoids overestimating the surface area if calculated directly from the voxel surfaces. Dividing this by the total volume gives the specific surface area. Total porosity, tortuosity, and specific surface area are listed in Table 3. All connected porosities were found to be within 3% of the total porosity.


Lattice Boltzmann simulations were carried out using Palabos (Latt, 2009). Palabos is an open-source software for computational fluid dynamics using the LBM. As it is highly flexible and parallel, it can handle a wide variety of fluid dynamics problems. The LBM itself has been thoroughly described and validated in the literature (Wolf-Gladrow, 2000; Succi, 2001; Chopard et al., 2002; Sukop and Thorne, 2006). In this study, a single phase gas flow is simulated through the obtained 3D images. Using the segmented images as input data, the walls of the porous medium are converted to bounce-back boundary conditions. The geometries were padded with walls in the direction perpendicular to flow in order to eliminate the need for periodic boundaries. The standard Bhatnagar-Gross-Krook (BGK) collision operator is applied with a D3Q19 lattice. Initially the fluid velocity is set to zero and fluid movement is induced by holding a fixed pressure gradient between the inlet and outlet. Steady state is reached when the standard deviation of the average energy, as measured over a fixed number of time steps, falls below a given threshold value (we chose 1000 time steps, and a threshold value of 10−4). To ensure the flow was in the laminar regime (Re << 1), it was checked if the permeability stayed constant when the applied pressure gradient was varied over several orders of magnitude. Furthermore, as discussed by Ferreol and Rothman (1995), the channels through which the fluid flows must be sufficiently resolved to avoid nonhydrodynamic effects. This was verified by refining the grid of parts of the analyzed volumes and checking that the permeability stayed constant for these refinements. The code is available at the Palabos website with a tutorial (Latt, 2009).

Steady-state results of the model are analyzed to obtain permeability. From the velocity distribution throughout the 3D volume Equation (1) is used together with the given pressure gradient and fluid viscosity to solve for permeability k in nondimensional lattice units. This value is easily converted to real world units by multiplying with the square of the effective length of a voxel side. For each volume, three permeability values are obtained, one in each orthogonal direction (in total 79 points; Table 3). Figure 2 shows the steps from a cropped grayscale image (A) to binary image (B) and finally to a velocity distribution simulated with the LBM (C) for a slice of sample AOK1A. Figure 3 shows the complex streamlines of the flow through frothy pumice (KPT32c). Alternatively, these streamlines can be used for tortuosity calculations (e.g., Matyka et al., 2008).

A large amount of memory needs to be allocated for the volumes that are analyzed. Therefore simulations were run on several computer clusters: Myrinet (University of Geneva, Switzerland, 32 nodes, two processors, and 2 Gb RAM per node), Atlas (Georgia Tech, Atlanta, USA, 128 nodes, eight processors, and 8 Gb RAM per node), Phoebus (CNRS, Orléans, France, 42 nodes, eight processors, and 8 Gb RAM per node), and Cadmos (EPFL, Lausanne, Switzerland, 4000 nodes, four processors, and 4 Gb RAM per node). Between 8 and 128 processors were used for each simulation, except on Cadmos where 2048 processors were used, and a run took between 10 and 30 h.


Effect of Resolution and Analyzed Volume on Measurements

The need for having an REV at high enough spatial resolution for adequate characterization of a porous medium is demonstrated on sample AOK1A following Zhang et al. (2000). To test spatial resolution effects, the segmented volume is rescaled by transforming each 2 × 2 × 2 voxels to 1 voxel with a certain rule. Two different rules are applied: if the 2 × 2 × 2 voxels have more than four white voxels, it is converted to one white voxel, otherwise it becomes black (rule 1); if the 2 × 2 × 2 voxels have more than seven white voxels, it is converted to one white voxel, otherwise it becomes black (rule 2). For rule 1 the wall becomes less resolved with increasing the effective voxel size, while for rule 2 it is the bubble network that loses resolution. Hereby we artificially increase the effective voxel size from 0.553 μm3, over 1.13 μm3, 2.23 μm3, and 4.43 μm3 to 8.83 μm3, while maintaining the physical size of the analyzed volume. The effect of losing resolution is demonstrated in Figure 4: rule 1 results in increasing permeability and connected porosity, while rule 2 has the opposite effect. For this example the decrease in spatial resolution becomes significant at an effective voxel size of 2.23 μm3 as the measurements of connected porosity and permeability deviate drastically from the ones done at an effective voxel size of 0.553 μm3. Note that the deviation of the connected porosity is solely due to misrepresentation of the solid-void space, while that of the permeability can arise from both the misrepresentation of the solid-void space and numerical errors induced by not resolving the flow in the channels properly (Ferreol and Rothman, 1995).

To find the REV the original segmented volume of 10243 voxels is divided into eight subvolumes of 5123 voxels, while maintaining the effective pixel size of 0.553 μm3. Then a subvolume of 5123 voxels is selected around the center point of the original volume, and is again divided into eight subvolumes of 2563 voxels. This procedure is repeated; eight volumes of 643 voxels are reached. The effect of losing volume on permeability and connected porosity is demonstrated in Figure 5. For this example the loss of volume becomes significant at a volume of 2563 voxels as the spread in measured connected porosity becomes larger than 0.1 and the spread in permeability becomes larger than one order of magnitude.

Validation of the Method on KPT Samples

The permeability measurements done with LBM simulations compare well with laboratory measurements done on KPT samples (Bouvet de Maisonneuve et al., 2009; Figs. 6A and 6B). Laboratory measurements were done on ~ (25 mm)3 volumes while volumes analyzed here varied between (0.37 mm)3 and (2.5 mm)3. Despite this vast volume difference of 4 to 6 orders of magnitude, similar spreads within the same pumice type and in between pumice types are found. This remarkable agreement between simulated and measured permeability shows that (1) the tomographic volumes qualify as REVs and (2) the pumices were imaged at a high enough spatial resolution. In fact, the points with the greatest mismatch (i.e., frothy type) are also the smallest grids analyzed (Table 2). This shows the necessity of doing large-scale simulations (> 5123 voxels) to be within the REV (Zhang et al., 2000).

Tortuosity has a clear influence on permeability of volcanic rocks (Le Pennec et al., 2001; Bernard et al., 2007; Wright et al., 2009). The trends in Figures 7A and 7C confirm qualitative conclusions of previous studies (Bouvet de Maisonneuve et al., 2009; Degruyter et al., 2010). The frothy pumice has a higher tortuosity leading to a lower permeability than the tube measured parallel to vesicle elongation despite its higher porosity. The anisotropic permeability within tube pumice is also clearly influenced by tortuosity differences caused by bubble elongation. The μVes pumice is characterized by having a relatively low tortuosity and low permeability. This is explained by the specific surface area measurements, which scale inversely with permeability (Saar and Manga, 1999), showing one order of magnitude difference between μVes and the other pumice types (Table 3; Fig. 8).

Other Deposits

Large volumes of 10243 voxels (greater than one billion grid points) were analyzed for samples from the other deposits (Table 2). Qualitative inspection of the scans ensured there were enough connected bubbles per analyzed volume except for the GNP sample. Although more data is needed for each deposit, these results shed a first light on the permeability development during these eruptions (Fig. 6C). A wide scatter can be observed in between the deposits as well as within the same deposit, similar to other studies (Mueller et al., 2005; Gonnermann and Manga, 2007). AOK and CAP show the influence of the bubble deformation as the permeability measured parallel to vesicle elongation is higher than the one measured perpendicular to vesicle elongation. Also GNP samples showed high deformation, but no connections perpendicular to vesicle elongation were found for these samples, confirming we did not scan a REV for this sample (Table 3). Lastly the ATO and MT sample show isotropic permeability similar to KPT frothy and μVes pumices. Similar tortuosity-permeability trends can be observed for samples from the other deposits (Figs. 7B and 7D), as the pumice with deformed bubble networks (AOK, CAP) has a larger spread than pumices with no distinct bubble deformation (ATO, MT). The GNP data point plots in the range of the KPT tube pumice measured parallel to vesicle elongation.


This study was designed to improve existing permeability models, which use a power-law relationship between porosity φ and permeability, 
with C and m fitting constants (e.g., Klug and Cashman, 1996; Blower, 2001; Mueller et al., 2005). More sophisticated models can be presented if the fitting constants can be understood better. Previous studies have shown that the cementation factor m can be correlated with tortuosity τ and connected porosity φc using Archie's law (τ2= φc−m), and C is related to the square of a characteristic length scale d of the porous medium (Costa, 2006; Bernard et al., 2007; Degruyter et al., 2010). A model used by Degruyter et al. (2010) replaces all fitting constants with measurable quantities: 
As this is an equivalent channel model, only the connected porosity φc is considered, and d is the equivalent channel diameter. Length scales such as the average diameter of bubble apertures (e.g., Blower, 2001) can be used to estimate the equivalent channel diameter. Here we use the specific surface area as a length scale, which is related to the channel diameter as d =c/s. Applied to Equation (3) we get 
A spread for m was obtained by combining the tortuosity and the connected porosity data (Archie's law; Fig. 5). Permeability predicted by equation 4 was compared with permeability obtained from the simulations. Ranges predicted by the empirical model compare relatively well with observed values for frothy and tube pumice parallel to elongation. The misfit with the other pumice types can have numerous causes as the empirical model contains many assumptions. Further refinement of the model can be achieved, e.g., by adding the influence of channel shape (Mortensen et al., 2005), or replacing Archie's law (Matyka et al., 2008). However, this will reduce its applicability within conduit flow models as more free parameters are introduced. The current model succeeds in separating different types of pumices as well as detecting anisotropic features within pumices, providing a reasonable way to estimate permeability, while keeping a necessary degree of simplicity.


Obtaining rock permeability data is fundamental in many areas of geosciences. 3D quantification of a porous medium is an absolute necessity to understand flow transport. As one might be able to estimate geometrical quantities such as specific surface area and porosity from 2D data, it is evident that important topological quantities of the medium like connectivity and tortuosity cannot be determined in such a way. As demonstrated by Figures 2 and 3 (and Animations 1 and 2), gas transport is governed by only a few dominant connected channels. Therefore only a small part of the porosity will govern gas loss in pumices. Moreover, the data in Figure 7 suggest an important control of tortuosity on permeability. X-ray microtomography imaging and LBM numerical simulations are becoming robust and accessible tools to deal with this issue as problems with spatial resolution and REV become negligible. The availability of open-source software such as 3dma-rock, Blob3d, ImageJ, Paraview, and Palabos will allow this method to become widely used and combine with lab measurements to build large permeability databases.

In order to obtain a realistic description of the development of permeability in the volcanic conduit, a thorough understanding of the controlling parameters, i.e., connected porosity, tortuosity, and specific surface area (or a similar length scale), is necessary, and it is clear that these properties are determined by complex dynamics of bubble growth, deformation, and coalescence (Gonnermann and Manga, 2007). Questions on (a) how the connected porosity is related to the total porosity (e.g., Burgisser and Gardner, 2004; Takeuchi et al., 2009), (b) how the tortuosity is related to deformation (e.g., Scholes et al., 2007), and (c) how dominant connected channels develop, need further investigation.

This work would not have been possible without the excellent assistance at the various synchrotron beamlines. We are very grateful to Lucia Mancini and Nicola Sodini (SYRMEP, Elettra), Christoph Hintermüller (TOMCAT, SLS), and Timm Weitkamp (ID-19, ESRF). We would also like to thank Margherita Polacci for supporting this project and assisting at the SYRMEP beamline. We counted on the expertise of Jonas Latt on using Palabos. For the use of the computer clusters we thank Josef Dufek (Atlas), Bastien Chopard and Christian Clemançon (Cadmos), and Nicolas Mayocourt (Myrinet). KPT samples and lab data were provided by Caroline Bouvet de Maisonneuve, CAP and GNP samples by Jean-Louis Bourdier, and the MT sample by Caroline Martel. W.D. was supported by the Swiss NSF grant #200021-111709/1. Work partially funded by the ERC grant 202844 under the European FP7. Reviews by Don Baker and two anonymous reviewers improved the manuscript.

Supplementary data