A novel method to visualize liquid distribution in snow: superimposition of MRI and X-ray CT images

Abstract

The relationship between the behavior of water in snow and its microstructure is crucial to improve the prediction of wet snow disasters. X-ray computed tomography (X-ray CT) is frequently used to observe snow microscopically. However, distinguishing between ice and water in the X-ray images is difficult because ice exhibits an X-ray absorption coefficient similar to that of water. In contrast, magnetic resonance imaging (MRI) acquires nuclear magnetic resonance (NMR) signals of protons in a liquid and visualizes the NMR signal intensity, enabling discrimination between water and ice signals. However, snow grains and pore spaces cannot be distinguished in MRI images because they do not generate NMR signals.

To investigate the relationship between the microstructure of snow and the distribution of liquids in snow, we developed a novel method that combines X-ray CT and MRI images to compensate for the disadvantages associated with each method. Using this method, we successfully visualized where liquid (C₁₂H₂₄) occupied pore spaces. We also showed the possibility of using C₁₂H₂₄ instead of water to obtain water retention curve of snow cover, which is a fundamental aspect of hydraulic properties. Although there is room for improvement in the visualization of water in snow, such as shortening the imaging time to escape snow metamorphosis and image superimposition methods, this method is expected to effectively elucidate the behavior of water in snow and clarify the characteristics of wet snow.

1. Introduction

An accurate description of meltwater movement in snow cover is essential to improve the prediction of wet snow disasters. Percolating water occasionally ponds at the boundary between snow layers with different characteristics caused by capillary barriers (e.g. Khire and others, 2000; Waldner and others, 2004; Avanzi and others, 2016). Water ponding above a capillary barrier is prone to have an unstable flow, subsequently resulting in a preferential flow path. Consequently, heterogeneous water discharge occurs from the snow (e.g. Marsh and Woo, 1984, 1985; Nomura, 1994; Yamaguchi and others, 2018). This heterogeneous water discharge affects the timing of full-depth avalanches because it is occasionally triggered by water reaching the bottom of the snow cover. Moreover, snow grain characteristics change rapidly under wet conditions, and the snowpack mechanical properties, stability, and optical properties alter with liquid water content (Marshall and others, 1999; Baggi and Schweizer, 2008; Mitterer and others, 2011; Techel and Pielmeier, 2011; Dietz and others, 2012; Mitterer and Schweizer, 2013; Schmid and others, 2015). Therefore, knowledge of the distribution of liquid water content in snow cover resulting from water movement is also crucial for understanding snow properties.

Two crucial physical quantities are considered in the study of the movement of water in snow. The first is the water retention curve (WRC), which is also known as the ‘water content–matric potential curve’, or ‘capillary pressure–saturation relationship’, which shows the relationship between the liquid water content and matric potential, and is a fundamental aspect of hydraulic properties in porous media (Klute, 1986). The second is the saturated hydraulic conductivity (Ks), which indicates the ability of porous materials such as soil and rock to transmit water. For example, previous studies (Yamaguchi and others, 2010, 2012) have measured and analyzed the snow WRC using bulk methods. Based on their data, they proposed models to characterize the WRC of snow in terms of grain size and density. To study Ks using bulk methods, certain studies (Shimizu, 1970; Calonne and others, 2012) have measured the air permeability of snow and then modeled the Ks of snow in terms of grain size and density. Combining these models with the bulk method has successfully reproduced the 1-D water movement in the snow to a certain extent (Hirashima and others, 2010; Wever and others, 2014, 2015). However, water in natural snow infiltrates selectively and not uniformly, resulting in a nonuniform distribution of water in the snow. Consequently, runoff from the bottom of the snowpack occasionally exhibits heterogeneity (Yamaguchi and others, 2018). To describe these complex water movements in snow, Hirashima and others (2014) developed a multidimensional water transport model and successfully reproduced heterogeneous water infiltration in the snow. Their model is innovative because it represents water movement in snow in 3D; however, the equations used in the model are based on the bulk method, which does not directly account for the effects of snow microstructure and porosity. As an example of the direct observation of detailed water movement in snow, Walter and others (2013) measured pore-scale water flow through snow using Fluorescent Particle Tracking Velocimetry. Their results indicated that water in snow exhibits highly complex movements and can flow upward against gravity, even when saturated. Katsushima and others (2020) showed that the finger flow occurring during penetration into dry snow is a crucial factor in channel development and that the penetration process is strongly influenced by differences in the vertical microporous structure. Therefore, understanding water movement in a natural snowpack requires a microscale perspective that clarifies both the pore structure and the layer structure. To advance our understanding in this field, methods need to be developed to directly measure the relationship between water movement and snow microstructure. Recently, several novel methods, such as using Raman spectroscopy (Maggiore and others, 2022) and near-infrared hyperspectral reflectance (Donahue and others, 2022), have been applied to measure detailed water distribution in snow; however, these methods can not directly obtain information on 3D snow microstructure properties, such as pore size distribution, that control water movement.

There are two methods available for examining snow structure at high resolution. One method uses X-ray computed tomography (X-ray CT), and the other uses magnetic resonance imaging (MRI). X-ray CT has a very high spatial resolution, ranging from tens to several micrometers, and is applicable in various fields, including the human body and electronic components. Therefore, it is the most versatile, nondestructive and high-resolution method for measuring 3D structures. Several studies have been conducted on snow imaging, analysis and the developing of experimental tools using X-ray CT (e.g. Coleou and others, 2003; Schneebeli, 2004; Heggli and others., 2011; Hagenmuller and others, 2013; Wiese and Shneebeli, 2017), and various studies have discussed snow characteristics from a microscopic perspective (e.g. Kaempfer and Schneebeli, 2007; Löwe and others, 2011; Calonne and others, 2012; Pinzer and others, 2012). However, these studies were limited to dry snow because the difference in CT values between water and ice is marginal, and the boundary between water and ice is submerged in noise and cannot be clearly distinguished. In contrast, MRI is a method of acquiring 3D internal information regarding an imaging target using nuclear magnetic resonance (NMR), which is a phenomenon that interacts with the electromagnetic waves of a frequency specific to atomic nuclei in a static magnetic field; the NMR signal of hydrogen atoms in liquid water has been acquired and imaged. Therefore, when a wet snow sample is imaged using MRI, no NMR signal is acquired from snow grains (solid water) or pore spaces, and the NMR signal is detected from only the liquid water contained in the sample. Recently several wet snow studies have been conducted using MRI (Adachi and others, 2009, 2017, 2019, 2020; Katsushima and others, 2020).

Considering the strengths of X-ray CT and MRI, we developed a novel method to directly observe pore spaces that are filled with liquid. The advantage of our new method is that we can directly measure the relationship between the hydraulic properties and the microstructure at a high spatial resolution. In this study, we first introduce the methodology of our novel method and discuss the measurement errors at this stage. Next, we present preliminary results, including actual measurements using our method, followed by the discussion of the results. Finally, we discuss strategies for improving the utilization of water in our methods.

2. Methods

2.1 Devices

In this study, we used two nondestructive 3D imaging techniques. One was an X-ray CT system used to analyze the microstructure of a snow sample before wetting. The other was an MRI to analyze the distribution of liquid in snow after wetting. The X-ray CT system was μCT 35 (SCANCO Medical, Switzerland), and its sample size depended on its resolution; the maximum scan size encompassed a height of 120 mm and a diameter of 37.9 mm. The resolution varied from 3.5 to 72 μm. The MRI in this study was specifically developed for measurements at low temperatures with a static magnetic field intensity of 1.5 T produced by a permanent magnet of sintered NdFeB (Adachi and others, 2017, 2019). The primary component of the MRI system, including the permanent magnet, was located in a cold room. Generally, the static magnetic field of a permanent magnet is temperature dependent, which becomes nonuniform at low temperatures. This phenomenon causes distortion in the captured images. To avoid these effects, the measuring section containing the permanent magnet was kept strictly at 26 °C, where the static magnetic field uniformity area of the permanent magnet was optimized using heaters controlled by a thermostat. In contrast, the temperature of the area housing the measurement samples is maintained at the same temperature as that of the cold room, at approximately 0 ℃, through the circulation of air from the cold room (Adachi and others, 2017). The imaging space of the MRI was a 50 mm sphere, and the image resolution can range from 50 to 400 μm depending on the image time. In this study, we employed the spin echo (SE) imaging method, one of the predominant imaging methods in MRI studies (Bernstein and others, 2004), to analyze MRI images. Although SE imaging methods can display tissue anatomy with a high resolution, they require multiple pulses and echoes to generate a signal, resulting in relatively protracted imaging times.

2.2 Samples case and snow samples

To capture the same sample for X-ray CT and MRI, an identical sample case must be used to prevent sample breakage. However, a typical sample case for X-ray CT cannot be used for MRI imaging because certain parts consist of metal, which must not be used in MRI because of magnetism. Therefore, we developed a sample case comprised of acrylic, considering only its use within the strong magnetic field of the MRI (Fig. 1a). The case length for entering the sample was 75 mm, and the outer and inner diameters were 26 and 24 mm, respectively. Due to this sample case, the sample could be imaged using both MRI and X-ray CT. An additional feature of the sample case was that it had a pipe installed in its interior walls (Fig. 1b). The outer and inner diameters of the pipe were 3 and 1.8 mm, respectively. The pipe was used as a guide to introduce liquid into the sample and served as a landmark when superimposing X-ray CT and MRI images. The inner diameter of the pipe was insufficient to introduce the liquid directly due to interference from capillary forces. Therefore, we used a syringe equipped with a needle featuring an outer diameter of 1.27 mm and a length of 80 mm. For injecting the liquid into the snow sample, we installed the syringe needle in the pipe until its tip reached slightly above the bottom of the holder, and we then slowly injected the liquid (Fig. 1c). The amount of liquid to be injected varied depending on the sample. This was because the water head, which was target of our observation, required adjustment to align with the MRI field of view.

Figure 1. Images of the sample case. (a) Image of the sample case from the side. (b) Image of the sample case from above. (c) Schematic diagram illustrating the process of injecting liquid injection into a sample.

Snow samples were prepared using snow stored in a cold room (−15 ℃) for several months, and their grain size was controlled using sieves. Two snow samples were prepared: large grains (L-sample) and small grains (S-sample). The average grain size and bulk density of each snow sample were calculated from X-ray CT image analyses. In this study, the grain size was defined as the sphere-equivalent diameter obtained from the volume of each grain. The grain size and bulk density of the L-sample were 1.4 (SD = ± 0.7) mm and 409 kg m⁻³, whereas those of the S-sample were 1.2 (SD = ± 0.2) mm and 456 kg m⁻³, respectively.

2.3 Image procedures

The sample case was first filled with sieved snow in the cold room at −15 ℃, where the X-ray CT was located. The microstructure of the samples (snow grains and pore spaces) was imaged using the X-ray CT (Fig. 2a). X-ray CT images the sample at a higher resolution than MRI does, but the resolutions must be identical to superimpose the two images. Therefore, the image resolution of the X-ray CT and the MRI was fixed at 72 μm to align with the resolution of the MRI. The total imaging time of the sample using X-ray CT at this resolution was 2 h. This study used Otsu binarization (Otsu, 1979), a standard image analyses method, to binarize the X-ray CT images (snow grains and pore spaces).

Figure 2. Every image obtained during the imaging process. Red broken lines in the figures indicate the location of the pipe. (a) Images obtained using X-ray CT. White areas represent snow grains, and black areas indicate the pore spaces. (b) Images obtained using MRI. Blue areas indicate liquid (C₁₂H₂₆), and white areas indicate snow grains and pore spaces. (c) Image-superimposed MRI and X-ray CT images. White areas represent snow grains, black areas indicate the pore spaces, and blue areas indicate the liquid (C₁₂H₂₆). (d) Vertical cross section of the snow sample obtained from the reconstructed 3D image. Black areas indicate snow grains, white areas indicate the pore spaces, and blue areas indicate the liquid (C₁₂H₂₆).

After X-ray CT imaging, the sample case with the snow sample was moved to the cold room at −5 ℃, where the MRI was located. The snow case with the snow sample remained in the cold room until a constant sample temperature was reached under the cold room temperature (−5 ℃). As MRI imaging is time-consuming at the present stage with a resolution of 72 μm, n-Dodecane (C₁₂H₂₆) with a melting point of −9.7 ℃ was used instead of water to prevent the risk of snow metamorphism during the imaging. C₁₂H₂₄ was injected into the sample from the pipe using the syringe (Fig. 1c), and the C₁₂H₂₄ then permeated the sample owing to the capillary force of snow. An MRI image was subsequently captured, which provided C₁₂H₂₄ distribution information, under −5 ℃ (Fig. 2b). The total imaging duration of the sample using the MRI with a resolution of 72 μm was 8 h. Otsu binarization (Otsu, 1979) was also used to binarize each image in MRI images (C₁₂H₂₄ and others, including ice and pore space).

MRI and X-ray CT images at the same height (height resolution: 72 μm) were manually superimposed using ImageJ software (Collins, 2007) based on the landmark of the pipe (Fig. 2c). Similar analyses were performed at each height, and the superimposed images were reconstructed to obtain 3D images to cross- and vertical-section images (Fig. 2d).

2.4 Difference in physical properties between water and C₁₂h₂₄

Differences in the physical properties of water and C₁₂H₂₄ likely affect their behavior in snow. The crucial physical properties to consider when assessing the movement of liquids in snow are density, surface tension, ice contact angle and the viscosity coefficient of each liquid. The first three (density, surface tension and ice contact angle) are related to capillary forces that determine the WRC, while the viscosity coefficient is related to the velocity of liquid movement. There are several density measurements available for C₁₂H₂₄ at positive temperature (see Appendix table 2). These results show the temperature-dependent nature of C₁₂H₂₄ density (Fig. 3a). Using its temperature dependency, we estimated the density of C₁₂H₂₄ at −5 ℃ to be 769 kg m⁻³, which is approximately 23% smaller than the density of water (1000 kg m⁻³). Measurements of the surface tension of C₁₂H₂₄ at positive temperatures were also reported (see Appendix table 3). These results show that the surface tension of C₁₂H₂₄ has a nearly constant value (0.02 N m⁻¹) at the temperatures of 20 ℃–50 ℃. From these results, we estimated the surface tension of C₁₂H₂₄ at −5 ℃ to be 0.02 N m⁻¹, which is approximately one-third that of the surface tension of water (0.07 N m⁻¹). No studies have reported on the direct measurement of the contact angle of C₁₂H₂₄ on ice. Due to the lack of measured data, the contact angle of C₁₂H₂₄ with ice was assumed to be the same as that of water for simplicity in this study. However, this assumption should still be open to debate. Several studies have reported on viscosity coefficient measurements of C₁₂H₂₄ at positive temperatures (see Appendix Table 4). These results demonstrate the temperature-dependent nature of the viscosity coefficient of C₁₂H₂₄ (Fig. 3b). Using its temperature dependency, we estimated the viscosity coefficient of C₁₂H₂₄ at −5 ℃ to be approximately 2.0 × 10³  Pa s, which is of a similar magnitude to that of water at 0 ℃ (1.7 × 10³ Pa s).

Figure 3. Temperature-dependent nature of the physical properties of C₁₂H₂₄. (a) Temperature-dependent nature of the density of C₁₂H₂₄. Plotted data are shown in Appendix Table 2. The line in the figure is an approximate line represented by the formulas in the figure. (b) Temperature-dependent nature of the viscosity coefficient of C₁₂H₂₄. Plotted data are shown in Appendix Table 4. The curve in the figure is an approximate curve represented by the formulas in figure.

3. Results and discussion

3.1 Measurement error estimation

In this study, the same sample was analyzed using two images measured by entirely different methods: X-ray CT and MRI. Therefore, the size of snow grains and pore spaces may differ in each image due to differences in the measurement method. To investigate the error caused by this factor, we compared the total area of the pores measured via X-ray CT and that obtained through MRI at the same cross-sectional height. As already mentioned in the introduction, pores can be measured directly with X-ray CT but not with MRI. In this study, we selected several cross-sections near the bottom of the sample where the pores were nearly filled with C₁₂H₂₄. In calculating the total pore space using MRI, the space filled with C₁₂H₂₄ in each cross-section was considered as the total pore space. We analyzed 40 cross-sections and summarized the comparison results between X-ray CT images and MRI images (Fig. 4a). Essentially, the total pore area in each cross-section measured by both methods exhibited good agreement; however, in several cases, the total area of pores measured via X-ray CT appeared to surpass that obtained through MRI. This result indicates that the pore space determined through MRI may be underestimated. In fact, averaging over 40 cross sections, the pore space determined using the MRI images was 98 (SD = ± 3) % of the pore space determined using the X-ray images. One reason why the pore space determined using the MRI was smaller than that using X-ray CT can be attributed to the method employed for pore space estimation in this study. Filling all the pore spaces within a porous medium with liquid presents a challenge due to the entrapment of air within certain pores. In fact, field experiments conducted on soils have reported that several percentage points of the pores remained unfilled with water despite the apparent saturation of the soil (Kosugi and others, 2002). Furthermore, water content measurements have reported that, despite the snow appearing saturated with water, a percentage points of the snow pore spaces are not filled with water. (e.g. Yamaguchi and others, 2010). Therefore, it is possible that certain pores that did not contain C₁₂H₂₄ were not recognized as pores in the MRI image. Therefore, it is logical to conclude that the total pore space calculated using MRI images was slightly smaller than that using X-ray CT images. On the other hand, in ten cross sections, the total pore space calculated using MRI image was slightly larger than the that using X-ray CT images. This may be due to the effect of a slight change in pore distribution during the transfer of the sample cases from X-ray CT to MRI, but the difference is only 2% of the total pore area for each cross-section. Based on these reasons, we concluded that the effects of different measurement methods on the error do not have a significant impact on the analyses.

Figure 4. Analysis of the effects of discrepancies between X-ray CT and MRI images. (a) Comparison results of pore spaces between images of X-ray CT and MRI images. Dashed lines show the 1:1 line. (b) Schematic diagram of the method for analyzing image mismatches between X-ray CT and MRI images. White areas are snow grain shapes determined through MRI, while the areas surrounded by orange lines are snow grain shapes determined using X-ray CT. Red arrows show examples of areas where images are mismatched.

The originality of our method lies in the superimposition of MRI and X-ray CT images. During image superimposition, we manually adjusted both images based on the position of the pipe. However, if the grains are not perfectly aligned with each other, the distribution of liquid-filled pores may be misaligned, which could affect the measurement results. In addition, when moving the sample case from X-ray CT to MRI, the grains in the case could slightly move, which could affect the imaging results. To investigate the measurement error resulting from these factors, we compared the shape of individual snow grains determined using X-ray CT with those determined using MRI in 10 cross-sectional areas where it was presumed that the pore space had been saturated with C₁₂H₂₄. Figure 4b shows an example of the comparison of snow grain shapes using X-ray CT and those using MRI. In the figure, white areas represent snow grain shapes determined by MRI, while the areas surrounded by orange lines represent snow grain shapes determined by X-ray CT. As MRI cannot directly image snow grains, the white areas of each grain were estimated from spaces not filled with C₁₂H₂₄. Therefore, air-filled pores surrounding the grains may also be included in the white shape. To evaluate the effect of position mismatching, the positions of snow grains determined via X-ray CT were used as a reference. If the snow grain position was matched perfectly between the X-ray CT and MRI images, the areas surrounded by orange lines would be white. In practice, black areas are discernible within the areas enclosed by the orange line, as indicated by the red arrows in Figure 4b. This indicates that areas that should be determined as snow grain areas in MRI images are instead determined as pores due to reasons such as the misalignment of the grain. Therefore, the unmatched ratio was calculated as the ratio of black areas in the area surrounded by orange lines. We investigated the unmatched ratio for grains in the ten cross sectional areas; on average, the unmatched area accounted for 17 (SD ± 15) % of the area across the grains in these ten cross sections. The current unmatched ratio may prove too substantial to overlook in the context of a detailed analysis. Therefore, there is potential for improvement in the matching of X-ray CT and MRI images. Specifically, it is necessary to develop algorithms that automatically perform positioning using image processing technology rather than performing positioning manually.

3.2 Visualization of liquid distribution within the microstructure of snow

Figure 5 shows the cross sectional views of the L-sample. Each cross-sectional view was obtained at a different height from the bottom of the sample (Fig. 6). The white, black and blue areas in the figures indicate snow grains, pore spaces and liquid (C₁₂H₂₄), respectively. Very few black areas are in the cross section near the bottom (height: 3.24 mm) (Fig. 5a). This indicates that the pore spaces were nearly entirely filled with liquid. In the cross sectional views, up to a certain height (Fig. 5a, b, c), the pore spaces surrounding the snow grains increased, but most of the pore space remained filled with liquids. As C₁₂H₂₄ was injected directly into the sample case in this study, defining a free water surface was impossible. Therefore, in a strict sense, the data obtained in this study cannot be called a WRC. However, this area may be considered equivalent to the air entry suction area identified in the WRC of the wetting process (boundary wetting curve). The black area rapidly increased because the liquid could not penetrate the large pore spaces (Fig. 5d, e, f). This trend indicates that smaller pore spaces exhibit greater capillary forces. Only small pore spaces between snow grains indicated the presence of liquid (Fig. 5g), and ultimately, all pore spaces were devoid of liquid (Fig. 5h).

Figure 5. Cross sectional views of the L-sample. White areas represent snow grains, black areas indicate the pore spaces and blue areas indicate the liquid (C₁₂H₂₆). (a) Height is 3.24 mm. (b) Height is 8.28 mm. (c) Height is 14.18 mm. (d) Height is 15.48 mm. (e) Height is 16.20 mm. (f) Height is 16.92 mm. (g) Height is 17.64 mm. (h) Height is 20.45 mm.

Figure 6. Vertical section views of the L-sample (distance of 12.53 mm). Distance is from the left side of the cross-sectional view. The left axis shows the height position of the cross sections in Figure 5.

Figure 6 shows the vertical cross-sectional view of the L-sample. In the lower-height region, as shown in the cross sectional view (Fig 5a, b, and c), small pore spaces were observed surrounding certain snow grains; however, the pore spaces were nearly entirely filled with liquid. From the vertical cross sections, the liquid in the snow was essentially connected from the bottom to the top, although some connections twisted and turned because of the influence of snow grains and capillary forces. The upper level of the liquid height (hydraulic head) was not constant across all vertical cross-sections. This can be attributed to the nonuniformity of the pore space. The cross sectional view (Fig. 5) reveals various pore space sizes at the same height. These results demonstrate the significance of the influence of snow microstructure inhomogeneities on the behavior of liquids in snow.

3.3 Analyses of the relationship between size distribution of pore space and liquid distribution

We investigated the influence of pore space distribution on liquid behavior in snow (Figs. 5, and 6). The shapes of the pore spaces were highly complex, with the pore spaces in contact and connected to each other (Figs. 5 and 6). They were developed in 3D. Therefore, segmenting each individual pore space was necessary. In this study, we adopted 3D watershed plugins in Classic Watershed of ImageJ/Fiji to determine the boundaries between the pore spaces. The concept of a Classic Watershed in the ImageJ/Fiji plugin is to perform watershed segmentation of grayscale 2D/3D images using the flood simulation described by Soille and Vincent (1990). Figure 7 shows an example of pore space analysis using the 3D watershed scheme. First, the 3D distribution of the pore space was reconstructed based on the X-ray images (Fig. 7a), and each pore space was then segmented using a 3D watershed scheme (Fig. 7b). Figure 7b suggests that the shape of the pore space segmented by the 3D watershed scheme is not simple, but exhibits various complex shapes. Therefore, the pore space segmented using the 3D watershed scheme was considered to be more natural than considering pores as simple tubes. Based on the determination of each pore space shape, the pore space in each X-ray CT cross-sectional image was separated (Fig. 8a). The MRI images were then superimposed on the X-ray images (Fig. 8b), as described in Section 2.3. Using the images in Figure 8a and b, the distribution of pore space areas and the ratio occupied by liquid for each pore space in each cross section were calculated. The walls of the sample case may have affected the behavior of the liquid in the sample. To avoid these influences, we only analyzed data in the center area of the sample with 11.52 × 11.52 mm (160 × 160 pixels). In the analysis, only pores with size (>0.003 mm²) were used in order to avoid the effects of image mismatch discussed in Section 3.1.

Figure 7. Example of pore space analyses using a 3D watershed scheme. (a) 3D image of the reconstructed pore space based on X-ray CT images. (b) 3D image of segmented pore spaces with a 3D watershed scheme. Each color represents an individual pore space.

Figure 8. Images used in the analyses. (a) Image of the segmented pore spaces. Black areas indicate snow grains, and white areas indicate the pore spaces. The black lines in the figure indicate the boundaries between the pore spaces determined by the 3D watershed. (b) Superimposed MRI and segmented pore-space images. Black areas indicate snow grains, white areas indicate the pore spaces, and gray areas indicate liquids.

Figures 9a and 9d show the analysis results, including that the characteristics of the pore space distribution averaged over the entire cross section differed between the L-sample and S-sample. The L-sample exhibited a moderate peak, and the pore spaces were distributed over a wide range (Fig. 9a), whereas the S-sample exhibited a sharp peak, and the pore spaces were concentrated in a narrow range (Fig. 9d). The grain size of the L-sample was larger than that of the S-sample while the density of the S-sample was higher than that of the L-sample. Yamaguchi and others (2012) indicated that the smaller grain size, the more uniform the pore size if they have similar density, and larger the density, the more uniform the pore size if they have an identical grain size. Therefore, the small grain size and high density contributed to the sharp peak in the pore space, as in the S-sample. In addition, the larger SD of grain size in the L-sample also affects the larger spread of pore space area. Considering the liquid occupancy according to the size of each pore space, the smaller the pore space, the higher the degree of liquid occupation up to a certain height range (Fig. 9b and 9e). However, above a certain height, the liquid occupancy decreased sharply from the pore spaces with a larger area. Figures 9c and 9f show the ratio of the pore space occupied by the liquid to the total pore space at each height, and their trend can be considered similar to the trend in the volumetric liquid content. As previously mentioned, the data obtained in this study is not a true WRC of the wetting process (boundary wetting curve) because the free water surface could not be determined. However, the trends in Figure 9c and 9f are similar to the measured WRC (boundary wetting curve) of snow using water (Adachi and others, 2020); that is, the nearly constant value areas are maintained at a certain height, and their values subsequently rapidly decrease with height above a certain point.

Figure 9. Analyses results. (a and d) Distribution of the area size of each pore space. The x-axis represents the pore space area, and the y-axis represents the ratio of the number of pore-space sizes divided into each class to the total number. a shows data of L-sample and d shows data of S-sample. (b and e) Distribution of the ratio of occupied pore space size by liquid divided into each class by its total number. The x-axis represents the pore space, and the y-axis represents the height of the cross-section from the bottom of the sample case (height resolution is 0.072 mm). The ratio of the pore space size of the liquid divided into each class to its total number is shown in the blue scale in the legend. b shows data of L-sample and e shows data of S-sample. (c and f) Ratio of pore spaces occupied by liquid to the total area. The x-axis is the ratio of the pore spaces occupied by liquid to the total area, and the y-axis is the height of the cross-section from the bottom of the sample case (height resolution is 0.072 mm). c shows data of L-sample and f shows data of S-sample.

3.4 Differences between WRCs of water and C₁₂h₂₄

As shown in Section 2.4, there are differences between the physical properties of C₁₂H₂₄ and water. We qualitatively investigated the influence of these differences on the behavior of the WRC. Adachi and others (2020) examined the boundary wetting curve of snow (the grain size was 1.2 mm and the density was 492 kg m⁻³), similar to the snow sample called the S-sample in this study. Figure 10 shows the comparison results between the results of Adachi and others (2020) and the S-sample. To obtain the WRC of the S-sample, we estimated the liquid content using the porosity calculated from density and the ratio of the pore space occupied by the liquid to the total pore space at each height in Figure 9f. In the figure, we also fitted all data with the van Genuchten model (van Genuchten, 1980) using the Retention Curve (RETC) software (van Genuchten and others, 1994). The VG model is as follows:

(1)$$S_e = \displaystyle{{\theta -\theta _r} \over {\theta _{s-\theta _r}}} = ( {1 + {\lfloor{\alpha h} \rfloor }^n} ) ^{{-}m}$$

where S_e is the effective water content, θ is the volumetric water content, θ_r and θ_s are the residual and saturated volumetric water contents, respectively, and h is the suction. α, n and m are parameters that affect the shape of the WRC. In this study, we assumed that m = 1–1/n as referred to by Yamaguchi and others (2010, 2012). As already mentioned, it is not possible to determine the free water surface in this experiment. Therefore, the free water surface was defined as the bottom of the sample case for simplicity of the analysis. Figure 10 indicates large differences in the WRC behavior between water and C₁₂H₂₄: The height of air entry suction of C₁₂H₂₄ is substantially smaller than that of water, and the change ratio of θ with the h of C₁₂H₂₄ is substantially larger than that of water.

Figure 10. Comparison of the WRC between the data of Adachi and others (2020) and the S-sample. The blue circles are data from the S sample, and the black triangles are data from Adachi and others (2020). Blue and black broken lines are the fitting curves of the VG model with each parameter in Table 1.

Table 1. Values of α and n in the VG model in Figure 10

To investigate the cause of the difference between the WRCs of water and C₁₂H₂₄, we focused on the difference in h, which is related to the behavior of the WRC. h is described by the Laplace equation as follows:

(2)$$h = {-}\displaystyle{{2\gamma \ \cos \beta } \over {\rho gr}}$$

where, r is the capillary tube radius, ρ is the density, γ is the surface tension, and β is the contact angle between liquid and ice particles. The values of the h of water and the suction of C₁₂H₂₄ are denoted by h_w and h_c, respectively. Assuming that the contact angles of water and C₁₂H₂₄ are identical, the ratio (R_h) of h_c to h_w (h_c/h_w) is shown as follows based on Eqn (2):

(3)$$R_h = \displaystyle{{h_c} \over {h_w}} = \displaystyle{{\gamma _c} \over {\gamma _w}} \times \displaystyle{{\rho _w} \over {\rho _c}} $$

where the subscript c indicates the value of the C₁₂H₂₄, and the subscript w indicates the value of the water. The calculated value of R_h was 0.36, which shows that the h of C₁₂H₂₄ was only 0.36-fold stronger than that of water under the same r. Table 1 lists the obtained values of α and n in the VG model for data of Adachi and others (2020) and the S-sample used in this study. The values of α by Adachi and others (2020) and the S-sample are shown as α_w and α_c, respectively. The ratio (R_α) of α_c to α_w (α_c/α_w) was calculated to be 2.6, a value of the same magnitude as the inverse of R_h (2.8). Generally, α in the VG model is related to the inverse of the air entry suction. Thus, our result, that the difference in α values between water and C₁₂H₂₄ can be explained by the difference in the strength of h, is logical. But the influence of the contact angle of R_h is particularly significant. In this analysis, the contact angle of C₁₂H₂₄ with ice was assumed to be the same as water, but the contact angle of C₁₂H₂₄ with ice should be properly measured for quantitative discussion. Conversely, for n, which is related to the slope of θ vs h in the VG model, we fail to explain why the value of n was different for water and C₁₂H₂₄. Therefore, to convert the WRC of C₁₂H₂₄ to that of water, additional experiments are required to quantitatively explain the difference in n between C₁₂H₂₄ and water.

4. Possibility of imaging with water based on this method

As this study aimed to establish a method for superimposing X-ray CT and MRI, C₁₂H₂₄ was utilized instead of water as the liquid to avoid affecting the snow metamorphism. In contrast, as shown in Section 3.4, the results of this study may not be directly applicable to the behavior of water in snow. Therefore, to achieve the final objective of this study, namely, to understand the influence of snow microstructure on the behavior of water in snow, there are two approaches: The first is to understand how the different physical properties of water and C₁₂H₂₄ affect the behavior of liquids in snow; based on this knowledge, we continue to use C₁₂H₂₄ instead of using water. The second is the use of water as a liquid. Hereafter, we discuss the possibility of using water in our system.

When water is used, it is essential that the sample environment is stabilized at 0 °C during the imaging time. To maintain a stable temperature of 0 °C, we are currently developing a new sample case that will enable conventional sample cases to be placed inside a larger sample case to provide additional space between walls. By filling the wall-to-wall space with ice and water, it is anticipated that a temperature of 0 °C can be maintained throughout the duration of image times (Fig. 11). MRI imaging time with C₁₂H₂₄ required 8 h for the resolution used in this study, and a similar or even longer imaging time would be required if water was used instead of C₁₂H₂₄. Therefore, snow metamorphosis caused by extended MRI imaging time is expected to be a major issue with the current method, even if the temperature is kept at 0 °C.

Figure 11. Image and schematic diagram of the new sample case that we are developing. Case 2 is the same sample case shown in Figure 1 and case 2 is one size larger than case 1.

There are two possible approaches to shorten the MRI imaging time. The first is to decrease the longitudinal relaxation time (T1). T1 is the time required for a substance to return from its excited state to its original state. Therefore, the longer the T1, the longer the imaging time (Bernstein and others, 2004). T1 for pure water was approximately 1200 ms (the T1 of C₁₂H₂₄ was approximately 600 ms) under the condition used in this study (1.5 T, ０℃). Therefore, to obtain sufficient proton density-weighted images with the SE method used for MRI imaging, several hours of imaging time were required with the same resolution used in this study. One possible approach to shorten the T1 is to use a low-concentration (5 mmol/L) aqueous copper sulfate (CuSO₄) solution (Katsushima and others, 2020). This shortens T1 by approximately 100 ms, which is 1/12 that of pure water. Therefore, the imaging time will become 1/12 of the imaging time required when using pure water. It is assumed that the presence of low concentrations of copper sulfate does not significantly affect the behavior of water in snow; however, this effect may need to be considered when the data are analyzed in detail.

The second method involves adopting a slightly coarser image resolution in MRI compared to that utilized in X-ray CT. Currently, the resolution of X-rays and MRI is unified at 72 μm; if the imaging resolution of MRI is doubled to 144 μm, the imaging time can be decreased by 1/4 of its current duration. In this case, however, the question remains whether two images with different imaging resolutions can be adequately superimposed and whether the relationship between pore spaces and the presence of water can be sufficiently analyzed.

These two methods have their advantages and disadvantages, and certain issues must be resolved. To achieve our ultimate goal of understanding the influence of snow microstructure on the behavior of water in snow, we are studying the above two methods and combining them appropriately to decrease the time required for imaging, which can avoid snow metamorphism.

5. Conclusions

We established a novel method to visualize the liquid distribution in snow by superimposing MRI and X-ray CT images. X-ray CT can visualize the microstructure of snow but cannot discriminate water in snow. In contrast, MRI can visualize water in snow but cannot discriminate snow microstructures, such as snow grains and pore spaces. Therefore, our innovative method successfully combines the advantages of X-ray CT and MRI. Using our method, the relationships between the microstructure of snow and the liquid (C₁₂H₂₄) distribution in snow were visualized. For example, it is possible to analyze information on the spatial distribution of pore spaces and to determine the extent to which liquid occupies pore spaces of varying sizes in a 3D context.

Although we have established a novel method, there remains potential for improvement in experiments using water as the liquid instead of C₁₂H₂₄. Specifically, the MRI imaging time should be decreased such that the snow does not metamorphose during the imaging. To solve this problem, we are currently researching several solutions, including decreasing the longitudinal relaxation time using a copper sulfate solution and decreasing the imaging time by increasing the resolution roughness. Through these studies, we hope to establish a method to clarify the relationship between the snow microstructure and water behavior in snow, and to better understand heterogeneous water movement in snow.