X-ray tomography is a nondestructive technique that can be used to study rocks and other materials in three dimensions over a wide range of sizes. Samples that range from decimeters to micrometers in size can be analyzed, and micrometer- to centimeter-sized crystals, vesicles, and other particles can be identified and quantified. In many applications, quantification of a large spectrum of sizes is important, but this cannot be easily accomplished using a single tomogram due to a common trade-off between sample size and image resolution. This problem can be circumvented by combining tomograms acquired for a single sample at a variety of resolutions. We have successfully applied this method to obtain crystal size distributions (CSDs) for magnetite, pyroxene + biotite, and quartz + feldspar in Bishop Tuff pumice. Five cylinders of systematically varying size (1–10 mm diameter and height) were analyzed from each of five pumice clasts. Cylinder size is inversely proportional to image resolution, such that resolution ranges from 2.5 to 17 μm/voxel with increasing sample size. This allows quantification of crystals 10–1000 μm in size. We obtained CSDs for each phase in each sample by combining information from all resolutions, each size bin containing data from the resolution that best characterizes crystals of that size. CSDs for magnetite and pyroxene + biotite in late-erupted Bishop pumice obtained using this method are fractal, but do not seem to result from crystal fragmentation. CSDs for quartz + feldspar reveal a population of abundant crystals <35 μm in size, and a population of crystals >50 μm in size, which will be the focus of a separate publication.


X-ray tomography is a powerful technique for qualitative and quantitative studies of rocks and other materials in three dimensions (for recent reviews, see Stock, 1999, 2008, 2009; Jerram and Higgins, 2007; Jerram et al., 2009). One of the key advantages of X-ray tomography is that it is relatively nondestructive, which allows for objects to be studied in situ. X-ray tomography results in a three-dimensional map of linear attenuation coefficient (herein referred to as a tomogram) that allows direct observation and measurement of properties of objects in three dimensions, entirely avoiding stereological corrections needed for measurements made in two dimensions (see Ketcham and Carlson, 2001; Gualda and Rivers, 2006; Jerram and Higgins, 2007).

Image resolution is, in principle, only limited by the characteristics of the radiation and instrumentation used. In practice, however, the image resolution that is attainable is often constrained by sample size, due to limitations on the size of area detectors or available instrument time (see Stock, 2009). This is a significant limitation, particularly for size distribution studies, which require quantitative information on objects that span a wide range of sizes (e.g., Jerram et al., 2009). For example, in this study we are interested in quantifying crystals as large as 1 mm, requiring large volumes (cubic centimeters) to be analyzed. At the same time, we are interested in populations of small crystals (∼10–20 μm), which require analysis at high resolution (∼2 μm or better). However, given the current constraints on technology and time, we are unable to analyze cubic centimeter volumes at micrometer resolution, making it impractical to obtain quantitative information for the desired range of object sizes using a single tomogram.

We present here a simple method to reduce this problem, which makes use of tomograms obtained at various resolutions (2.5–17 μm per volume element, voxel), to allow study of objects spanning approximately two orders of magnitude in size (10–1000 μm). This is a stepwise approach that allows us to obtain information for large crystals from large samples imaged at low resolution and for small crystals in small samples imaged at high resolution. The method is essentially an extension to three dimensions of nested scanning electron microscope images used in many vesicle size distribution studies (e.g., Adams et al., 2006), but differs in that we must analyze several pieces from each sample due to the technological limitations on generating high-resolution tomograms within a single large piece.

In Gualda and Rivers (2006), X-ray tomography was employed to yield quantitative textural information, particularly crystal size distributions (CSD), for pumice clasts of the Bishop Tuff (California) (Bailey et al., 1976; Hildreth, 1979, 2004; Wilson and Hildreth, 1997, among many others). In that study the tomograms used were of relatively low resolution (8.5–17 μm/voxel), such that populations of small crystals (<35 μm) could not be reliably quantified. We are currently detailing textures of late-erupted Bishop Tuff pumice and using a significantly wider range of tomogram resolutions (2.5–17 μm/voxel), which has allowed us to document and quantify a petrologically important population of small crystals (10–35 μm). This previously undocumented population of quartz and feldspar crystals in Bishop pumice will be the subject of a separate publication.


A variety of tomographic imaging systems is commercially available and can be used to image objects of geological interest, including pumice and other rock types (Stock, 1999, 2008, 2009). Image resolution achievable with X-ray tomography has increased significantly in recent years and continues to grow. At present, synchrotron-based systems can routinely yield resolutions ∼1–20 μm/voxel, while desktop systems can achieve a wider range of resolutions (e.g., ∼1–250 μm).

Resolution of tomographic images is typically limited by the volume of the sample, either due to limitations on the physical size of area detectors or by constraints on available analytical time. Consequently, although some desktop systems can be used to analyze very large samples, they tend to result in relatively low resolution data. In turn, synchrotron-based X-ray tomography yields relatively high resolution data of great quality, but places significant constraints on the maximum sample size (see Rivers et al., 1999; Ketcham and Carlson, 2001; Sutton et al., 2002; Stock, 1999, 2008, 2009).

A key goal of this study was to obtain data for a wide range of crystal sizes while keeping uncertainties under reasonable levels. Minimum errors in CSDs are due to counting uncertainties, which follow Poisson statistics (Gualda, 2006). As such, while increasing the total number of crystals counted from 100 to 1000 decreases the relative uncertainty from 10% to 3%, a further increase from 1000–10,000 crystals only decreases uncertainties to 1% relative. These small counting uncertainties will easily be overwhelmed by other sources of error, such that there is little to be gained from counting more than 1000 crystals in a given size bin. In addition, image processing becomes much more challenging when more than a few thousand crystals are involved. Mock and Jerram (2005) and Gualda (2006) showed that CSDs with at least 300 objects adequately describe CSDs, such that the ideal scenario is one in which the total number of objects being counted is on the order of a few thousand, but with no size bin including fewer than ∼100 objects (10% relative uncertainty). Given that rocks tend to include an increasing number of crystals and vesicles with decreasing object size (see Marsh, 1988, 1998), a strategy is needed to manage the number of objects in the various size bins and achieve the compromise sought between data quality and feasibility of imaging and processing.

We give an example to make the problem more apparent. In our work, a high-resolution (2.5 μm/voxel) tomogram of a small cylinder (∼1 mm diameter) yielded 84 crystals 8.5–17 μm in size. If a large cylinder (∼1 cm diameter) could be analyzed at this same resolution, it would yield ∼84,000 crystals in the 8.5–17 μm size range. A more extreme example is provided by another ∼1 mm cylinder analyzed at high resolution, which yielded ∼1500 crystals 8.5–17 μm in size; a 1 cm cylinder of this sample analyzed at high resolution would contain 1.5 × 106 crystals!

If it were practical to image a 1 cm sample at 1–2 μm/voxel, subvolumes could be used to manage the total number of objects per bin size. However, we are limited by the characteristics of the imaging systems or by analytical time, such that this approach is currently impractical. We suggest here an alternative approach: We analyze various pieces from the same sample, such that a relatively small sample is imaged at relatively high resolution to yield information on smaller crystals; a larger sample is imaged at somewhat lower resolution but is only used to yield data on crystals that appear in sufficient numbers and that are large enough to be adequately imaged at that resolution. As a result, we are able to keep the total number of objects per bin size in the desired 100–1000 range.

Here we illustrate the application of this method to determine CSDs spanning a wide range of sizes. The method can, however, be used for any number of applications. For example, detailed quantitative and qualitative examination of other textural features, such as spatial distributions and morphologies of objects at different scales, is possible with this method. In addition, while we focus here on pumice clasts, which are particularly well suited for tomographic studies of crystal and vesicle textures, the method is applicable to a variety of other rock types and materials.


In this study we analyzed five pumice clasts from late-erupted Bishop Tuff. Samples were collected from the locality known as the Aeolian Buttes (Universal Transverse Mercator grid 17.4E/92.1N, section 14 of T1S, R26E), which are part of the northern deposits of the Bishop Tuff (Ig2NWb; Wilson and Hildreth, 1997).

Sample Preparation

Five to seven cylinders were cut from each pumice clast. Sizes of the cylinders varied systematically. Size was reduced by twofold steps in diameter; while tomographic image resolution varies linearly with sample diameter, sample volume is proportional to the cube of the diameter, resulting in eightfold steps in volume. This range of sizes corresponded with a range in tomogram resolution from 2.5 to 17 μm/voxel (Fig. 1). In earlier studies (see Gualda and Rivers, 2006), it was demonstrated that objects larger than 4 voxels in linear dimension (>35 voxels per object) can be reliably quantified in any particular tomogram. Given that we previously used only 2 cylinder sizes (A, B) with resolutions of 17 μm/voxel and 8.5 μm/voxel, respectively, we were able to quantify crystals down to 35 μm. The use of 2 more resolutions in this study has reduced the quantifiable size to ∼10 μm.



All tomographic imaging was performed on the bending magnet beamline of the GeoSoilEnvironCARS (GSECARS) sector of the Advanced Photon Source at Argonne National Laboratory. Detailed descriptions and illustrations of the experimental setup at GSECARS were given in Rivers et al. (1999), Sutton et al. (2002), and Gualda and Rivers (2006).

Cylinders were placed on a stage and rotated 180° in 0.25° steps. At each step, a radiograph of the sample was obtained, resulting in 720 radiographs per tomogram. Reconstruction of this radiograph sequence was performed using the program ‘tomo_display’ (Rivers and Gualda, 2009), and results in a three-dimensional (3D) map of linear attenuation coefficient (for details, see Rivers et al., 1999).


Image processing was performed using ‘Blob3D’ (Ketcham, 2005) and ‘vol_tools’ (Rivers and Gualda, 2009), following procedures detailed in Gualda and Rivers (2006), with two major adjustments.

1. In this work, we did not go through the time-consuming process of separating touching grains. Pumice is characterized by isolated crystals in a lower attenuation glassy matrix, such that touching crystals are very rare and can be easily identified in the output of Blob3D.

2. In each clast, we analyzed for three mineral components: magnetite, pyroxene + biotite, and quartz + feldspar. Distinctions within these groups, particularly between quartz, sanidine, and plagioclase, are challenging and very time consuming (see Gualda and Rivers, 2006). Quartz and pyroxene are much more abundant in Bishop pumice than feldspars and biotite, respectively, and our results for early-erupted Bishop pumice (Gualda and Rivers, 2006) demonstrate that CSDs for sanidine are similar in shape to quartz CSDs. Therefore, we opted to treat quartz + feldspar and pyroxene + biotite in combination. Caution is necessary when assessing CSDs of multiple phases in combination because independent CSDs of phases may show signals of different petrologic processes. If a CSD is determined for these phases in combination, indicators of these processes may be lost due to overprinting or in the combination of different CSD shapes. We believe that the considerable abundance of quartz and pyroxene relative to feldspar and biotite makes this problem relatively minor in this case.

These two adjustments allowed us to study a much larger number of samples and tomograms per sample than possible in our previous work (Gualda and Rivers, 2006), yielding quantitative results over a significantly wider crystal size spectrum.



The use of tomograms obtained at various resolutions allows us to examine the textures of pumice over a range of spatial scales. Tomograms from each resolution show qualitative and quantifiable differences in textures of glass and vesicles, as well as the extent to which fragmented crystals exist at different scales (Fig. 1).

For example, the low-resolution tomogram shown in Figure 1C (Run A: 17 μm/voxel) contains a large fragmented quartz crystal. The high-resolution tomogram (Run D: 2.5 μm/voxel), however, shows the texture of the surrounding glass in great detail, permitting visualization of the fine vesicularity typical of Bishop pumice.

In addition, 3D reconstructions of classified phases can be used to assess crystal and bubble morphologies, as well as spatial distributions of crystals at various scales (Fig. 1D). For example, at the highest resolution (Run D), there are very few crystals, but the morphology of the magnetite crystals is clear. Furthermore, the fragmented crystal seen in the tomogram at low resolution is clearly apparent in the animation of this cylinder. The difference in crystal size and spatial distribution from the high-resolution to the low-resolution runs is also clear in these reconstructions.


CSDs are typically shown as semi-logarithmic plots of crystal size versus population density (Marsh, 1988, 1998). In this work we have used the maximum crystal dimension (in μm) as the measure of crystal size, which was estimated by the maximum dimension of a best-fit ellipsoid determined by Blob3D for each classified crystal.

The population density in a given bin size is the number of crystals present/bin size/sample mass (in units of μm−1 mg−1). We use mass for normalization instead of volume (as is often done in CSD studies) due to the high and variable vesicularity of pumice. As such, sample volume is not as meaningful as in holocrystalline rocks and is also not as easily determined. Glass volume has also been suggested as a normalizing quantity (Proussevitch et al., 2007), but is difficult to determine when vesicularity is at a scale smaller than tomogram resolution. Mass determination is trivial, however, and can be used if the whole sample is imaged (i.e., the whole sample fits in the tomogram; see also Gualda, 2006).

Uncertainties on all plots are based on counting statistics and are equal to forumla where N is the number of crystals counted (Gualda, 2006).

Combination of Information at Various Resolutions

For each sample, at least one cylinder was imaged at each of four different resolutions (17, 8.5, 4.25, 2.5 μm/voxel; see Figs. 1A. 1B), with cylinder diameters and heights of ∼10, 5, 2.5, and 1 mm, respectively. Two different cylinders were analyzed at the lowest resolution (Run A), and one cylinder was imaged for each of the three higher resolutions (Runs B–D). Each cylinder was labeled uniquely, with a letter and a number and referred to as a “run.” The letter corresponds to the resolution of the resulting tomogram and the number label denotes multiple runs at that resolution (e.g., the label AB_5301_RunA indicates pumice clast AB_5301, 17 μm/voxel, first cylinder imaged at this resolution, while AB_5301_RunA1 corresponds to the second cylinder analyzed at the same resolution).

Accordingly, the CSD obtained from each cylinder provides useful information in different size ranges. Combining all five runs from one pumice clast thereby allows us to obtain CSDs over a large range of sizes, spanning nearly two orders of magnitude (∼10–1000 μm), a significant improvement over previous studies (e.g., Gualda and Rivers, 2006).

The largest cylinders (Runs A and A1) were analyzed at 17 μm/voxel. In these runs, a drop in the abundance of crystals smaller than ∼60 μm (4 voxels in linear dimension) is observed, while in higher resolution runs (Runs B–D), higher numbers of crystals are found in this size range (Figs. 2A, 2C). It is clear that the low-resolution runs cannot adequately quantify crystals smaller than 60 μm (see Gualda and Rivers, 2006). In contrast, high-resolution runs are adequate for small crystals, but they lack large crystals (Figs. 2A, 2C). In this case, the imaged sample is small, such that there is a lower probability, if not an impossibility, of the cylinders including crystals of these larger sizes.

In all runs, a drop in the number of crystals smaller than 4 voxels in diameter is apparent (Figs. 2A, 2C). This demonstrates that, regardless of the resolution, a minimum number of ∼30 voxels is necessary for adequate quantification of crystals in a tomogram (Gualda and Rivers, 2006). In our CSD plots, we designated bin sizes such that each bin boundary is equal to an integer number of voxels (Fig. 2). Using this system, this quantification limit is easily recognized for the various runs.

It is interesting that, in a few cases, a higher population density is observed in a higher resolution tomogram (smaller analyzed volume) than in lower resolution tomograms (larger analyzed volume). We attribute this to the so-called “nugget effect” commonly observed in the mining industry, whereby rare nuggets of the ore of interest may strongly bias quantitative estimates of the bulk abundance of the ore in the deposit (U.S. Bureau of Mines, 1996). In our case, the occasional presence of a small number of unusually large crystals or vesicles can lead to apparent population densities that are overestimated. Compare, for example, the 140–280 μm bin range for Run C in magnetite with Runs A and B (Fig. 2A); what may not be immediately apparent from the way the data are plotted is that both Run C and Run B include only 2 crystals in that size bin, while Run A includes 43 crystals. In order to minimize this effect, for each bin size we use information from the largest cylinder that can appropriately quantify that particular bin size, such that larger and more representative sample volumes are used when possible (see Fig. 2).

In order to obtain CSDs spanning a large range of sizes, the results obtained for the four independent resolutions are combined into a single CSD for a given sample (Figs. 2B, 2D). Data from multiple runs analyzed at the same resolution (e.g., A and A1) can be simply combined by adding the number of crystals and dividing by the total mass analyzed when obtaining population densities. However, a combination of data from different resolutions is impractical. Because the value of the population density is normalized to the mass of the analyzed cylinder, the results obtained can be integrated seamlessly. Therefore, for each bin, data are taken from the run that most confidently resolves crystals of that size (Figs. 2B, 2D): the largest bin sizes (70–140 μm, 140–280 μm, 280–560 μm, and 560–1140 μm) come from A runs, while the smaller bin sizes (35–70 μm, 17.5–35 μm, and 8.75–17.5 μm) are taken from B, C, and D runs, respectively.


The CSDs we obtained using this method reveal important aspects of the crystallization and evolution of the Bishop magma. Magnetite and pyroxene + biotite size distributions show a power law, or fractal, distribution (Fig. 3), as evidenced by the high correlation coefficients observed (0.90–0.99; see Figs. 3B, 3D). The fractal dimensions (D) vary over a limited range (D = 2.3–3.0 for magnetite, Fig. 3B; D = 2.5–3.2 for pyroxene + biotite, Fig. 3D). Fractal size distributions like these are often taken to be the result of crystal fragmentation (Bindeman, 2005). However, magnetite, pyroxene, and biotite crystals all appear euhedral in our tomograms (Fig. 1), showing no signs of extensive fragmentation, leading to the inference that these power law distributions were generated by another process. Baker et al. (2006) discussed a number of mechanisms by which a power law distribution may arise; while some of these theories are relevant only for vesicle size distributions (e.g., bubble coalescence; Gaonac'h et al., 1996), others have implications for crystals (e.g., continuous nucleation; Blower et al., 2002). However, Baker et al. (2006) concluded that power law bubble size distributions likely result from bubble coalescence; plus, they demonstrate experimentally that bubble size distributions evolve toward exponential distributions. It is interesting that nonfragmental CSDs are almost invariably exponential, and Marsh (1998) argued that, in the absence of coalescence, exponential distributions are to be expected from kinetic theory. In this sense, magnetite and pyroxene + biotite fractal size distributions obtained here are unexpected and difficult to explain based on theoretical considerations and experimental evidence.

The quartz + feldspar distributions are particularly interesting, as they indicate the existence of two populations of crystals: a population of large (>50 μm) crystals and a population of more abundant small (<35 μm) crystals (Fig. 2). This conclusion is robust regardless of the exact choice of tomogram from which to extract information for a particular bin size; in other words, if data for the range 140–1140 μm were taken from Run A, 70–140 from Run B, 35–70 from Run C, and 8.75–35 from Run D, the CSD would be similar, and the existence of these two crystal populations would still be apparent. We interpret this result to demonstrate that a major event took place that changed the crystallization regime from growth dominated (continually growing few existing crystals) to nucleation dominated (nucleation of many new crystals). In addition, this event likely took place near to the time of eruption, as the many small crystals in our distributions indicate there was little time for crystals to grow before eruption. Detailed discussion of these results and their significance for the evolution of the Bishop magma are beyond the scope of this presentation and will be the focus of a separate publication.


X-ray tomography is an effective method to study crystals and vesicles in rocks (and particles in materials in general). In many cases, it is important to be able to document sizes of particles spanning more than one order of magnitude. In studies of volcanic rocks, for example, small crystals (i.e., <50 μm) may record the ascent and eruption history, while large crystals (>50 μm) record earlier stages of crystallization, such that it becomes necessary to quantify textures over at least 2 orders of magnitude in crystal size (e.g., ∼10–1000 μm) in order to properly understand the record preserved in crystal populations.

A significant barrier in obtaining data over such a wide range of crystal sizes is a usual trade-off between sample size and tomogram resolution, which makes it impractical to adequately quantify objects of a large range of sizes from a single tomogram. This limitation can be overcome by combining data from various cylinders for each sample, imaged at multiple resolutions (i.e., larger pieces, lower resolution; smaller pieces, higher resolution). We successfully employed such an approach in our studies of CSDs in pumice from the Bishop Tuff.

We use 5 cylinders that differ from each other by a factor of 2 in diameter (∼1–10 mm) and resulting tomogram resolution (2.5–17 μm/voxel) and combine data from all cylinders (or runs) to yield CSDs for each sample. For each bin size, we select data from the run that best quantifies crystals of that size; this allows us to obtain an overall crystal size distribution for each sample that spans 2 orders of magnitude in crystal size (10–1000 μm).

Application of this method to pumice clasts from the late-erupted Bishop Tuff reveals the following.

  1. Magnetite and pyroxene + biotite show fractal crystal size distributions, but appear to be euhedral in our tomograms. This suggests that fragmentation, which is often considered the cause of fractal distributions, is unlikely to be the explanation for the observed distributions in this case, and another mechanism needs to be sought.

  2. It is interesting that quartz + feldspar distributions show the existence of two crystal populations: (1) a population of numerous small (<35 μm) crystals, not recognized in previous studies (Gualda et al., 2004; Gualda and Rivers, 2006); and (2) a population of less abundant large (>50 μm) crystals. Recognition of these two populations has important implications to the evolution of the Bishop magma, which will be the focus of future publications.

We emphasize that the method we suggest is not limited to studies of crystals in pumice; this method can be used in almost any quantitative textural study of rocks, as well as for many other kinds of materials. It is ideally suited for cases in which the number of objects is inversely correlated with their size, a common feature in natural systems.

We are indebted to Mark Rivers for his continued guidance, support, and interest in our work at GSECARS (GeoSoilEnviroCARS, Argonne National Laboratory), and for many thought-provoking conversations and suggestions. We also thank Alfred Anderson for his thoughts and discussions about our work on the Bishop Tuff. Portions of this work were performed at GSECARS (Sector 13), Advanced Photon Source (APS), Argonne National Laboratory. GSECARS is supported by the National Science Foundation-Earth Sciences (grant EAR-0622171) and Department of Energy-Geosciences (DE-FG02-94ER14466). Use of the APS was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract DE-AC02-06CH11357.

Supplementary data