A casting method using contrast-enhanced diethylphthalate for micro-computed tomography of snow

Abstract

Casting snow is necessary to prevent metamorphism and deformation prior to X-ray micro-computed tomography (μCT) imaging. Current methods are insufficient for large-scale field sampling of snow due to safety considerations associated with the casting medium and/or lengthy sample preparation times. Here, a casting method using contrast-enhanced diethylphthalate (DEP) for μCT of snow is presented. The X-ray contrast of DEP is enhanced with barium titanate nanoparticles (BaTiO₃) and iodine (I₂). A partially unsupervised, three-phase segmentation method utilizing traditional Gaussian smoothing followed by a three-step process to address transition voxels is also presented. Synthetic images derived from real snow samples are used to evaluate the segmentation method with various configurations of trapped air bubbles. Real snow samples spanning a range of specific surface areas (SSAs) (8–28 m² kg⁻¹) and densities (135–463 kg m⁻³) are used to assess the performance of the segmentation method on real, cast samples. The method yields SSA, density and correlation length errors of less than 10% for synthetic images with air bubble surface areas less than 333 m⁻¹ per sample volume for eight of the nine snow samples. For eight of the nine cast samples, the method yields errors of less than 10% for all three parameters.

1. Introduction

The microstructure of snow can be used to calculate physical properties such as mechanical strength (Hagenmuller and others, 2014), thermal conductivity (Kaempfer and others, 2005), albedo (Haussener and others, 2012; Ishimoto and others, 2018) and permeability (Courville and others, 2010). These physical properties are relevant to a number of fields including avalanche forecasting (Schweizer and others, 2003), climate modeling (Henderson-Sellers and Wilson, 1983) and remote sensing (Kokhanovsky and others, 2019; Rutter and others, 2019). The principle quantities needed to calculate many of the relevant physical properties are specific surface area (SSA) (surface area per mass of snow (Legagneux and others, 2003)), density and exponential correlation length (Mätzler, 2002). For example, SSA and density can be used to calculate thermal conductivity and air permeability (Calonne and others, 2012), while the correlation length can be used for determining the anisotropy of permittivity (Leinss and others, 2020) or thermal conductivity (Löwe and others, 2013), and for forcing microwave models (Wiesmann and Mätzler, 1999). Currently, several methods such as manual density cutting and optical techniques (Matzl and Schneebeli, 2006; Gallet and others, 2009; Montpetit and others, 2012) can be used to obtain the density and SSA, respectively, but cannot provide the correlation length. To obtain all three quantities with one technique, high-resolution penetrometry (HRP) (Proksch and others, 2015), serial sectioning (Perla and others, 1986) or micro-computed tomography (μCT) (Coléou and others, 2001) must be used. While HRP, serial sectioning and μCT all provide the SSA, density and correlation length, only serial sectioning and μCT provide the structure at the level of detail required by recent physical models (Picard and others, 2018). As serial sectioning is prohibitively slow (ten sections per hour (Perla and others, 1986)), μCT is the most practical technique for obtaining the complete snow microstructure.

One of the main challenges associated with μCT imaging of snow is the relative remoteness of field sites where samples are routinely taken (e.g. Antarctica and Greenland). As it is often impracticable to bring a μCT scanner to these remote locations, the samples must be transported from the field site to a laboratory. During the transportation process (taking days to months), metamorphism of the snow structure occurs and delicate samples may break. Both of these issues can be solved by casting the snow sample, which inhibits metamorphism by displacing the air within the sample and simultaneously provides structural support. Casting therefore allows samples to be transported and stored for long periods of time prior to scanning, which facilitates the inclusion of microstructure data in large, field-campaign data sets (Lemmetyinen and others, 2016; Maslanka and others, 2016; Calonne and others, 2017).

Current casting mediums suffer from drawbacks with respect to safety and sample preparation time. For example, Flin and others (2003) developed a casting technique utilizing 1-chloronaphthalene (chloronaphthalene) as the casting medium. Chloronaphthalene has the advantage of providing sufficient X-ray contrast with ice for direct processing of the μCT images, but is difficult to use under field conditions due to its hazard rating and environmental toxicity. Heggli and others (2009) developed a casting method utilizing diethylphthalate (DEP) as a safer alternative to choloronaphthalene. Unfortunately, DEP does not provide sufficient contrast with ice for direct processing of the μCT images. To work around this, the method developed by Heggli and others (2009) utilizes two scans of every sample: one in the cast state and one after sublimation of the ice phase. Unfortunately, the multiple scans and long sublimation step (7 d) result in long processing times; the processing of a single, 1 m snow profile takes approximately 1 month. As the most common application for snow casting is remote field work in which the snowpack is sampled at high spatial resolutions both vertically (depth profile) and horizontally (location dependence), the limitations of the methods described above are prohibitive. Therefore, a new casting medium is needed.

An ideal casting medium would have the following properties: a freezing point between 0°C and approximately − 10°C, low water solubility, the ability to penetrate the snow pore space, high X-ray contrast with ice, safe for use in the field, low cost and sufficient thermodynamic stability over time. As DEP already meets the requirements listed above, except contrast, it provides a logical starting point for developing a new casting medium. The challenge is improving the contrast of DEP without compromising any of these physical, chemical or safety requirements.

In μCT, contrast is determined by differences in X-ray attenuation. For a given X-ray energy, attenuation is dependent on the atomic number and density of a compound, which is why iodine containing compounds are often used as contrast agents in clinical CT of soft tissue (Hrvoje and Greenstaff, 2014), for example. Unfortunately, the contrast agents used in medical CT imaging are generally incompatible with organic solvents such as DEP, because medical contrast agents are designed to be polar for effective dissolution or dispersion in blood (Hrvoje and Greenstaff, 2014). Related to solvent compatibility is the issue of thermodynamic stability. The contrast agent should be stable in DEP for a minimum of approximately 1 week to be practical for field use. Certain compounds, such as colloids, are dispersed and are not necessarily stable over long periods of time, depending on their formulation. Solutions, on the other hand, are stable, but affect colligative properties such as freezing point depression. Finally, the contrast agent must be sufficiently inexpensive as a single 1 m snow profile can require multiple liters of casting medium, depending on the snow density and sample size.

After casting, a segmentation technique is necessary to process the μCT images. Similar to the casting medium, the segmentation technique must also meet certain requirements; it should be accurate, fast and unsupervised. For the intended use (snow casting), unsupervised means being able to reproducibly process multiple μCT images of various snow types without needing to adjust or tune parameters. For cast snow samples, this requires a technique capable of handling the three-phase system that results from air bubbles becoming trapped within the sample during casting. Three-phase segmentation techniques are typically more complex and require longer computation times than their binary counterparts. The primary challenge is addressing the presence of mixed voxels (3D pixels). Mixed voxels are created by the partial volume effect, which is caused by the finite resolution of the μCT detector, an X-ray beam that is not a perfect point, and the mismatch of the detector pixels with the phase boundaries within the sample. Therefore, voxels on phase boundaries end up with intensity values between the intensities of the pure phases making up the boundary. For cast snow, this results in a population of ‘transition voxels’ between air bubbles and the casting medium, which appear to be ice, but are actually just mixed voxels. These transition voxels cannot be properly segmented with a simple threshold-based technique, such as that used by Heggli and others (2009).

To tackle the transition voxel problem, a number of segmentation techniques have been developed. Hagenmuller and others (2013) developed an energy-based segmentation, which assigns the transition voxels based on minimization of interfacial area. This segmentation technique is accurate and improves upon the sequential segmentation technique by Flin and others (2003). The method used by Heggli and others (2009) circumvents the transition voxel issue altogether by only performing binary, threshold-based segmentations. As the X-ray attenuation of ice and DEP is similar, the images from both scans in this method are effectively two-phase; the cast image consists of an air phase and a combined DEP/ice phase, while the post-sublimation image consists of only air and DEP. However, as discussed above, this method is not practical due to the long sublimation step. Additional three-phase segmentation techniques are found in other fields, such as medicine (Bromiley and Thacker, 2008), geomaterials (Bhattad and others, 2010; Hashemi and others, 2014; Taniguchi and others, 2014; Shah and others, 2016), food processing (Pinzer and others, 2012; Masselot and others, 2021), composites (Bale and others, 2012; Moroni and others, 2016) and sea ice (Lieb-Lappen and others, 2017). These techniques also have their own drawbacks. For example, Pinzer and others (2012) used an iterative approach to correctly segment transition voxels based on the expected air volume fraction. This requires prior knowledge of the volume fraction of air, which is not possible for cast snow samples. The technique proposed by Hashemi and others (2014) uses simultaneous region growing, which requires several manually set parameters. This introduces operator error and also precludes automated processing of multiple images. Therefore, the key challenge to any new, three-phase segmentation technique is the handling of transition voxels.

In this work, we demonstrate a casting method using contrast-enhanced diethylphthalate for μCT of snow. Barium titanate (BaTiO₃) and iodine (I ₂) are used to improve the contrast of DEP without compromising the necessary physical, chemical and safety requirements discussed above. Furthermore, they are sufficiently cheap and stable in DEP for large-scale field use. The accompanying three-phase segmentation method used for image analysis utilizes a traditional Gaussian filter followed by a three-step procedure to deal with transition voxels; the resulting process is accurate, fast and partially unsupervised. In the following sections, we detail our experiments for assessing the proposed casting and segmentation method. First, synthetic images derived from real snow samples are used to validate the segmentation method as well as to assess the performance with various configurations of trapped air bubbles. Second, the casting and segmentation method is used to compare the SSA, density and correlation length of real snow samples prior to and after casting.

2. Materials and methods

2.1. Experimental

2.1.1. Contrast agents

BaTiO₃ nanoparticles (99.9 wt%, 100 nm diameter) were purchased from GetNanoMaterials (St-Cannat, France) in powder form. I ₂ (≥99.5 wt%, round particles) and DEP (≥99 wt%) were purchased from Sigma Aldrich (Buchs, Switzerland). For the BaTiO₃/DEP dispersion, 25.19 g of BaTiO₃ was added to 560.52 g of DEP at room temperature, resulting in a 4.3 wt% dispersion, and mixed with a Kinematica (Luzern, Switzerland) Polytron PT 10/35 mixer for 20 min at approximately half of maximum speed (exact speed is unknown due to age of instrument). For the I ₂/DEP solution, 16.81 g of I ₂ was added to 562.77 g of DEP at room temperature, resulting in a 2.9 wt% solution, and mixed with a Kinematica Polytron PT 10/35 mixer for a few minutes until no solid iodine was visible in the bottom of the flask. The mixing speed for the iodine solution was lower than for the barium titanate disperion because the mixing was only needed to accelerate iodine dissolution, as opposed to dispersing the barium titanate nanoparticles. Both mixtures were transferred to a − 5°C chamber and allowed to equilibrate prior to use. For the purposes of this paper, the term DEPc will be used to generically describe DEP with contrast agent; when necessary, the contrast agent will be specified.

2.1.2. Snow samples

Samples in this study were collected from six snow sources. S1, S2, S5 and S6 were collected from natural snow samples, including a buried surface hoar layer (S5). S3 and S4 were produced in the laboratory using the device described by Schleef and others (2014a). Of the six types, three (S1, S2 and S6) were cast with both BaTiO₃ and I ₂. The iodine castings of the remaining samples (S3, S4 and S5) could not be used due to issues encountered during the casting process. For the S5_I sample, the weak layer broke prior to casting and for S4_I and S3_I, I ₂/DEP segregation rendered the samples unusable. This resulted in a total of nine samples evaluated here. Their 3D structures are provided in Figure 1 and a summary of their properties is provided in Table 1. For completeness, a subjective snow classification according to Fierz and others (2009) is also provided in Table 1.

Fig. 1. 3D images of the snow samples used: (a) S1_I, (b) S1_B, (c) S2_I, (d) S2_B, (e) S3_B, (f) S4_B, (g) S5_B, (h) S6_I and (i) S6_B. All volumes are 300 × 300 × 300 voxels (300 voxels = 5.4 mm), except for g), which is 500 × 500 × 500 voxels (500 voxels = 9.0 mm) in order to capture the large, vertical structures.

Table 1. Summary of snow samples used and their properties

Errors are absolute.

Sample volumes were 600 × 600 × 600 voxels for all samples except S5_B and S6_I. For the synthetic experiment, both S5_B and S6_I had volumes of 500 × 500 × 500 voxels, while they had volumes of 551 × 551 × 501 voxels and 500 × 500 × 600 voxels, for the cast experiment, respectively. The differences in sample volumes for S5_B and S6_I were due to issues encountered during casting.

2.1.3. Casting

Snow samples were cast following the method by Heggli and others (2009). Snow that had been stored at − 23°C was placed into a polypropylene cylinder (diameter of 30 mm) with a piece of metal on top to prevent the sample from floating (density of DEP is 1200 kg m;⁻³, density of ice is 917 kg m;⁻³) or rotating. The sample holder was then placed into a μCT holder (diameter of 37 mm) and secured in place. Next, the sample was cast with DEPc that had been precooled to − 5°C (room temperature also − 5°C). The casting in this method was performed in batches, as opposed to continuously as described by Heggli and others (2009), but a similar result is expected. In both methods, the aim is to provide sufficient time for the casting medium to penetrate the pore space via capillary action. This is done in an attempt to mitigate the trapping of air bubbles within the structure (Perla and others, 1986). Therefore, the height of the liquid outside of the snow structure was maintained below the height of the liquid within the pores of the snow structure. Once sufficient DEPc had been added to fill the snow structure, the sample was sealed with tape and moved to − 60°C where it remained until frozen.

For the comparison of the uncast (ice/air) sample with the cast (ice/air/DEPc) sample, a baseline μCT scan was first taken of the uncast sample. The entire sample holder was then carefully removed from the μCT scanner and the snow sample was cast directly within the holder as described above. Care was taken to avoid movements that could have affected the snow structure. After the DEPc completely froze, the sample holder was returned to the scanner. Due to the design of the sample holder, the sample retained its orientation for both scans. As snow constantly metamorphoses in air, the sample holder was covered with tape during scanning of the uncast samples and the samples were cast as quickly as possible after the initial scan (within 6–12 h). In between scanning and casting, samples remained covered with tape and were stored at − 15°C. In addition to preventing metamorphism, sealing the samples also prevented iodine from sublimating and depositing on surfaces within the scanner. Representative images of uncast and cast samples are shown in Figure 2.

Fig. 2. Representative images of the S1 samples prior to and after casting: (a) S1_B prior to casting, (b) S1_B after casting with BaTiO₃, (c) S1_I prior to casting and (d) S1_I after casting with I ₂. Note that the respective pre- and post-casting structures are slightly different even though the same region of all samples is shown. This is due to small shifts in the samples between scans. The axes’ units are voxels.

2.1.4. X-ray tomography

μCT scans were performed on a μCT 40 by SCANCO Medical AG (Brüttisellen, Switzerland). X-ray tube settings of 70 kVp and 8 W were used for all samples. For all samples, a total exposure time of 600 ms was used with 1000 projections per slice (180° of rotation). The resolution was 18 μm.

2.2. Image processing

2.2.1. General

Image processing scripts were written in Python 2.7 using standard Python libraries as well as the Simple ITK library from Kitware, Inc. (Clifton Park, NY, USA) (Lowekamp and others, 2013; Yaniv and others, 2018). Initial volume selections were made with SCANCO Medical AG's in-house software to reduce the file size. All further steps were performed in Python. SSA, density and (directional) correlation length were calculated from the correlation function as described by Debye and others (Debye and others, 1957; Mätzler, 2002; Proksch and others, 2015; Krol and Löwe, 2016). The air bubble surface area per sample volume (ASAV), calculated as the surface area of the air bubbles divided by the total sample volume, is used to characterize the air phase in the cast and synthetic images. In this paper, the terms SSA, density and correlation length always refer to the ice phase, while the ASAV is used exclusively with respect to the air phase. The image processing scripts were run on an 8-core, 2-processor, 2.9 GHz machine with 128 GB of RAM.

2.2.2. Segmentation

In the following description, all steps were performed on 3D volumes. For uncast samples (i.e. ice/air), the grayscale image was blurred with a Gaussian filter (width σ = 1.2 and support= 2) and then segmented using the threshold-based technique presented by Hagenmuller and others (2013), which utilizes a mixture model with three Gaussian curves to determine the binary segmentation threshold.

For cast samples, the following proposed method (outlined in Fig. 3) was used. First, the grayscale image was segmented as described above for the uncast samples, except the Gaussian filter had a σ value of 2 instead of 1.2. This was necessary because the cast images were noisier than the uncast images. Next, a three-step procedure consisting of distance mapping, gradient thresholding and component labeling was used to address transition voxels. Note that, as shown in Figure 3, the distance mapping and gradient thresholding was performed on the blurred grayscale image, but the reassignment of voxels was performed on versions of the segmented image. For distance mapping, the blurred grayscale image was segmented into two separate images using the approach by Hagenmuller and others (2013). The binary segmentation threshold for one image was μ _ice–2σ _ice and μ _ice + 2σ _ice for the other image, where μ _ice and σ _ice are the mean value and standard deviation of the Gaussian curve fit to the ice phase, respectively. Both segmented images were then passed through the Simple ITK filter SignedDanielssonDistanceMapImageFilter() and voxels with absolute Danielson distances less than or equal to 1.4 in both images were reassigned (performed on the ‘segmented’ image as referenced in Fig. 3) to the DEPc phase; whether the voxels were assigned to the air or DEPc phases was irrelevant since both phases have the same value in the final, binary image. The goal with this condition was to identify voxels that were in between an air and DEPc phase, which is the case for the transition voxels around air bubbles. A Danielson distance of 1.4 was found to provide maximum removal of air bubble ring structures without removing other voxels that were already correctly assigned (i.e. those of the ice structure). Although this approach does not take into account direction or distance to each specific phase (such as by Hagenmuller and others (2013)), it was found to be sufficiently accurate and fast. Parallel to the distance mapping procedure, the blurred grayscale image was passed to the Simple ITK filter GradientMagnitudeFilter(), which calculates the maximum gradient of each voxel. Voxels with gradients greater than the 99.5 percentile as determined by the cumulative distribution function (CDF) were reassigned (performed on the ‘segmented after distance mapping adjustment’ image as referenced in Fig. 3) to the DEPc phase. This was done because the largest gradients occur at the boundary between air and DEPc. The 99.5 percentile was experimentally determined to provide maximum removal of air bubble rings without removing large numbers of ice phase voxels. Finally, the segmented image, now with both distance mapping and gradient thresholding adjustments, was passed to the Simple ITK filter ConnectedComponentImageFilter(), which was used to identify small ice-like volumes that were not connected to the main ice structure. The threshold for removing these objects was one half of the volume of a sphere with diameter equal to the optical diameter. The optical diameter was calculated as

(1)$$D_{\rm optical} = \left({6\over \rho\, SSA} \right)$$

where D _optical is the optical diameter, ρ is the density of ice and SSA is the SSA of the sample (Leppänen and others, 2015). The one half was chosen as a compromise between removing real ice structures and capturing larger regions of artifacts. The component labeling was only performed on voxels that were not touching the volume border, because it is not clear if these voxels are partial ice grains or true artifacts. After component labeling, the segmentation procedure was complete. For the cast samples, computation time was found to be approximately 15 min for volumes of 600 × 600 × 600 voxels.

Fig. 3. A flow chart of the segmentation process. The images in full lines are all done ‘in-place’ while the distance mapping and gradient thresholding processes are performed in parallel (dashed lines). Where the dashed arrows rejoin the full lines is where the adjustments from these processes take place. The double line around the bottom box indicates that this is the final version of the segmented image.

2.2.3. Method verification with simulated cast images

In order to determine the accuracy of the proposed image processing method, a test was performed comparing a truth image and a semi-synthetic image (referred to as the ‘synthetic image’ from now on). The goal was to validate the segmentation algorithm without the influence of errors caused by the casting process. The synthetic images were also used to simulate different outcomes of the casting process with respect to the distribution of air trapped within the ice structure. A visual overview of the synthetic image generation process is provided in Figure 4 and the general steps are outlined in the following list.

Fig. 4. An overview of the steps used to create the synthetic images: (a) uncast image of S1_I, (b) segmented version of the uncast grayscale image (this is the truth image), (c) air phase generated by the boolean model after correction for the existing ice structure, (d) ideal three-phase image, (e) blurred three-phase image and (f) final synthetic image after addition of scaled μCT noise. The axes’ units are voxels.

1. (1) Begin with an uncast grayscale image (Fig. 4a)

2. (2) Segment the uncast grayscale image (Fig. 4b)

3. (3) Generate air bubbles based on volume fraction and radius (Fig. 4c)

4. (4) Assign the remaining void space (i.e. not already ice or air bubbles) as DEPc and set the grayscale intensity value for each of the three phases (air, DEPc, ice) (Fig. 4d)

5. (5) Apply a Gaussian filter to blur the image (Fig. 4e)

6. (6) Add noise from a scan of an empty μCT sample holder (Fig. 4f)

Truth images were created by segmenting uncast grayscale images. Air bubbles were artificially added to the truth images by generating a microstructure from a boolean model of monodisperse spheres with a fixed seed value (Torquato, 2002). Seventeen synthetic images were generated using nominal air volume fractions of 0, 0.001, 0.01, 0.02 and 0.1, and air bubble radii of 90, 180, 360 and 720 μm. Where the air and ice structures overlapped, the ice structure was retained, which is why the air volume fractions are nominal values. After addition of the air bubbles, the remaining void space was assigned to the DEPc phase, resulting in an image consisting of three phases representing air, ice and DEPc. Each phase was then assigned a grayscale value based on the known mean values from cast samples. The chosen values depended on a number of factors including the contrast agent being mimicked (BaTiO₃ or I ₂), its concentration in DEP and the μCT settings. The result was a ‘ideal’ three-phase image with perfectly defined phases. The process used by Hagenmuller and others (2013) was then used to convert this ideal image into a synthetic grayscale image. First, the ideal three-phase image was blurred using a Gaussian filter (width σ = 1.2 and support= 2), after which noise from a scan of an empty μCT sample holder was added. Noise from a real scan was chosen over uncorrelated Gaussian (white) noise, because it more accurately represents the images produced by the μCT scanner; the μCT noise exhibits short range spatial correlations. The noise was scaled such that the resulting histogram closely matched a representative cast sample as shown in Figure 5. The final synthetic images were then passed through the segmentation steps described above and the resulting SSA, density and correlation length were compared to the truth values. Synthetic images were generated for each of the nine samples listed in Table 1. An overview showing nine of the 17 synthetic configurations generated for the S1_I sample is provided in Figure 6.

Fig. 5. A comparison of (a) the cast image of S1_I, (b) the synthetic image of S1_I with a nominal air volume fraction of 0.02 and radius of 40 voxels and (c) the resulting histograms. The axes’ units in (a) and (b) are voxels.

Fig. 6. An overview of the synthetic air bubble configurations generated. All images are based on the S1_I sample and the ideal three-phase image is used for clarity. Structures in (a), (b) and (c) were generated with an air bubble radius of 5 voxels and nominal air volume fractions of 0.001, 0.02 and 0.1, respectively. Structures in (d), (e) and (f) were generated with an air bubble radius of 40 voxels and nominal air volume fractions of 0.001, 0.02 and 0.1, respectively. The axes’ units are voxels.

3. Results

A step-wise example of the segmentation process is provided in Figure 7 for a synthetic image based on the uncast S1_I sample with an air bubbles radius of 40 voxels and a nominal air volume fraction of 0.02. Figures 7a, b show the uncast grayscale and segmented images, respectively. The progressive removal of the transition voxels by the distance mapping, gradient thresholding and component labeling can be seen in Figures 7c–e, respectively. Figures 7f1–f4 show zoomed in views of one air bubble corresponding to Figures 7b–e, respectively.

Fig. 7. An example of the segmentation steps for a synthetic image of S1_I with an air bubble radius of 40 voxels and a nominal air volume fraction of 0.02. The steps are shown as (a) starting grayscale synthetic image, (b) segmented image after the Gaussian filter, (c) segmented image after distance mapping, (d) segmented image after gradient thresholding, (e) final segmented image (i.e. after component labeling) and (f1–f4) zoomed in views of an air bubble from (b), (c), (d) and (e), respectively. The axes’ units are voxels.

3.1. Synthetic images

Comparisons of the SSA, density and exponential correlation length between truth values and those calculated from the synthetic images are provided in Figures 8a–c, respectively. A 1:1 line is provided in each figure as a visual guide and a zoomed in view is provided in Figure 8a to better differentiate the samples with similar SSAs. The colorbars indicate the nominal ASAV, which is nominal due to the use of the nominal air volume fraction as described above; ‘nominal’ will be omitted for the sake of brevity in the remainder of this paper, except when necessary to distinguish between cast and synthetic values. The 0 m⁻¹ ASAV points are not shown on the plots; this was necessary to use a logarithmic scale in the colorbar to differentiate between points at the lower end of the ASAV range; the values for these points are provided separately in Table 1. Average and median errors, as well as standard deviations for the data in Figure 8 (including 0 m⁻¹ ASAV points) are also provided in Table 1. An example correlation function, from which SSA, density and correlation lengths are derived, is provided in Figure 9.

Fig. 8. Results of the synthetic comparison between the truth and synthetic images for (a) SSA, (b) density and (c) correlation length. The color bar in the lower right applies to all three plots.

Fig. 9. The correlation function for synthetic and truth images of the S1_I sample with air bubble volume fraction of 0.02 and radius of 40 voxels.

As shown by the standard deviations in Table 1, and confirmed visually in Figure 8, the error distribution is largest for SSA, followed by density, and finally, the correlation length. The largest errors for all three parameters occur at the largest ASAV values.

The SSA generally decreases between ASAV values of 0 m⁻¹ and approximately 167 m⁻¹, after which it increases. There are a number of points that do not follow the trend rigidly, but these points typically occur near the transition point and are small (SSA only changes a few percent in the unexpected direction). The method provides SSA errors less than 10% for ASAV values less than or equal to 333 m⁻¹ for all samples.

In contrast to SSA, we do not see any decrease in density with increasing ASAV. Instead, the density generally increases with increasing ASAV. Similar to the SSA trend, several samples have small, out of place increases or decreases, where the density changes by a few percent. The S1_I, S1_B and S4_B samples all show larger outliers at 417 m⁻¹; the aforementioned trend holds prior to and after the 417 m⁻¹ points for these samples. Ignoring S4_B (because all points have errors greater than 10%), the method provides density errors less than 10% for ASAV values less than or equal to 833 m⁻¹.

For correlation length, all samples except S1_I, S1_B and S5_B show an increasing correlation length with increasing ASAV at low ASAV values, followed by a direction change to decreasing correlation length with increasing ASAV at large ASAV values. The point at which this direction change occurs is not as consistent across all samples as was observed for SSA. For S5_B, the correlation length generally increases with increasing ASAV, while for S1_I and S1_B, the correlation length generally decreases. Ignoring S4_B (because it is the only sample with errors greater than 10%), the method provides density errors less than 10% for all ASAV values tested here (up to 3333 m⁻¹). Including S4_B, the 10% threshold is valid for ASAV values less than or equal to 42 m⁻¹.

3.2. Cast images

Figure 10 provides a comparison of the SSA, density and correlation length of real snow samples prior to and after casting. The comparison is analogous to Figure 8. Errors and calculated ASAV values are provided in Table 1. As each sample could only be cast once, there is only one point per sample.

Fig. 10. Results of the comparison between the uncast and cast images for (a) SSA, (b) density and (c) correlation length.

S5_B shows the largest SSA error and is the only sample with an SSA error above 10%. S4_B shows the largest error for both density and correlation length and, for both parameters, is the only sample with an error above 10%. All other samples have errors less than 10% for all three parameters.

For SSA, all samples except S4_B are above the 1:1 line. For density, all samples with uncast densities above approximately 254 kg m⁻³ are below the 1:1 line, while samples with lower densities are above. For correlation length, all samples except S5_B are above the 1:1 line. Points above the 1:1 line indicate that the cast value was overestimated with respect to the uncast value, while points below the 1:1 line indicate an underestimation.

Table 1 provides a comparison between the cast samples and the synthetic image with the closest nominal ASAV. For some nominal ASAV values, two synthetic configurations led to the same nominal ASAV. In these situations, the two values were averaged. This averaging had little effect on the results, because, for these points, the values (SSA, density, correlation length) differed by a maximum of 0.6% (average of 0.1%). The largest calculated ASAV value for the cast samples was 384 m⁻¹ for S3_B, while the other samples had values less than or equal to 172 m⁻¹. Table 1 also shows the difference between the cast and synthetic values for SSA, density and correlation length, in percent (with respect to the synthetic value). This was calculated as

(2)$$ Error ( \% ) = \left({{Cast-Synthetic}\over Synthetic} \right)\times \! 100 $$

where Cast and Synthetic refer to the cast and synthetic values of each of the three parameters, respectively. Therefore, a negative error indicates that the synthetic value was greater than the cast value. This comparison of the cast samples and the synthetic point at the closest nominal ASAV is used to determine if the segmentation method is consistent in how it under- or overestimates the three parameters. For SSA, all samples except S1_I and S1_B are consistent. For density and correlation length, all samples except S5_B are consistent.

3.3. Error summary

Figure 11 shows the absolute error of SSA, density and correlation length for every synthetic and cast sample evaluated as a function of ASAV (nominal for synthetic values and calculated for cast samples). The dashed line shows the 10% error threshold. This plot provides a visual representation of many of the trends and results discussed above. For example, it shows that the correlation length is relatively insensitive to ASAV compared with SSA and density, which are more sensitive. It also shows that, for both SSA and density, the error generally increases with increasing ASAV. Furthermore, we can see details such as the large density errors associated S4_B as the open, red circles above the 10% error line at low ASAV.

Fig. 11. A summary of the absolute error of SSA, density and correlation length for all synthetic and cast images as a function of ASAV. The ASAV is the nominal value for the synthetic images and calculated for the cast images. The error is with respect to the truth and uncast values for the synthetic and cast images, respectively.

4. Discussion

4.1. Casting

The concentrations of BaTiO₃ (4 wt%) and I ₂ (3 wt%) were selected to provide sufficient contrast while minimizing the amount of materials used. This is advantageous from a cost perspective as well as to minimize changes in the physical properties of the DEP after addition of the contrast agents. For example, because solubility decreases with decreasing temperature, an iodine concentration below the maximum solubility limit (experimentally determined to lie between approximately 5 and 14 wt% at room temperature) was chosen in order to avoid precipitation of iodine when the solution is cooled. For BaTiO₃, on the other hand, freezing point depression is not expected, because the mixture is a dispersion (You and others, 2016). Instead, the rheological properties of dispersions are known to be affected by high loadings, such as those used in the semiconductor industry (Mikeska and Cannon, 1988; Moulson and Herbert, 2003). This could, for example, impact the casting medium's ability to penetrate pore spaces. Because both contrast agents are shown here to be viable for snow casting, the contrast agent used for a specific experiment or field campaign may be selected based on the available equipment, resources and contrast required.

As discussed in the Materials and Methods section, the DEPc was precooled to − 5°C prior to casting. Although this is below the literature value for the freezing point of pure DEP (− 3.25°C (Chang and others, 1967)), the DEP and DEPc used in these experiments were not found to freeze at − 5°C, over several months for DEPc and on the order of years for DEP. As DEPc with both BaTiO₃ and I ₂ remained liquid at − 5°C, this is unrelated to the expected freezing point depression caused by the iodine. DEPc samples were found to freeze at temperatures up to − 15°C, and once frozen, were found to remain frozen up to − 10°C. Samples frozen at − 15 °C were found to take longer (additional days to weeks) compared to those at − 60°C. The freezing process at both − 15°C and − 60°C could be initiated by adding a small piece of frozen DEP or DEPc to the subcooled DEPc (i.e. − 15°C or colder). After addition, the freezing process was rapid; crystals immediately began forming near the frozen piece of DEP or DEPc. Samples frozen in this way exhibited segregation of the contrast agents from the DEP. This segregation is consistent with the behavior observed for the freezing of supercooled liquid solutions and dispersions (You and others, 2016), but is undesirable for this application, because it removes contrast from the bulk DEPc phase. The segregated structures were observed to occur most readily in areas of the samples with little or no snow; the segregation appeared to be inhibited by the presence of snow. This bodes well for field applications where the samples are often cast in larger cross-sections than the 30 mm diameters used in this study, because only the middle of the sample is ultimately scanned and it is hypothesized that any segregation would likely stop before penetrating the core of the sample. This could still be an issue for low density snow (new snow), because the density may be insufficient to prevent segregation. Fortunately, the use of DEPc does not preclude processing the sample by the method proposed by Heggli and others (2009). It should be noted that some samples also exhibited segregation even without addition of frozen DEP or DEPc. For example, this occurred for the S3_I and S4_I samples, which is why they were not included in this work. It is assumed that some type of contamination induced the freezing process in a similar way, but this could not be fully investigated.

As mentioned above, samples took longer to freeze at higher temperatures and were therefore frozen at − 60°C. Although − 60°C freezers are not always available during field work, dry ice (sublimation temperature of − 78.5°C (Perry and others, 2000)) and liquid nitrogen (melting point of − 210°C (Perry and others, 2000)) are relatively widely available and already used during field campaigns. We therefore believe that the freezing conditions in this study are representative of those achievable in the field with dry ice or liquid nitrogen.

The challenges associated with freezing were not investigated in detail, but the literature suggests that DEP may be relatively resistant to nucleation. Chang and others (1967) showed that homogeneous nucleation in pure DEP only occurs below − 53°C and we found that some samples of DEPc with BaTiO₃ (no snow) could remain liquid at − 15°C for multiple months. The DEP purity used by Chang and others (1967) (measured as 99.88 wt% assuming the impurities have the same molecular weight as DEP) is similar to the DEP used in this study (≥99 wt%) (Chang and others, 1967). However, due to expected contamination of the snow samples by dust and the presence of snow itself, heterogeneous, not homogeneous, nucleation should be expected. This is confirmed by samples that froze at − 15°C, but the nucleation process was slow (days to weeks) at this temperature. Once the freezing process began (i.e. a portion of the sample was visibly frozen), the remainder of the sample generally froze within a few hours. Further investigation is needed to better understand the freezing behavior of DEP, but this was outside of the scope of this work.

During the casting process, precautions were taken to avoid changing the snow structure or its orientation within the sample holder. Unfortunately, even with the precautions taken, small shifts of the samples during casting were inevitable. This was likely partially due to handling of the sample and partially due to the buoyancy of the snow in DEP, which can cause weakly bound portions of the structure to float or shift. Additionally, it is assumed that some degree of mechanical deformation or destruction of the snow structure is also caused by the surface tension of the casting liquid. Investigation of this effect, however, is outside of the scope of this work.

As mentioned above, the small casting volume used in this work differs from that typically used in the field. In addition to potentially reducing the impact of contrast agent segregation, the larger field volumes may also alleviate issues related to sensitive snow layers, such as the weak layer in S5. As the cast cross-section becomes smaller, each broken bond between grains represents a larger fraction of the total layer strength and is therefore more critical for preventing collapse. Larger volumes reduce this effect and make the layer less likely to collapse.

4.2. μCT

The μCT settings play an important role in the quality of the images and the subsequent performance of the segmentation method. Increasing the X-ray tube voltage not only decreases the noise, but also decreases the contrast. Increasing the power has the same effect. Longer integration times and averaging may be used to decrease the noise, but this increases the scan time. A balance was chosen to minimize scan time while reducing the noise to a level for which the method performed well. It should be noted that the X-ray attenuation of a compound is highly dependent on the X-ray energy distribution due to the presence of absorption edges (large changes in absorption at energies equal to an electron shell binding energy). As we did not observe any unusually large changes in contrast by changing the X-ray tube settings, we assume that the energy distributions were similar enough near the absorption edges of the contrast agents to inhibit dramatic changes. Optimizing the X-ray energies for particular contrast agents could further improve contrast.

With regard to noise, all barium titanate containing samples were more noisy than their iodine counterparts. For every snow type (except S3, S4 and S5 due to issues described above), one casting was obtained with BaTiO₃ and one with I ₂. The standard deviation of the noise in each barium titanate sample was approximately 1.7 times as large as the standard deviation of the noise in the corresponding iodine samples. This likely has to do with the interaction of the contrast agents and the X-rays. Generally speaking, higher attenuation materials increase the noise of the image, because the image is based on fewer photons reaching the detector (for a given flux and exposure time). This is why the cast samples are noisier than the uncast samples. However, for the contrast agents used here, we would have expected the I ₂ containing images to be noisier, because the I ₂ samples showed greater attenuation than the BaTiO₃ samples. Since the opposite was observed, something else must be at play. For example, it is known that the presence of high attenuation substances can lead to additional image artifacts through beam hardening and scattering (Joseph and Spital, 1982). Additionally, the nanoparticles may have been in a partially aggregated state due to imperfect dispersion with the mixer or due to aggregation post mixing, but this was not investigated.

It should be additionally noted that the noise used for the generation of the synthetic images contained ring artifacts from the μCT reconstruction. These occur in the center of the detector and appear as concentric circles. In the synthetic images, they appear in the bottom right corner (e.g. Fig. 4). This adds not only additional error to the synthetic images, but also makes them more realistic than if uncorrelated Gaussian noise had been used. Ring artifacts may be corrected in post-processing, but this was not performed here.

4.3. Segmentation method

The segmentation method is described as partially unsupervised because multiple values (blurring parameters, gradient threshold and Danielson distance) are indeed manually set and may not be generally applicable for all systems. However, for a fixed setup (i.e. scanner setup, sample volume, casting agent concentration), these values do not need to be changed. As a typical laboratory will use the same scanner, once these parameters are set, multiple cast samples may be analyzed without further adjustments.

In general, the performance of the segmentation method is affected by a number of factors. First, the σ value used for the Gaussian blurring of the synthetic and cast samples is larger than the σ value (1.2) used for the uncast (also truth) samples. This was necessary to compensate for the additional noise present in the samples with DEPc compared to the ice/air samples (ranging from 1.1 to 2.9 times larger). As the cast samples were noisier than the uncast samples, using the same σ value resulted in noticeably more jagged boundaries for the cast samples compared to the uncast. The larger σ value helped smooth the boundaries, but in doing so overblurred some areas. This is particularly important for the air bubble boundaries, which, when overblurred, become too thick to fully remove with the distance mapping threshold of 1.4. This is why the additional gradient thresholding step was necessary.

The component labeling step was found to be more effective for transition voxels associated with smaller air bubbles and those farther from ice structures. This is because only disconnected structures larger than one half the volume of an optically equivalent sphere were removed by component labeling; larger bubbles more easily exceed this threshold and transition voxels near ice structures are more likely to remain connected after distance mapping and gradient thresholding. If these first two steps (distance mapping and gradient thresholding) remove more transition voxels and break up the air bubble shells into smaller pieces, the component labeling will subsequently remove more voxels. In this way, the performance of the distance mapping and gradient thresholding impacts the performance of the component labeling. It should be noted that the calculation of the optical diameter was based on the truth and uncast images for the synthetic and cast samples, respectively. For a real segmentation, where only a cast image is available, the optical diameter can be calculated using the segmented image after gradient thresholding. For the samples tested in this work, the component labeling was not shown to have a large influence on the final SSA, density or correlation length. Other component labeling techniques (Schleef and others, 2014b) take only the largest component, which assumes that the entire ice structure is connected.

The calculation time of the proposed segmentation method was approximately 15 min for a 600 × 600 × 600 voxel volume on the machine listed in the Methods section. For a 1 m snow profile using a 600 × 600 pixel cross-section, the technique by Heggli and others (2009) takes approximately 1 month (including scan time). The proposed method would take approximately 24 h for the image processing only. Scan time for the proposed method is approximately 5–10 d depending on the concentrations used and scan parameters. In all of the μCT methods discussed here, the acquisition of μCT scans is either the rate limiting step or is of the same order of magnitude as the image processing itself. Therefore, future improvements for obtaining snow microstructures will come from using faster μCT devices or alternative measurement techniques.

The use of other segmentation methods (e.g. (Hagenmuller and others, 2013) could be used in place of the proposed method, but a comprehensive comparison of the large number of methods was outside the scope of this paper.

4.4. Error threshold

The 10% error threshold used to assess the accuracy of the casting and segmentation method is based on errors reported in the literature. For SSA, errors range from about 10 to 15% for optical methods (Matzl and Schneebeli, 2006; Gallet and others, 2009; Arnaud and others, 2011; Montpetit and others, 2012), while methane adsorption results in errors around 12% (Legagneux and others, 2002). For density, density cutting is known to provide errors of 6.2 to 11% (Conger and McClung, 2009), while Heggli and others (2009) reported errors up to 10%. For correlation length, Proksch and others reported errors of 16.4% for HRP when comparing to μCT (Proksch and others, 2015). The method by Heggli and others (2009) provides an excellent comparison for this work due to the similarity of the casting method and the method discussed here. Heggli and others (2009) reported errors up to 7 and 10% for SSA and density, respectively. It has also been shown that differences in SSA and density of up to 5% are obtained with different calculation methods, even when the same binary image is used (Hagenmuller and others, 2016). We therefore believe that the 10% error threshold used here is a reasonable benchmark for SSA, density and correlation length.

4.5. Synthetic images

For SSA, the large errors associated with large ASAV values, as well as the direction change from low ASAV values to high ASAV values, can likely be at least partially attributed to the fixed gradient threshold of 99.5% of the CDF. At small ASAV values, the gradient threshold incorrectly removes real ice voxels because the fraction of high-gradient transition voxels is limited by the amount of air/DEPc interface. As more air is added (i.e. increasing ASAV), the number of transition voxels removed increases, while the number of incorrectly removed ice voxels decreases. This process tends to decrease the SSA. After a certain point, there are more transition voxels present than can be removed, so the SSA increases with further increasing ASAV. Using a dynamic gradient threshold could therefore increase the range of ASAV values for which the method provides accurate results.

Compared to SSA, the density is less sensitive to ASAV, as expected. This is likely because, for the air bubble distributions chosen in this work, the total number of transition voxels is relatively small compared to the total number of ice voxels. Since the density is only dependent on the total number of ice voxels, there are not enough transition voxels to induce large density errors except for the very highest ASAV values. The trend of increasing density with increasing ASAV is again likely due to the gradient threshold. At low ASAV values, more ice voxels are erroneously removed and the density is lower. At high ASAV values, there are more transition voxels than the gradient threshold can remove, so more voxels remain in the final structure and the density is higher. It is unclear why 417 m⁻¹ yields such substantial outliers in this trend for S1_I, S1_B and S4_B.

Compared to SSA and density, the correlation length is the least sensitive to ASAV. The trends for increasing or decreasing correlation length with increasing ASAV do not seem to follow any clear pattern. Specifically, it is unclear why S1_I, S1_B and S5_B behave differently than the other samples. S5_B is a weak layer formed of large plate-like grains with significant void space and is the lowest density sample tested. Additionally, this sample volume was slightly smaller than the others (551 × 551 × 501 voxels instead of 600 × 600 × 600 voxels) due to partial collapse of the weak layer during sample preparation. The combination of small sample volume with low density might exacerbate errors by the segmentation method and contribute to the observed behavior. S1_I and S1_B are the most dense samples, but it is unclear why this would lead to the described behavior.

4.6. Cast images

4.6.1. SSA

For SSA, all cast samples except S5_B have errors below 10%. The large SSA errors associated with the S5_B sample could be due to the smaller sample volume and low density, which would make the sample more sensitive to small shifts during the casting process. This could potentially be corrected for with image registration as performed in Heggli and others (2009), but was not done here. The calculated ASAV values for all the samples with an error of less than 10% are below the ASAV threshold of 333 m⁻¹ determined from the synthetic analysis, except for S3_B. While the calculated ASAV for S3_B (384 m⁻¹) is higher than the nominal ASAV threshold (333 m⁻¹), it is still consistent with the synthetic data because we cannot determine exactly where the 10% threshold falls between 333 m⁻¹ and the next ASAV value of 417 m⁻¹.

For SSA, all samples except S4_B are above the 1:1 line. This indicates that the casting and segmentation method systematically overestimated the SSA for all samples except S4_B. S4_B may be an outlier, or there may be an SSA dependence on the performance of the segmentation method. The synthetic data supports the idea of an SSA dependence since all 0 m⁻¹ ASAV points of samples with truth SSA values below 10 m² kg⁻¹ are found above the 1:1 line, while the 0 m⁻¹ ASAV points of samples with larger SSAs are found below. An SSA dependence could be explained by differences in how the Gaussian filter affects low and high SSA samples; the filter is more likely to oversmooth and eliminate smaller features on the thinner structures found in the high SSA samples compared to the larger, rounder grains in the low SSA samples. However, the 10 m² kg⁻¹ threshold identified in the synthetic data is not consistent with the cast data, where S3_B is found above the 1:1 line with an uncast SSA of 21.6 m² kg⁻¹. It is difficult to determine the cause of this apparent inconsistency because there is only one point in this SSA range.

Finally, S1_I and S1_B were the only samples that were inconsistent when comparing the cast data to the synthetic data at the closest nominal ASAV. The SSA of both samples was overestimated in the cast samples, while it was underestimated in the synthetic data. This is likely due to the small SSA errors for these samples. As the points are already close to the 1:1 line, even small variations in the calculated SSA can move them from one side to the other. The differences between the cast SSAs and synthetic SSAs at the matching ASAVs are 1.2 and 2.5% for S1_I and S1_B, respectively, while the absolute errors of the cast and synthetic SSAs with respect to the truth values are less than 1% for both samples. Additionally, the nominal ASAV does not represent the true ASAV of the synthetic data, so there could be additional effects due to deviations from the nominal ASAV.

4.6.2. Density

For density, all cast samples except S4_B have errors below 10%. This is consistent with the synthetic data, where density errors less than 10% were observed for nominal ASAV values less than or equal to 833 m⁻¹ (again excluding S4_B). It is unclear why S4_B performs well for SSA, but poorly for density. The large error associated with S4_B cannot be purely based on density because it has neither the highest nor the lowest density. S4_B does, however, have the highest SSA. As described above for SSA, it is possible that the oversmoothing of high SSA structures leads to an increase in the number of ice phase voxels, which artificially increases the density. Additionally, visual comparison of the truth and synthetic segmented images for S4_B show that the ice structure boundaries are noticeably more jagged in the synthetic image compared to the truth image. Therefore, it is possible that the artificially increased density (through overblurring and feature loss as described above) leads to large changes in density, but, for the SSA calculation, is counteracted by the additional surface area produced by the jagged boundaries.

For density, all samples with uncast densities above approximately 254 kg m⁻³ are below the 1:1 line, while samples with lower densities are above. This is somewhat supported by the synthetic data, where the 0 m⁻¹ ASAV points for all samples except S4_B are under the 1:1 line. This means that only one sample, S5_B, is inconsistent between the cast and synthetic data with respect to under- and overestimation. At the equivalent ASAV point, the density of S5_B is overestimated in the cast data and underestimated in the synthetic data. As SSA is an intrinsic property and does not significantly change if the total amount of ice in the sample (i.e. density) changes, shifts in the structure during casting are not expected to dramatically change the SSA, but would have effects on the density. This is exactly what is observed; S5_B is consistent for SSA, but not density.

4.6.3. Correlation length

For correlation length, all samples except S4_B have errors less than 10%. This is consistent with the synthetic data, where correlation length errors less than 10% were observed for all nominal ASAV values (up to 3333 m⁻¹), again excluding S4_B. This is similar to the results for density and it is likely that the factors described above for the poor performance S4_B with respect to density are also responsible for the poor performance with respect to correlation length.

For correlation length, all samples except S5_B are above the 1:1 line, which is supported by the location of the 0 m⁻¹ ASAV points in the synthetic data; all 0 m⁻¹ ASAV points (including S5_B) in the synthetic data are above the 1:1 line. The reason for the inconsistency of S5_B between the cast and synthetic data is likely due to the small sample size and shifts during casting. Since there was relatively little ice in the uncast volume (i.e. low density) to begin with, even small changes could cause large differences between the uncast and cast samples.

5. Conclusions

This work demonstrates a casting method using contrast-enhanced DEP for μCT of snow. It is shown that both BaTiO₃ and I ₂ can be used to enhance the contrast of DEP to provide sufficient contrast with ice for direct processing of μCT images, and both are safe for field use. A three-phase segmentation technique using traditional Gaussian smoothing followed by a three-step method for handling transition voxels is shown to be fast (15 min for 600 × 600 × 600 voxel volumes), partially unsupervised and able to accurately calculate the SSA, density and correlation length of various snow types with errors of less than 10% for ASAV values below 333 m⁻¹. The presented method is a significant improvement on previous snow casting methods for large-scale sampling of snow.