Quantitative measurement of bubble textures in pumice clasts using a digital stereo microscope with low-angled ring illumination

Quantitative analysis of bubble textures in a large number of volcanic pyroclasts is critical to investigating the eruption dynamics in a volcanic conduit. Here, we used a digital stereo microscope with low-angled ring illumination (DSM-LaRI) to measure bubble textures on unpolished cutting surfaces of pumice clasts. As the DSM-LaRI enhances brightness contrast between the bubbles (pores) and the matrix, we easily obtained the two-dimensional data on the size and shape of bubbles by image analysis. The DSM-LaRI imaging provided the distributions of size and shape of bubbles at least 50 µm across. We applied the DSM-LaRI to analyze more than 1000 pumice clasts from the 232 AD Taupo eruption and measured the mean bubble radius (R¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{R}$$\end{document}) and the mean deformation degree (D¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{D}$$\end{document}) in the individual clasts. The distribution of R¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{R}$$\end{document} and D¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{D}$$\end{document} in each layer showed a distinctive difference between the fallout and the flow deposits. These quantitative data are consistent with a qualitative classification in a previous study. Although the new DSM-LaRI method has the disadvantage of the low spatial resolution, it allows for the analysis of a large number of pumice clasts in a short time, which can address larger scale heterogeneity, by efficiently generating a large representative suite of bubble size and shape data to link bubble textures to conduit processes. This provides vital information for quantitatively modeling eruption dynamics.


Introduction
Bubble textures in eruptive products have been used to infer magmatic processes occurring before, during, and immediately after the eruption. For instance, the bubble number density can be converted to the decompression rate during magma ascent (Toramaru 1995(Toramaru , 2006, and the bubble shape can preserve records of strain rate and shear stress (Rust et al. 2003;Rust and Cashman 2007). The analysis of bubble textures has become an essential research method to reveal the intra-conduit processes during an eruption.
The most common approach to characterize bubble textures is two-dimensional (2D) observation on thin sections and grain mounts (Klug and Cashman 1994;Klug et al. 2002;Shea et al. 2012). Thin-section images taken by a scanning electron microscope (SEM) provide textural information that can be later used to obtain bubble number density and bubble size, in 2D. Then, the data are stereologically converted into 3D textures. This approach has the advantage of high spatial image resolution and can accurately capture the size and shape of microscale bubbles with circular to elliptical shapes. However, the stereological conversion of 2D to 3D data is challenging for elongated polydisperse bubbles, because a cutting section rarely passes through the bubbles' most volumetrically representative parts and the probabilities that a section intersects the individual bubbles cannot be derived analytically (Shea et al. 2010). Recently, the use of 3D imaging by X-ray computed tomography (CT) is increasing in bubble texture analyses on pyroclasts (e.g., Gurioli et al. 2008;Polacci et al. 2010;Giachetti et al. 2011). This approach can accurately visualize bubble textures in 3D. The recent development of high-resolution tomography technique allows us to resolve very thin glass walls of around 1 μm (Degruyter et al. 2010;Dingwell et al. 2016;Fauria et al. 2017). However, subsequent complex coalescence of bubbles, and relaxation of glass can produce irregular bubble networks inappropriate for bubble size distribution analysis (Kennedy et al. 2016).
Along with the microscale bubbles described above, pheno-bubbles (Toramaru 2014), which are large bubbles visible to human eyes, have been studied in recent decades. The macroscopic classification of the pyroclasts based on the visual observation of the pheno-bubbles is helpful to get a general view of bubble textures in an eruption. For example, the stratigraphic change of the relative abundance of tube pumice (consisting of tubelike bubbles) was interpreted to be the result of the temporal evolution of magma ascent velocity and the conduit diameter (Taddeucci and Wohletz 2001). The coexistence of tube pumice and spherical pumice (consisting of almost spherical bubbles) in the same unit has some implications for the velocity distribution across a conduit (Polacci et al. 2001(Polacci et al. , 2003 and for the conduit diameter (Mitchell et al. 2019). The relative abundance of glass shards, i.e., fragments of large bubble walls, is a useful clue to distinguish the fragmentation mechanism between magmatic and phreatomagmatic eruptions (Hiroi and Miyamoto 2016).
Although the macroscopic classification is now commonly used to observe bubble textures in eruption deposits (Taddeucci and Wohletz 2001;Polacci et al. 2003;Houghton et al. 2010), it has the disadvantage of relying on qualitative classification, which may vary from person to person. In addition, the boundaries between pyroclast types are frequently gradational. Digital image analyses are useful for objective classification, but it is technically unrealistic to apply imaging by SEM or X-ray CT to a large number of pyroclasts for obtaining a distribution of bubble textures. Making thin sections for SEM is timeconsuming work, and analyzing bubble textures by X-ray CT is also time consuming especially when resolving thin bubbles walls. The latter technique is not cheap and not overly common, either. For density measurement, it has been recommended to analyze at least 100 pumice clasts to characterize a layer (Shea et al. 2010;Houghton and Carey 2015). To get a quantitative distribution of bubble textures of the whole deposit, an efficient method of digital image analysis is required.
In this study, we develop a method of quickly analyzing visible-size bubbles from unpolished cutting surfaces using a digital stereo microscope with a low-angled ring illumination (hereafter termed DSM-LaRI). It allows us to analyze bubble textures of a statistically sufficient number of pumices (an order of a thousand). We apply the technique to the 232 AD Taupo eruption products, consisting of a variety of bubble textures . Comparing the 2D data by DSM-LaRI with those by SEM and the 3D data by X-ray CT, we evaluate the significance and limitations of the new method.

Materials
The 232 AD Taupo eruption The 232 AD Taupo eruption in New Zealand is a good test case for measuring a wide variety of bubble textures statistically and stratigraphically. The reasons are that its pyroclasts are well preserved and the eruption stratigraphy/chronology is well constrained (Wilson and Walker 1985). This eruption, which ejected phenocryst-poor rhyolite of 35 km 3 DRE, was one of the largest rhyolitic eruptions in the world in the last 5000 years (Wilson and Walker 1985). It formed seven units of eruption products from at least three vents within Lake Taupo (Smith and Houghton 1995). The eruptive products include three phreatomagmatic fall deposits (units 1, 3, 4), Hatepe plinian fall deposit (unit 2), Taupo plinian fall deposit (unit 5a), pyroclastic flow deposit (unit 5b), Taupo ignimbrite (unit 6), and the late-stage dome (unit 7) (Houghton et al. 2003). While magma composition appears homogeneous throughout, the great diversity of eruption styles is attributed to purely physical processes in the eruption, such as magma-water interaction and vent migration (Houghton et al. 2003. Houghton et al. (2003Houghton et al. ( , 2010 reported microscopic bubble textures of pumice clasts from the 232 AD eruption in detail, based on SEM analysis of thin sections. They used textural parameters, such as bubble number densities and bubble size distributions, to constrain the magma ascent history leading to the abrupt shift of eruption styles. Houghton et al. (2010) also conducted the macroscopic classification based on visual observation. They found that the two large dry eruptions (units 5a, 5b, and 6), particularly unit 6, were characterized by an increase in tube pumice and highly vesiculated pumice ( Fig. 4b of Houghton et al. 2010). An increase in tube pumice toward the flow stage has been confirmed in other caldera-forming eruptions (Taddeucci and Wohletz 2001;Polacci et al. 2003). Although its formation mechanism is debated, the stratigraphic difference of abundances of tube pumice may indicate the temporal and spatial variations of shear deformation in conduit flow (Bouvet de Maisonneuve et al. 2009;Dingwell et al. 2016;Ohashi et al. 2018). Therefore, for the purpose of a direct comparison with the data of Houghton et al. (2010), we use pumice clasts from Taupo plinian fall deposits (unit 5a), Taupo ignimbrite (unit 6), and those in the fall deposits of Hatepe plinian (unit 2; a sustained, moderately powerful plinian phase). We do not use samples from the intra-plinian flow (unit 5b).

Field sampling
Samples from units 2, 5a, and 6 were collected from different outcrops (Fig. 1a). The fall deposits of Hatepe plinian (unit 2) were sampled at location A about 25 km away from vent 1. This outcrop shows the whole sequence of the 232 AD Taupo eruption and has been investigated by previous studies (Walker 1980(Walker , 1981. The thickness of unit 2 at location A is about 150 cm (Fig. 1b).
To investigate the transition from a plinian fall to ignimbrite, we conducted detailed sampling of unit 5a at location B, which is about 20 km away from vent 3. This outcrop, located on the main dispersal axis of unit 3, is outside the area covered by intra-plinian flow (unit 5b). In location B, the Taupo plinian fall deposit (unit 5a) can be divided into two sub-units according to sharp changes in color, content of lithic clasts, and grain size of pumice clasts: (1) a fine lower unit (30 cm) containing brown pumice clasts and rich in lithic fragments, and (2) a coarse upper unit (180 cm), characterized by vesiculated white pumice clasts and poor lithic content (Fig. 1b). We collected pyroclasts from the bottom to the top of unit 5a at 20-cm intervals. Only the lowermost section was 10-cm thick. Pumice clasts from the Taupo ignimbrite (unit 6) were collected at location C 48 km away from vent 3.
For measuring the distribution of vesicularity, Houghton and Wilson (1989) and Shea et al. (2010) recommend using at least 100 pumice clasts with the diameter range from 16 to 32 mm . This range is also recommended by Mitchell et al. (2018), who measured bulk vesicularities of a wide range of pumice clasts from the Hatepe plinian (unit 2).
In this study, we collected 100 clasts with a diameter of 16-32 mm in unit 2, but due to the restriction of the baggage weight, we could not bring a sufficient number of clasts for units 5a and 6. When the number of pumices with 16-32 mm was less than 100 , we compensated the measurements with pumice clasts of other sizes. We used 25 clasts from the 8-16 mm size range for unit 5a, and 32-64 mm clasts for unit 6.

Textural analyses using a digital stereo microscope with low-angled ring illumination (DSM-LaRI)
We cut pumice samples along the plane parallel to the visually identified direction of bubble elongation. A standard cutting machine (MC-110, MARUTO) was used. Images of the cutting surfaces were taken with a digital stereo microscope (VHX-1000, Keyence) at  Smith and Houghton (1995). Location A, B, and C are outcrops where we collected samples of Hatepe plinian fall deposit (unit 2), Taupo plinian fall deposit (unit 5a), and Taupo ignimbrite (unit 6), respectively. The isopach of fall deposit of Taupo plinian fall deposit (unit 5a) is from Walker (1980). The values are in centimeters. The area of intra-plinian flow units of Taupo plinian phase (unit 5b) is from Wilson and Walker (1985). b Stratigraphic logs for location A, B, and C. Pumice clasts shown as open symbols; wall rock clasts as closed symbols; fine ash as dashed lines a magnification of 20×. The spatial resolution of an image taken by the DSM was 9.5 µm/pixel (Table 1). The boundaries between a bubble and matrix glass were not visible using the conventional illumination normal to the surface (Figs. 2a and 3a). We illuminated the surface obliquely ( Fig. 2b) using a low-angled ring light attachment (VH-K20, Keyence). We set the illumination angle at the lowest of the adjustable range (10-50° from horizontal according to the personal communication with the manufacturer). The low-angled illumination generated shadows inside the bubbles and clarified their boundaries (Fig. 3b). We adjusted the light intensity to optimize the contrast between bubbles and matrix glasses. The light intensity measured with an illuminometer (HP-881D, HoldPeak) was 4000 ± 500 lux.
Image processing was conducted by Matlab. First, we converted the images taken by the DSM-LaRI to grayscale images and then transformed them into binary images. The ring illumination from the lens's circumference made the center of the image brighter than the edges. To recognize bubbles correctly in the binarization, we used an adaptive threshold algorithm, which computes the threshold of brightness values locally (Bradley and Roth 2007). Second, bubbles smaller than 20 pixels (equivalent to a circular area with a radius of 0.024 mm) were erased as they were indistinguishable from noise in brightness values. Bubbles contacting with the image borders were also removed because they do not retain their original shapes. Third, the trace holes of crystals, which were made in the cutting process, were removed manually. The manual removal of trace holes takes about 20 s per image. Figure 3c shows an example of the obtained binary images. Large bubbles have been segmented correctly, but small bubbles less than a few tens of microns have been uncaptured. Raising the sensitivity of the adaptive threshold of brightness helps to recognize those small bubbles but fails in segmenting bubbles. We empirically determined the optimal value of the sensitivity, so that we can segment large bubbles correctly and recognize small bubbles as much as possible. Nevertheless, the image processing regarded small bubbles immediately adjacent to the large bubbles as a matrix phase. This method should be applied to analyze relatively large bubbles (larger than tens of microns) in a sample.
We used the binary images to measure each bubble area, which we converted to the equivalent bubble radius, R (the radius of a circle whose area is equal to that of the bubble). Each bubble was approximated by an ellipse having the same center of gravity (Fig. 3d). The shape of the fitted ellipse was used for calculating the bubble deformation degree (Taylor 1934).
where a and c are the major and minor axes of the ellipse, respectively. D = 0 corresponds to a circle and D = 1 , an extremely elongated one. The bubble deformation degrees within pyroclasts have been measured for estimating capillary number during magma ascent (Rust et al. 2003;Rust and Cashman 2007).
As a calibration of the DSM-LaRI technique, we measured the size of holes drilled in a test sample. The error in the measurement of the hole radius by DSM-LaRI was less than 0.01 mm. See Additional file 1 for more information.

Thin-section imagery
We performed microscopic analyses of thin sections by SEM to verify the new method by the DSM-LaRI. We selected four representative pumice samples from unit 5a, which exhibited distinct bubble characteristics (Fig. 4). Thin sections were made from the same cutting surfaces captured by the DSM-LaRI. We followed the imaging technique described by Shea et al. (2010). Backscattered electron images of the thin sections were collected with an electron probe microanalyzer (EPMA) with a field emission gun (JEOL, JXA-8530F Plus at the Earthquake Research Institute of the University of Tokyo) at an accelerating voltage of 12 kV. The whole image of images, which were scanned at a range of magnifications from 1× to 500×, were transformed into binary images by setting grayscale thresholds, with spatial resolutions ranging from 0.1 to 30 µm/pixel (Table 1). Bubble walls in pumice are often too thin to be resolved, and are sometimes destroyed while making thin sections. Missing bubble walls were reconstructed by drawing them manually. The final binary images were analyzed by FOAMS (Shea et al. 2010) to obtain bubble sizes and bubble deformation degrees.

3D analyses by X-ray CT
We performed 3D imaging by the X-ray CT with two resolutions. We selected nine pumice clasts from each unit (units 2, 5a, 6) to be measured by low-resolution X-ray CT (METROTOM 800, Carl Zeiss, at Saitama Industrial Technology Center). A sample shaped in a cuboid ( ∼ 1000 mm 3 ) was mounted on the rotational stage in the X-ray CT system and rotated 360 • in steps of 0.45 • ( 800 projections). The configuration used for the tomographic scan was a tube voltage of 130 kV , a tube current of 60 µA, and a voxel size with edge lengths in each direction of 8 µm. The sizes of scanned images are shown in Additional file 1: Table A1. The spatial resolution of the X-ray scanner (METROTOM 800) is not enough to resolve very thin glass walls (a few microns). To reconstruct missing bubble walls, we adopted the marker-controlled watershed algorithm (Meyer and Beucher 1990).
The detailed procedure of the image processing using Fiji (Schindelin et al. 2012) is described in Additional file 1.
To resolve the fine textures, we also performed highresolution X-ray CT for some samples by SkyScan 1272 (Bruker) at Saitama Industrial Technology Center. As the X-ray from this CT scanner is too weak to penetrate a dense pumice clast, we selected a few highly vesiculated pumices from those analyzed by the low-resolution X-ray CT. The samples were cut into small cuboids with a square base of about 1 mm × 1 mm and a height of about 2 mm . They were imaged by a 50 kV-200 µA energy source and rotated 360 • in steps of 0.1 • ( 3600 projections). The spatial resolution is 0.45 µm/voxel. Threedimensional images were reconstructed by the NRecon software. The method of image processing of bubbles is the same as that for the low-resolution X-ray CT. Fig. 3 Digital stereo microscope images. Images a with the conventional light (Fig. 2a) and b with the adjustable low-angled ring illumination (Fig. 2b). c A green/transparent binary image overlying the image of b. The green-colored parts are bubbles recognized by image processing. 'T' indicates trace holes of crystals. d A white/black binary image of the recognized bubbles. Ellipses fitting the bubbles are shown in red. The white bar shows a 5-mm scale in all panels   Table 2). The left column is the binary images of the DSM-LaRI like the one in Fig. 3c. The middle and right are the SEM images with bubbles appearing black. The scale is common within each column SEM because we have polished the cutting surface when making thin sections.

Representative samples
The 2D surface vesicularities measured by the DSM-LaRI are shown in Table 2. They are 41.7-67.1% less than the 3D vesicularities measured by the Archimedes method (Houghton and Wilson 1989).
For the representative samples, Fig. 5 shows the area fractions of bubbles having the individual ranges of the equivalent bubble radius. The inset histograms show the data of all bubbles. The apparent difference between the histograms from DSM-LaRI and SEM is that the former data lack small bubbles ( < 0.05 mm ) seen in the later data. As explained in the methods, the analysis with the DSM has been optimized for the bigger bubbles and the DSM resolution of 9.5 µm is not sufficient to recognize thin bubble walls (microns).
The main histograms show the distribution only for bubbles larger than 0.05 mm to compare the results in the same bubble size. The difference between DSM-LaRI and SEM seems to be subtle in microvesicular pumices (Fig. 5c, d), but non-negligible in coarsely vesiculated pumices (Fig. 5a, b). To quantify the difference, we adopted an edit distance, d edit proposed for comparing two histograms (Cha and Srihari 2002). It is defined as follows: where H DSM-LRI j and H SEM j are the area fractions of the j-th bin in the DSM-LaRI histogram and in the SEM histogram, respectively, and n is the total number of the bins. The edit distance represents the minimum amount of necessary movements to transform one histogram to the other. For a fixed number of bins, the greater the edit distance is, the more dissimilar the two histograms are. Table 2 shows the edit distances between the DSM-LaRI and SEM histograms. The coarsely vesiculated samples show larger edit distances than the microvesicular samples.
In all the four cases, the biggest bubble measured from the DSM-LaRI data is smaller than that from the SEM data, which is caused by the over-segmentation of large bubbles. Although the low-angled illumination is expected to darken bubbles' interiors, it sometimes creates bright parts within bubbles, especially inside large bubbles (Fig. 4a, b). The image analysis identifies the illuminated internal parts as bubble walls and segments a single bubble into smaller ones. We believe the SEM images to exhibit more reliable bubble shapes, representative of the sample sections.
The area fraction of bubbles with the equivalent radii of < 0.1 mm is systematically smaller with the DSM-LaRI compared to the SEM. This difference is particularly noticeable in the coarsely vesiculated samples (Fig. 5a, b). In the microvesicular pumice samples, bubbles between 0.05 and 0.1 mm are surrounded by much smaller bubbles (Fig. 4c, d). Therefore, these 0.05-0.1-mm-diameter bubbles are slightly easier to detect without resolving bubble walls. On the other hand, in the coarsely vesiculated pumice samples, those bubbles often appear to contact each other, and bubble walls are thinner (Fig. 4a,  b). These contacting bubbles are recognized as a single, larger one in the DSM-LaRI; while, the SEM can resolve bubble walls and separate them. For the above reason, the edit distance of the bubble radius histograms from DSM-LaRI and SEM is larger for coarsely vesiculated pumice samples than for the microvesicular samples. Figure 6 shows the area fractions of the bubble deformation degrees. As in Fig. 5, the main graphs show the results for bubbles larger than 0.05 mm. Compared to the bubble size distributions (Fig. 5), the bubble elongation distributions are similar between the DSM-LaRI and  Table 2. It is noted that we can reasonably compare the edit distances calculated for the same set of bins. In the samples with elongated bubbles (Fig. 6b, d), the DSM-LaRI data tend to shift to the smaller values of D than the SEM data, which again is due to the over-segmentation of long bubbles in the DSM-LaRI image processing (Fig. 4b, d).
The X-ray CT makes it possible to measure D for highly deformed bubbles in 3D. Figure 7 shows the relationships between D and bubble size obtained by the three imaging methods for the same pumice sample from Taupo ignimbrite (unit 6). With an increasing bubble radius ( R ), D also increases, which is similar between the DSM-LaRI and low-resolution X-ray CT. We can check the validity of these results using the high-resolution X-ray CT images. The cross-sectional image of the high-resolution X-ray CT shows thin bubble walls of a few microns (Fig. 7c). The relationship between D and R shows an almost constant value of D(∼ 0.65) for bubbles with 0.025 < R < 0.15 mm ( Fig. 7c; green area). On the other hand, the DSM-LaRI and the low-resolution X-ray CT exhibit the smaller D ( 0.4 < D < 0.65 ) in the corresponding range of R (Fig. 7a, b). These results suggest that the DSM-LaRI and the low-resolution X-ray CT may not accurately measure the shape of bubbles with R < 0.15 mm.
The accuracy of the DSM-LaRI also depends on the deviation of the cutting surface from the true direction of bubble elongation. This problem is described in Additional file 1.

Statistical characterization of bubble textures in the outcrop scale
As described earlier, by measuring bubble sizes and shapes through the stratigraphic sequence, we can infer the temporal and spatial variations in magma-filled conduit flow. An underlying assumption is that magmas in individual fragments have ascended along different portions of a conduit, and that they preserve the corresponding deformation histories in their bubble textures. The above assumption is based on a common understanding of laminar flow mechanics. To discuss the dynamics in a conduit, we define textural parameters characterizing a pumice clast.
We employ the average bubble radius and the average bubble deformation degree obtained by the DSM-LaRI for each pumice clast. First, we use sufficiently large bubbles ( R > 0.15 mm ) for the DSM-LaRI resolution (Fig. 7). Second, we derive a weighted average, biased to larger ( > ∼ 0.5 mm) bubbles. Namely, the average bubble radius, R , and the average bubble deformation degree, D , are given by: Fig. 7 Analysis results of tube pumice from Taupo ignimbrite (unit 6). The upper images are the cross-sectional images obtained with a DSM-LaRI, b the low-resolution X-ray CT, and c the high-resolution X-ray CT. The lower graphs show the equivalent bubble radius versus the bubble deformation degree. The blue crosses are all the data, and the yellow dots are the average of D with every radius interval of 2 × 10 −2 mm for a, b and 5 × 10 −3 mm for c. The green-highlighted area shows the range of 0.025 < R < 0.15 mm where R i , D i , and A i are the equivalent bubble radius, the bubble deformation degree, and the bubble area, respectively, of the i-th bubble in one pumice clast measured by the DSM-LaRI. Figure 8 shows the relationships between R and D in units 2, 5a, and 6. The data points corresponding to the representative samples (Figs. 4, 5, 6 and Table 2) are indicated with images. We can see no systematic correlation between R and D in the samples from the 232 AD Taupo eruption. It also demonstrates no apparent clustering. Figure 9 shows the histograms of R and D in units 2, 5a, and 6, as well as the abundances of pumice types reported by Houghton et al. (2010). Although we analyzed the pumice clasts in unit 5a at 20-cm intervals, we found no stratigraphic change in bubble textures ( R and D ). Therefore, we show the sum of all the sub-layers in unit 5a in Fig. 9c. The numbers of the analyzed samples are 100 for units 2 and 6, and 1100 for unit 5a. The measurement results by the DSM-LaRI are consistent with the visual observation by Houghton et al. (2010). Units 5a and 6 include large bubbles, which leads to their broader size distributions than unit 2. The histogram of D in the Taupo ignimbrite (unit 6) is markedly different from the other plinian fall units. The highest peak in unit 6 lies at 0.5 < D < 0.58 , and the distribution is broad from 0.28 to 0.73 . On the other hand, the plinian fall units 2 and 5a are mono-disperse with a peak of around D = 0.35 , which corresponds to less , elongated bubbles. The arithmetic means and standard deviations of R and D are summarized in Table 3.

Integrated methodology for measuring bubble textures
Imaging by the DSM-LaRI provides quantitative bubble textures in the scale of > 0.05 mm . It is not easy to obtain similar amounts of data by the existing 2D approach making thin sections or the 3D approach by X-ray CT because of the long preparation and acquisition time. On the other hand, the developed method with DSM-LaRI cannot resolve microscale bubbles separated by very thin walls. Measuring a wide range of bubble sizes is necessary for accurate number density and size distribution of bubbles. The bubble number density can be converted into the decompression rate in an explosive eruption (Toramaru 2006), and the bubble size distribution is useful to infer the nature of nucleation and/or coalescence events (Klug and Cashman 1994).
For the SEM and X-ray CT analyses, it is necessary to select a small number of representative clasts from each stratigraphic unit. Previous studies analyzed the apparent density and selected several clasts representing low, moderate, and high vesicularity groups (Carey et al. 2009;Shea et al. 2010;Houghton et al. 2010). The selection was based on the assumption that the eruption had lateral vesiculation gradients in a conduit (Shea et al. 2012). When considering the lateral variations of shearing in a conduit (Polacci et al. 2003), the DSM-LaRI analysis will provide the bubble elongation dataset useful for the selection of representative clasts containing spherical, moderately, and highly elongated bubbles. They may have come through different paths with different shear deformation histories in a conduit and recorded the  Table 2. b The number density of crosses in a shown by colors corresponding textures in various scales and shapes. The methodology integrating the DSM-LaRI with the SEM and X-ray analyses could contribute to the investigation of lateral variations of magma ascent.

Textural heterogeneities and implications for magma ascent
We have quantified the heterogeneous bubble textures of the pumice clasts from the 232 AD Taupo eruption and the textural transition from the plinian fall deposits to the ignimbrite (Fig. 9). We will now discuss the development of these heterogeneous textures in terms of bubble size and bubble elongation.
The data obtained using the DSM-LaRI demonstrate that R increases from units 2 to 5a and more to 6 ( Table 3). Post-fragmentation bubble expansion could be a key factor during explosive volcanic eruptions (e.g., Mitchell et al. 2018). Mitchell et al. (2018) performed a detailed investigation of pumice bulk density variations over a size range between 4 and 128 mm . They focussed their study on the Hatepe eruption (unit 2), revealing that bulk density markedly decreased as clast sizes exceeded 32 mm. Although most pumice clasts analyzed in our study are less than 32 mm, pumice clasts might have expanded in the pyroclastic flow (unit 6) during transport and emplacement because pumices were kept at a temperature sufficiently high to allow viscous growth for longer in the pyroclastic flow than in the eruption column. The increase of mean vesicularity from units 5a to 6 is correlated with the increases of R (Table 3). However, bubbles in unit 6 are more elongated than unit 5a,  though post-fragmentation expansion tends to return a deformed bubble to a spherical shape. Pumice clasts, particularly those with deformed bubbles, have fully connected porosities (Spieler et al. 2004;Okumura et al. 2009;Mueller et al. 2011;Colombier et al. 2017). Based on the above considerations, post-fragmentation expansion may have had minimal effect on bubble textures of pumice in the Taupo ignimbrite. A possible error that changes bubble size from units 5a to 6 is the direction of the cutting surface of a pumice. In this study, we cut a pumice clast in a direction parallel to bubble elongation, which apparently increases a measured bubble size. This effect may be remarkable in unit 6 whose bubbles are relatively elongated. The increase in R might have also been caused by bubble coalescence during magma ascent and/ or the initial bubble size at the conduit inlet. Ultimately, the cause of the change in bubble size requires further study.
The data of DSM-LaRI, as well as that of Houghton et al. (2010), show the increase of tube pumice, composed of highly elongated bubbles, from the fallout phases to the flow stage (Fig. 9). Alternatively, an increase in tube pumice abundance within the flow phase has been observed in the Minoan eruption deposits (Taddeucci and Wohletz 2001). Previous studies attributed the formation of tube pumice to the increasing rate of simple shear at the conduit walls (Taddeucci and Wohletz 2001;Polacci et al. 2003;Dingwell et al. 2016). This interpretation explains the coexistence of tube pumice and pumice with spherical bubbles. However, it is still unknown why tube pumice is rare in plinian falls but abundant in ignimbrites that are generally thought associated with wider conduits. The conduit radius for unit 6 was numerically estimated to be 500-600 m (Legros et al. 2000). Supposing a parabolic velocity profile and the same ascent velocity between the plinian fall phase and ignimbrite phase, the average shear rate in the ignimbrite phase should have been smaller than that in the plinian fall phases. Nevertheless, a lot of tube pumices are found in the ignimbrite. The relation between the tube pumice formation and the conduit radius needs to be investigated in the future. The current results from analyzing a large number of samples ( > 100 samples) will enable us to link the dispersion of the bubble elongations to the conduit-scale flow field for each phase of an eruption. In addition to that, the current method based on DSM-LaRI can statistically test the proportions of various textures to determine their relative volumes without sampling bias, leading to a better understanding of conduit processes.

Future work
There is room for improvement in the image acquisition process of the DSM-LaRI method. The current method cannot resolve bubbles glass walls smaller than the DSM resolution of 9.5 µm. Even if we raise the magnification of the lens, it is difficult to resolve them because the cutting surface is rough and the contrast is not sharp between bubbles and glasses. If we could cut a pumice clast more smoothly and illuminate at shallower angles, we could prevent the unexpected illumination of bubble interiors even for microscale bubbles. The image processing of the DSM-LaRI method can also be improved. In the present technique, bright parts inside bubbles, especially large bubbles, are recognized as bubble walls, which results in over-segmentation. This problem could potentially be solved by improving algorithms, with further analysis.
The main advantage of the DSM-LaRI is that it provides us with a large number of digital images from several hundred pyroclasts, which will be useful for developing a statistical analysis method for characterizing disparate sub-units with an eruption deposit. Similar approaches have been developed to characterize volcanic ash particles (e.g., Liu et al. 2015;Avery et al. 2017). For example, Shoji et al. (2018) used a machine learning method to classify ash particles into several clusters, which correspond to eruption styles. The application of the machine learning technique for bubble textures will contribute to the classification of eruption styles and improve our understanding of vesiculation processes in a conduit. An efficient image acquisition technique is essential for this approach.

Conclusions
We measured bubble textures from unpolished cutting surfaces of pumice using a digital stereo microscope with low-angled ring illumination (DSM-LaRI), which produces a good luminance contrast between the flat surfaces and bubbles. The DSM-LaRI significantly reduces the time taken to measure 2D bubble textures compared to the conventional method with SEM. We first compared the results obtained by the DSM-LaRI with those by SEM. We found that the DSM-LaRI method provided the distributions of size and elongation of bubbles consistent with SEM data for bubbles larger than 0.05 mm, though it could not measure microscale bubbles ( < 0.05 mm ). We also performed the X-ray CT analyses and confirmed that the DSM-LaRI showed reasonable bubble deformation degrees for bubbles larger than 0.15 mm . To characterize various bubble textures in an eruption, we define the average bubble radius, R , and bubble deformation degree, D , for each pumice clast. Investigating the histograms of R and D for the individual eruption phases revealed the textural evolution in the 232 AD Taupo eruption. Large bubbles increased from the Hatepe plinian fall to the Taupo plinian fall, and the amount of pumice with