We have developed a novel method of measuring bubble size distributions from their remains expressed on the surfaces of volcanic ash particles. The morphology of the ash fragments retains a record of bubble size at the time of fragmentation in the curvature of the convex surfaces on the ash fragments. This curvature can be measured using stereo scanning electron microscopy (SSEM), and morphology can be represented using a digital elevation model (DEM) of ash particle surfaces. Due to the vagaries of the sensitivity of SSEM imagery to surface roughness, a three-point fit technique produces more robust results for curvature than a least-squares approach for curve-fitting of ellipsoids of revolution to ash surfaces. The method allows measurement of vesicles within a size range from one to over a million cubic microns. The inferred bubble size distributions so obtained can potentially provide valuable insights regarding magma dynamics and vesiculation that lead to the explosive eruptions that produce ash. An error analysis of the methodology indicates reasonably accurate reconstruction of bubble geometries and bubble size distributions (BSD). Accuracy is constrained primarily by the size of ash particles themselves since the mode of the BSD should be at least a standard deviation smaller than the dominant dimensions of the particles for robust results.


Explosive volcanic eruptions pose serious natural hazards, but erupted materials can provide insights regarding the processes that drive and control the nature of such energetic eruptions. Plinian eruptions are the most energetic (and thus devastating) class of eruption, and the rapid expansion of volatile-rich magmatic foam leads to fragmentation of the magma into fine particles. In such events, most of the erupted material fragments into particles that are smaller than 2 mm, and 30% to over 50% of them are fine ashes less than 100 μm in size (Rose and Durant, 2009). Because the growth of bubbles of gas (mostly water near the vent) by exsolution and decompressive expansion drives eruptions (Proussevitch and Sahagian, 1998), studies of the vesicles preserved in volcanic products can shed light on eruption processes. However, most of the material erupted is in the form of fine ash, which has never yet been analyzed to provide information regarding the bubbles that burst to create the ash.

The goal of this study is to develop a technique to measure bubble sizes from the fragments of vesicles reflected in volcanic ash particles primarily for the purpose of reconstructing bubble size distributions (BSD, see acronyms in Table 1) in the ash source magmas at the time of their fragmentation during eruption. When these BSDs are determined, they can be used for characterization of the vesiculation parameters of the host melts (bubble nucleation density, growth rates, growth factors such as diffusion and/or decompression, coalescence, etc.) to better understand the driving mechanisms of volcanic eruptions.1

Volcanic eruption products of all sizes are affected by bubbles of gas that create vesicles. In large contiguous particles (pumice clasts, bombs, etc.) many complete vesicles can be observed and measured. Numerous vesicles can also be observed on the surfaces and interiors of large ash particles (Riley et al., 2003). However, with smaller ash particles that are on the order of the size of the vesicles (5–10 μm), it is less clear if a concave surface reflects part of a broken vesicle, or if it reflects a concoidal fracture of glass. As such, the lower limit of ash particle size that can be practically used in this analysis is ∼5–10 μm.

The opportunity thus exists to measure the curvature of these craters to reconstruct their prefragmentation bubble sizes. These craters are believed to faithfully represent original bubble wall curvature because time scales of ash cooling in air (minutes, according to Thomas and Sparks, 1992) are much smaller than the surface tension strain relaxation (days, as per Rust and Manga, 2002) in the postfragmentation melt phase. When an unbiased and statistically representative population of vesicle sizes could be collected from multiple ash particles in a given sample, it is possible to reconstruct their bubble size distribution. This BSD can represent a preexisting magmatic foam status at the time it has undergone brittle fragmentation during explosive eruptions (Spieler et al., 2004).

In this paper we describe a new measurement technique, based on commonly available instruments and software that can be readily reproduced in any research laboratory. The BSDs obtained through this technique can then be used in many volcanology applications that address volcanic product characterization and physical volcanology problems related to magma vesiculation and brittle fragmentation, such as tensile strength of silicic melts above glassification temperatures and fragmentation topology that is relevant to ash particle sizes. The ability to determine prefragmentation bubble size distributions from ash samples may have additional applications beyond eruption dynamics, and these are left for future studies and the broader research community.

Technical Aspects

The technique is based on stereo-pair imaging of ash fragments by scanning electron microscopy (SEM) with sample tilting capabilities, which we refer to as stereo SEM or SSEM. Traditionally, SEMs have been used to acquire high contrast images of micron and submicron samples. These images have been used to determine two-dimensional (2D) parameters of the objects of interest. With recent advances in image analysis, it is now possible to determine the three-dimensional (3D) parameters from stereoscopic images captured using a SSEM as illustrated in Figure 1. The parallax inherent in the stereoscopic images is then used to create a digital elevation model (DEM), or “topography,” of a sample particle's surface. Subsequently, the DEM can be analyzed to determine the 3D parameters of the object, such as curvature of vesicle wall surfaces (Stampfl et al., 1996).

In this study we used a JSM-6400 scanning electron microscope and S-4700 field emission scanning electron microscope at the Michigan Technological University (courtesy of Professor W. Rose). Stereo pair images were processed at the University of New Hampshire using MeX, a commercial analysis package (Alicona Imaging GmbH, Graz, Austria).

Volcanological Applications

Bubble size distributions (BSD) in fine volcanic ash determined from SSEM can provide insights regarding the physical conditions of the vesiculation processes in explosive eruptions (Blower et al., 2003; Cashman and Mangan, 1994; Gaonac'h et al., 2005; Proussevitch et al., 2007a; Toramaru, 1990). This is important as it can shed light on the extreme vesiculation conditions that lead to brittle fragmentation of the magmas that produce observed ashes (Alidibirov and Dingwell, 2000; Gardner et al., 1996; Koyaguchi and Mitani, 2005; Mader, 1998; Papale, 1999; Sahagian, 1999; Spieler et al., 2004; Zhang, 1999).

We consider two classes of ash particles that result from magma fragmentation during eruption. The first is a “simple” particle that represents the plateau border between three or more bubbles and contains no complete or partial vesicles within it. The full length of each side of its exterior is an arc that represents part of a bubble wall. By determining the radius of curvature of the arc in multiple dimensions, the volume of the original bubble can be determined. Actual reconstruction of a complete vesicle from “matching” ash particles like a jigsaw puzzle is impossible due to the broad areal scatter of particles in an eruption, so we assume that each particle we measure is the sole representative of its adjacent bubbles. The second class of ash particles is larger, and contains numerous indentations, or craters, on its exterior surfaces, each representing part of a bubble. These “compound” ash particles may form from the fragmentation of several surrounding large bubbles, and may even contain much smaller, complete vesicles in their interiors. Although we cannot measure internal vesicles with SSEM (but this is possible with XuM; Kiely et al., 2007), these compound particles are ideally suited for populating size distributions with the numerous vesicles reflected on their surfaces. We consider both types of particles in our analysis.

Ash samples for the study were provided by the Alaska Volcano Observatory, U.S. Geological Survey (courtesy of Drs. J. Larsen and K. Wallace). We processed and analyzed fine ashes from Hayes volcano eruption (1600 B.C.), Katmai eruption (1912 A.D.), and Augustine eruption (2006 A.D.). All three volcanoes are located relatively close to each other in South Central Alaska just north of the base of Aleutian arch.

Hayes volcano is located in the Tordrillo Mountains between the Alaska Range and Cook Inlet of the Gulf of Alaska. The samples examined in this study are from the last major eruption 4400–3600 yr B.P. (Begét et al., 1991; Riehle et al., 1990). Waythomas and Miller (2002) suggest that the event “was probably a Plinian-style eruption characterized by multiple, explosively generated ash clouds that likely extended high into the atmosphere, possibly reaching stratospheric levels.” The eruptions generated widespread tephras with a distinctive biotite, amphibole, and pyroxene mineralogy (Riehle et al., 1990).

The Katmai samples are from the eruption of 1912 near Mount Katmai that resulted in the most voluminous volcanic event of the twentieth century. The erupted magma ranged from andesitic to rhyolitic compositions (Hildreth and Fierstein, 2000).

The Augustine ash samples are from the eruption event in 2005–2006. There were 13 explosive episodes in 20 days that resulted in up to 14 km ash columns (Power et al., 2006). The erupted material has predominantly andesitic composition.

While the ashes from these eruptions have a wide range of particle sizes that mostly depend on the collection distance from the volcanoes, we used samples of fine ashes with dominant particle sizes from 10 to 20 to 50–60 μm since those are too small to study with other known techniques (X-ray tomography, for example, see Polacci et al., 2006) and large enough to have multiple bubble imprints (vesicles) on a surface of a single ash particle.


In this section we describe a step-by-step procedure developed for bubble size measurement from their imprints on volcanic ash particles.

SSEM Imaging (Step No. 1)

The first step involves creating stereo pair images using a SSEM with eucentric tilting (Piazzesi, 1973) (Fig. 1). The number of ash particles and, therefore, stereo pairs, depends on the total number of bubbles (partial vesicles that define ash surfaces) needed to be measured in order to produce a statistically representative bubble population for each ash sample (at least 100 bubbles for populations with the ranges of sizes observed to date—see error analysis below). Measurement of more vesicles for statistically representative distributions would benefit from measurement automation software, but this is beyond the scope of this paper. We observed and created stereo pairs from SSEM of ∼10 ash particles, each containing ∼5–50 vesicles on their surfaces (Fig. 2). Every vesicle on each ash particle was characterized to ensure an unbiased sample.

Building Digital Elevation Model (DEM) from Stereo Pairs (Step No. 2)

Each stereo pair of images was processed to build a DEM of the ash particle and create a surface topography map (see Fig. 1). In this step, we used the standard commercial Alicona MeX software package, which often comes bundled with tilting SEM hardware. It is necessary to manually input the working distance, pixel size, and the angle of tilt for each stereo pair to generate the DEM. The time taken to generate the DEMs is anywhere from a few seconds to a few minutes depending on the video processing capabilities of the host computer. The number of particles and their DEM maps per ash sample depends on the number of vesicles needed for a statistically representative population for BSD analysis.

Algorithms to extract depth detail from stereo images are commonly classified into two categories—area-based matching and feature-based matching (Barnard and Fischler, 1982). Alicona MeX uses an area-based algorithm by Scherer et al. (1999) that in turn employs a rank-based method to match homologous (also known as characteristic or reference) points (i.e., surface pixels representing the same location in each image) and calculates their z values. The 3D point cloud is networked using triangle mesh, based on Delaunay triangulation, to form the DEM.

It is important to note that DEM reconstruction works on homologous points only so that the z coordinates of these points are directly calculated from stereo-pair (x,y) coordinates. Normally there are forumla homologous points in an image of an irregular shaped object, where Npixels is the number of pixels in the image. So the z values for all other pixels are interpolated. This impacts the quality of the DEM reconstruction the most in those parts of the image where shading gradient (curvature) areas are lacking reference points (scratches, nicks, dust, etc.). This is quantified in the Error Analysis section below.

Building Bubble Cross Sections (Step No. 3)

The DEM of the ash particle surfaces enables the construction of cross-section profiles of the bubble wall fragments (craters) preserved as partial vesicles on ash particle surfaces (Fig. 2). These profiles are essential for subsequent calculations of vesicle curvature and eventually for reproduction of bubble sizes. There is no automated image processing software for vesicle size measurement from ash available as of yet, so it is quite laborious to measure hundreds of vesicles in multiple volcanic ash particles. The manual part of the technique (step No. 3) poses a major challenge in obtaining a statistically representative and unbiased population for BSD analysis.

Because the technique measures curvature of concave bubble imprints on ash particle surfaces, it is up to the discretion of each investigator to choose the optimal number and orientation of curvature profiles that would best characterize the geometry of each individual vesicle and facilitate calculation of its original volume. For instance, we used two orthogonal profiles across craters of bubbles that appear to have had spherical or elliptical shapes in order to calculate their volumes. The number of profiles needed depends on the shape of the crater (one for spherical; two for ellipsoid of revolution; more for irregular shapes). To make the profiles, we manually connected two opposite points on the crater perimeter so that each line goes through the deepest part of the crater center (Fig. 2).

It is very important to minimize bias in vesicle selection for the curvature profiles. All visible craters on an ash particle should be profiled. On the processed images we marked each measured vesicle on the exposed ash particle to ensure that none are missed or redundantly measured. This is the most laborious step of the process because curvature profiles must be taken manually. It may be possible to automate this step in the future. We have developed a Graphical User Interface (GUI) application called BubbleMaker for building profiles, circle or elliptical fit evaluation (see Step No. 4), and volume calculations (see Step No. 5). BubbleMaker creates session files with graphical illustration of the cross-sectional profile locations on an ash particle, and saves all profile and calculations data in XML format, so that the work of a technical operator (e.g., student) could be reviewed by his/her supervisor. A deployable application (BubbleMaker) is available from the authors upon request.

2D Circle and Ellipsoid Fit Analysis to Evaluate Curvature (Step No. 4)

In this step of the method, curvature of bubble wall fragments is calculated from their 2D cross-section profiles by ellipsoid functions using two function fit analysis approaches, but one is clearly more robust, as explained below.

The first approach, referenced here as the three-point method, involves fitting a circle to three characteristic points on the profile chosen by the investigator such that 
where R is the circle radius and (xc, zc) are coordinates of its center (Fig. 2).

The three user-selected points usually include two opposite points of the crater edges and a clearly identifiable middle point. There are two advantages of this method. The first is that the user can accurately identify the crater edge and bottom that are image homologous points (see section “Building Digital Elevation Model [DEM] from Stereo Pairs [Step No. 2]”) and thus minimize error in DEM reconstruction. Other points on the profile might be a result of interpolation of the DEM reconstruction and thus have much greater uncertainty. This is because of the way SEM stereo imaging works. The intensity (number) of electrons returning to the detector from the particle surface depends on the material and the angle of the surface. A smooth surface with uniform curvature will not show any rapid changes in reflected electrons as the beam scans across the surface. This results in a uniformly gray image on which there are no identifiable points to match between the two projections of the stereo image. If, however, there are small dust particles on the surface, these can be identified and matched for accurate 3D elevation mapping from two tilted images in a stereo pair. (How can you see the invisible man? Throw dust on him!) If there is a particularly large particle on the surface, the particle will be seen, but the smooth crater walls surrounding it will not be as visible to the process. This leads to a major disadvantage of the least-squares method described below. The second advantage of the three-point method is that it can be done on the DEM map directly and thus bypass the time-consuming profiling step, thus facilitating future automation.

The second curvature calculation approach, referenced as the least-squares method, utilizes all points of the cross-sectional profile for either the circle (Rx = Rz) or ellipsoid (RxRz) function, 
The task involves implementation of a nonlinear least-squares method to the profile by minimizing the geometric error (sum of squared distances from the points on the profile to the fitted circle). This analysis was done along with a chi-square function minimization fit analysis (Johnson and Kotz, 1970) with algorithms of Numerical Recipes (Press et al., 1992) as it was implemented in Proussevitch et al. (2007b). One problem in this approach is Cartesian coordinates of the DEM cross-sectional profiles. On the steep parts of the profiles, chi-square minimization is ambiguous as it does not adequately represent the distance between an observation point and the fitting function. Consequently, we used circle fit in Cartesian coordinates first—as shown by Equation (1)—to find the first-order circle center (xc,zc), which we then used as the center for spherical coordinate system transformation and initial iteration for the circle center in the subsequent function minimization of an ellipse (RxRz) or circle (Rx = Rz) in spherical coordinates 
where α is the angle of inclination from the vertical axis. We thus minimized the radial distance between fit circle and the points of the cross-sectional profile.

We find that the three-point method yields more reliable vesicle sizes than the least-squares method because it is based on DEM homologous points (benchmarks), while most of the profile z values for the least-squares method are the result of DEM interpolation. In essence, the least-squares method is fooled by the presence of large particles (dust) on the crater surface, seeing only the tops of these, and “draws” a curve along the tops to the crater edge. This results in a larger than actual radius for the relatively clean craters, and in a smaller radius for dust contaminated craters as more profile points over the dust particles are located inside the actual vesicle surface (Fig. 3). For this reason, the three-point method is more robust, as the investigator chooses an interior point judiciously and is not “fooled” by the presence of large dust particles on the crater surface. Indeed, it is the presence of small dust particles that enhances SSEM sensitivity (invisible man effect).

Volume Calculations (Step No. 5)

Bubble volume calculations depend on geometry of each individual vesicle. Because the most common (and easily calculated) shapes are spheres or ellipsoids, we normally used two orthogonal cross-sectional profiles for those in order to calculate their volumes from measured curvatures by applying appropriate volume equations (Fig. 2 illustrates an example of a spherical vesicle measured by curvatures of two orthogonal profiles). For some other shapes (e.g., pipe vesicles) volume calculations can involve measurements of axial length in addition to curvature of perimeter in one or more places along its channel. Other more complex shapes, such as hexagonal honeycomb textures of highly vesicular foams, might require curvature measurements in more than two profiles as they depend on shape-specific geometries.


Reconstructed volumes of collected bubble populations in the three processed ash samples were analyzed using statistical approaches formulated in Proussevitch et al. (2007b). The resulting bubble size distributions (BSDs) are presented in Table 2.

All three ash samples displayed monomodal log-normal BSDs with many similar characteristics, such as modal bubble sizes (5–12 μm in diameter) and standard deviations (Fig. 4). Log-normal distributions have been argued previously (Proussevitch et al., 2007a) to reflect simple, single-stage bubble nucleation and growth, so a test of the ability of our method to measure bubble sizes is to compare “observed” sizes calculated from constructed cross sections to a theoretical log-normal distribution. The very close fit to log-normal distributions (see χ2 in Table 2 and histogram match to the distribution curves on Fig. 4) indicates simple vesiculation history (Proussevitch et al., 2007a; Shea et al., 2010).

The bubble number density (BND) is a derived parameter emerging from BSD (Proussevitch et al., 2007b). BND must always be referenced to melt volume rather than bulk volume in order to be relevant to the bubble nucleation process. As indicated in Table 2, the samples we analyzed in this study have BND of at least 10,000 per mm3 up to ∼1 million per mm3 (1015 per m3). These BNDs are comparable to those in silicic pumices (Cashman and Mangan, 1994; Klug and Cashman, 1994; Rust and Cashman, 2004) suggesting that regardless of the total vesicularity (initial dissolved water content), the prefragmentation number of nucleation sites per melt volume is about the same for similar silicic magmas. This is consistent with eruption style being driven by postnucleation processes of bubble growth and expansion (Proussevitch and Sahagian, 2005). Consequently, it appears that it is indeed possible to measure bubble sizes from ash fragments, thus paving the way for future studies that explore eruption processes on the basis of eruption products.


The approach described in this work is not a direct method of bubble size measurement, so it is prone to potential error sources as described below. Here we will assess four of those in order of their potential effect on characterization of prefragmentation vesiculation of host volcanic melts. Of course, the best method of error assessment would be a comparison of ash particle BSDs with those directly observed in volcanic pumices of the same eruption. Unfortunately this approach cannot be used in this work due to difficulties of identifying host pumice formations and recognizing correspondence of nonfragmented materials with their fragmented counterparts. The question is whether the BSD of nonfragmented material could be the same as those that had been fragmented to ash. The answer is that they are likely to be different and thus cannot be used for the error analysis of this method. So we will try to assess indirect error sources of the method stemming from actual known component errors of DEM reconstruction, vesicle geometry, profile and vesicle selection human bias of an operator as this is done manually, and finally standard statistical errors relevant to the size of observed population and function-fit procedures.

Instrument and DEM Induced Errors

SEM imaging is well-known for its ability to produce high-quality images for a wide range of magnifications. DEMs are created from stereo-pairs captured using a SEM equipped with a eucentric tilting stage and using input data such as magnification, relative tilt angle, and the working distance.

Studies by Marinello et al. (2008) and Schroettner et al. (2006) have concluded that instrumental errors produced by improper calibration of the eucentric stage or not using an optimal tilt angle are critical to accurate 3D reconstruction and can go as high as 30% of z value (reference depth), so proper calibration of the instrument is very important and we must assume it is done properly on calibration standards within the facilities we use. Assuming the eucentric stage is properly calibrated, measurements by Marinello et al. (2008) on a specially made sample with a vertical step feature varied by 5% from the reference height.

Vesicle Geometry Errors

The next potential source of error stems from some limitation of the bubble imprints on ash particle surfaces to represent their actual geometry. First, we are using a two-axis rotational ellipsoid model of the vesicles since only two vesicle radii could be actually measured in the imprints—the longest and shortest of the imprint crater. Real vesicles may have a three-axis ellipsoid shape or complex warping due to confined conduit or surface flow dynamics (e.g., pipe vesicles). While there have been no systematic studies of vesicle shapes in actual volcanic rocks, the theoretical study of Rust and Manga (2002) of bubble warping by shear stresses of magma flow suggests a two-axis rotational ellipsoid shape. This has been later checked by studying actual vesicle shapes in three volcanic formations (Major Island, Rock Mesa, and Big Glass Mountain) in Rust et al. (2003). So the two-axis measurement method of this study would be adequate for most volcanic vesicles. In rare cases of nonrotational ellipsoidal shapes the error in bubble volume will be proportional to the difference between volumes of rotational and nonrotational ellipsoids: 
where R1R2R3. So, from Equations (4) we can calculate the nondimensional variance υ of volume measure as 
where the radii refer to mean actual and rotational vesicle geometry in the sample and the variance applies to a uniform distribution of these ratios.

A graphical form of Equation (5) is plotted on Figure 5. It indicates relatively small errors associated with the unlikely case of vesicles with three-axis rotational ellipsoids (Rust and Manga, 2002; Rust et al., 2003) measured as two-axis rotational ellipsoids within a reasonable range of these ratios. These errors are unlikely to exceed 25% of absolute bubble volumes that actually translate to even smaller numbers since all bubble size distribution analysis is based on Log10 transformation of vesicle volume units. These, in turn, are likely to become even smaller in subsequent statistical analysis since these errors are believed to be nonsystematic.

Vesicle Selection Bias

A human error or biased BSD could stem from the choice of vesicles to measure on ash particle surfaces by the operator. While the procedure requires measuring all (no exceptions) visible bubble imprints on the SSEM exposed side of the particle, some of these could be potentially missed by the operator or obscured by dust or other obstructions. In the industrial quality-control procedures these types of errors are commonly evaluated by comparing data from several operators over the same pool of samples. Since we have studied just three ash samples, we have carefully examined each stereo-pair image by the team of co-authors for missing (unmeasured) vesicles in order to minimize or eliminate the bias in vesicle selection for measurements.

Distribution type can also affect the measurement results. Based on the fact that ash particle sizes can conflict with actual vesicle dimensions (vesicle dimensions are the same or larger than the ash particle itself) there are two possible scenarios: (a) incomplete and (b) out of range bubble distributions.

Incomplete distributions are relatively easy to identify on the basis of a lacking mode on their distribution density graph. These could result from ash particles that are too small relative to vesicle sizes so that modal and larger vesicles will not be adequately characterized on concave ash surfaces, or from inadequate instrument magnification rendering it impossible to measure modal or smaller vesicle sizes. The ideal case of recoverable distribution truncation is illustrated on figure 5 of Proussevitch et al. (2007a) so that the truncation does not exceed one standard deviation from the mode of any given continuous distribution. All three samples studies in this work have a good range of nontruncated population (Fig. 4), so that BSD function fit could be confidently performed. In any case, it is essential to reconstruct a complete BSD for the application of determining vesiculation parameters.

Out of range distributions are possible if one or more modes of a multimodal distribution fall out of range of measurement capabilities due to instrumentation limitations or due to large differences between vesicle and ash particle sizes. For instance, vesicles of modal size at or around 3 mm cannot be measured on 10–20 μm ash particles. These modes could be completely lost and missing in the resulting BSDs. It should be possible to assess this by analyzing coarse pyroclastic materials of the same eruption episode, but such a study is beyond the scope of this introductory methods paper.

Standard Statistical Errors

Here we present an evaluation of standard statistical errors associated with analysis of measured vesicle population of size N in an assumed Gaussian distribution type over logarithmic units of volume (Proussevitch et al., 2007a, 2007b). This problem is separated into three independent issues: (a) optimization of histogram bin sizes (depends on observed population size); (b) distribution function fit (Gaussian distribution) by error minimization (truncated population histograms) (fig. 5 of Proussevitch et al. [2007a]); and (c) statistical precision (error) for the distribution mode and variance values (also depends on the measured population size).

Optimization of histogram bin sizes is critical in building a distribution histogram to be used in subsequent function fit analysis. Since observed vesicle populations in ash particles in this study are on the order of 100 per sample, a choice of wide bin sizes may translate to fewer points for function fit, while narrow bin sizes may result in poor bin counts and their Poisson probabilities. For a large number of counts in a bin the standard error of counted events is 
where ni is the number of counted invents in bin i. So the bin count error is (a) very sensitive to a bin size; and (b) this error propagates to chi-square values (precision) of function fit analysis in the next step of data processing. The number of histogram bins nb should be in the order of 
However, this method works only for nontruncated monomodal distributions. In order to optimize histogram bin sizes for an arbitrary mono- or multi-modal distribution, we used the method provided by Shimazaki and Shinomoto (2007) that is based on minimizing the cost function 
where Δ is bin width, and (k,υ) = f(Δ) are mean and variance of the bin count of events ki, 

Applying Equations (8) and (9) in turn minimizes bin count errors (6) as demonstrated in Shimazaki and Shimomoto (2007).

Error minimization function fit is performed using histogram bin count data and is described in detail in Proussevitch et al. (2007a). The function fit error is characterized by a χ2 (chi-square) parameter that is scaled over bin count errors of Equation (6) (see section 9 of Proussevitch et al., 2007a). These errors are always smaller than those of Equation (6) if the fit function is correctly chosen (otherwise the function fitting should be considered invalid).

Statistical precision errors for the distribution mode and sigma values are well-known for each type of distribution function. For the Gaussian function used in this work it is 
where N is the modal population of events to which Equation (10) is applied. Based on data in Table 2 and population size of measured vesicle sizes, the statistical precision is ∼0.1 Log10 units of cubic microns which translates to ∼5%–10% of absolute volume values.

Summary of Error Analysis

There are a few potential error sources involved in the analysis of bubble size distributions described above. While the error sources described above may introduce some potentially significant errors stemming from vesicle size measurement of an individual single vesicle (geometry assumptions), these may not significantly affect the final result for calculation of distribution function moments (coefficients) so long as the errors are not systematic. In the worst-case scenario of having a systematic measurement error (unlikely for this method), the mean of the distribution function is still presented in Log10 volume units that does not exceed 25% of absolute values or 0.6 of Log10 volume units. This is relatively small considering that different classes of volcanic vesicles in different natural environments vary in size and number density by many orders of magnitude (6–9 or more of Log10 vesicle volume units).


This work demonstrates the viability of measuring bubble size distributions (BSDs) from fine volcanic ashes. We hope that our new ability to measure vesicle sizes from their control on ash particle morphology will enable us and others to accurately determine bulk vesicularity just prior to the point of fragmentation.

While the primary goal of this study was to develop and demonstrate a technique for measuring vesicle sizes from their fragments in fine, micron-sized volcanic ash particles, we have managed to collect representative populations of vesicles from a few ash samples in hand (Augustine, Katmai, and Hayes). We found that BSDs in all three ash samples have the following common features.

(1) The edges of craters left on ash particle surfaces from disrupted bubbles in fragmenting magma appear not to be rounded or otherwise relaxed, suggesting that they faithfully represent the sizes and shapes of bubbles just prior to brittle fragmentation in an eruption column.

(2) There are two ways to measure vesicle sizes from ash fragments, namely (a) three-point fit, and (b) least-squares fit, of circles and ellipses in vesicle fragment (crater) profiles. These two methods yield consistently different values for both BSD modes and total number density. The least-squares fit leads to larger vesicles because it is “fooled” by particles on the surfaces of the crater remnants of vesicles, thus making them seem shallower than they really are. The three-point method allows judicious avoidance of the tops of such particles.

(3) Vesicles in the studied volcanic ash samples have monomodal log-normal size distributions.

(4) Modes of bubble sizes inferred from our study lie within a narrow range between ∼5 and 12 μm in diameter (100–1000 μm3).

(5) Bubble number density (BND) as referenced to melt volume is at least 10,000 per mm3 and can reach ∼1 million per mm3 (1015 per m3). These calculated BNDs from the samples analyzed suggest that regardless of initial dissolved water content, the number of nucleation sites per melt volume is comparable for similar silicic magmas. This is consistent with eruption style being driven by postnucleation processes of bubble growth and expansion (Proussevitch and Sahagian, 2005).

(6) The method described in this paper provides a new tool for exploring bubble size distributions and preeruptive conditions in energetic eruptions on the basis of the resulting ash fragments, and makes it possible to address a suite of previously intractable issues in subsequent studies.

The ability to determine BSD at the point of fragmentation sets the stage to answer a number of important questions pertaining to eruption dynamics.

(1) What are the vesiculation dynamics that lead to explosive silicic eruptions?

(2) What is the relation between ash particle size distribution and bubble number density?

(3) How do particle size distribution, bubble size distribution, and bubble number density vary with distance from the vent?

(4) Do large bubbles burst first, leading to magmatic foam fragmentation and preferential preservation of small bubbles within ash fragments?

(5) What would be the character of vesicle remnants preserved in eruptions that experienced a complex history of multiple nucleation events and episodic growth?

(6) To what extent and at what scale does heterogeneity within expanding magmatic foam determine the geometry, size, and surface morphology of ash fragments?

While the goal of this paper has been limited to demonstrating the viability of using SSEM to measure prefragmentation BSDs by analyzing ash particles, analysis of even the few ashes we had on hand has led to intriguing implications as well as questions for future research. It is hoped that when applied to specific volcanological problems, this new tool will provide new insights regarding eruptions processes in energetic ash-producing eruptions.

The authors are grateful to Jessica Larsen and Kristi Wallace (Alaska Volcano Observatory, USGS) for providing ash samples. We express special thanks to Professor W. Rose and Owen Mills (Michigan Technological University) for providing access to instruments for SSEM imaging and DEM reconstructions, and to Kevin Grens (Alicona Imaging GmbH) for MeX support. Guilherme Gualda, and two anonymous reviewers greatly improved the manuscript through their comments and thoughtful suggestions, including the addition of an Error Analysis section. This material is based upon work supported by NSF under award numbers EAR-0509859, EAR-0838292 (UNH), EAR-0838314, and EAR-0509856 (Lehigh University).

1We distinguish vesicles from bubbles in that bubbles are a gas phase in a liquid magma or lavas, while vesicles are voids in rock that can be observed and measured in the lab. While we can only observe vesicles, it is the bubbles in which we are interested for physical volcanology.