Quantitative three-dimensional ice roughness from scanning electron microscopy

2.1 SEM Imaging of Ice

A Hitachi S-3400N VPSEM (henceforth scanning electron microscopy, “SEM”) equipped with a backscattered electron detector and a Deben Ultra-Cool stage MK3 version Peltier cooling element is used to collect scanning electron micrographs, using a protocol similar to that described by the authors in previous papers [e.g., Pfalzgraff et al., 2010; Neshyba et al., 2013]. In all experiments, an accelerating voltage of 17 kV and a probe current of 70 μA are used. The typical experimental procedure is as follows. The specimen stub, made of rough-cut copper, is mounted on the room temperature cooling element. A few milliliters of deionized water are frozen and cooled to −15°C in an aluminum reservoir. The reservoir is placed in the chamber, and the chamber is closed and pumped down to a nominal operating pressure of 50 Pa, which corresponds to an ice vapor equilibrium temperature of −32°C. The temperature of the Peltier cooling element is then reduced to −31°C, and the specimen stub is allowed to equilibrate with the cooling element. The temperature is then slowly lowered to −39°C at a rate of 0.5°C per minute. This slow rate of cooling prevents crystals from preferentially freezing to the cold stage background and increases the quantity of viable crystals growing on the copper stub. At this temperature, crystals grow quickly and appear to present smooth, prismatic facets. Once several suitable hexagonal crystals are located and imaged for calibration, the temperature is increased to −33°C or above. Further images of the same crystals are then acquired as they develop rough surfaces.

The SEM detector geometry is such that four backscattered electron detectors are positioned symmetrically around the electron beam source; each detector occupies a quadrant of an annular disk (Figure 1a). The internal radius is 2 mm, and the external radius is 7 mm. The backscatter detector assembly is approximately 10 mm above the sample during imaging. The source passes through the midpoint of the detectors and scatters from the substrate surface. Using Hitachi's 3-D acquisition mode, returning electrons are captured by each detector independently at an interval of approximately 4 s per image, producing four near-simultaneous images of the surface. These images consist of pixels measuring ~1 µm across (depending on magnification), whose values are given in backscatter intensity units (BIU) in the range 0 (black) to 255 (white). As can be seen by the brightness variation in Figure 1b, detectors A and C are most sensitive to variations in tilt angles in the x direction, while B and D detectors are most sensitive to variations in the y direction. The dependence of backscattered intensity as a function of facet orientation is a critical part of our development and is described in section 5. Next, we review formalism related to roughness distributions.

2.2 Roughness Distribution Analysis

A variety of analytical PDFs, all azimuthally isotropic on a two-dimensional plane, have been reported in the literature in the context of ice distortion and roughness. Those analytical distributions are parameterized in various ways, e.g., through scale and shape parameters. Our goal in this section is to develop methods by which statistical properties of an observed roughness distribution can be used to assign values to those parameters.

Following Yang and Liou [1998], Shcherbakov et al. [2006a], and Neshyba et al. [2013], we assume a surface height function, z(x,y), defined relative to a reference plane (usually a crystal facet), that describes the mesoscopic structure of the surface. With any given point on this surface we associate a small, planar microsurface. The square of the gradient of this microsurface is given by the following:

(1)

Associated with a given microsurface is therefore also a zenith angle, θ, which permits us to define a roughness variable:

(2)

The roughness variable is related to Z² according to

(3)

Z² and r are equal to zero when a microsurface is flat; otherwise, they are positive quantities. When considering an ensemble of microsurfaces, therefore, we expect statistical properties of Z² and r to be useful in characterizing the roughness of the ensemble. In terms of the roughness variable, a key measure is its mean, ⟨r⟩. In terms of Z, key measures are the standard deviation of Z² (std(Z²), here designated σ²) and the moments of Z² (designated < Z² >, < Z⁴ >, etc.), We expect all these measures to increase with increasing roughness of a given surface. Although they are distinct, one useful relationship (achieved through a Taylor expansion of equation 3) is as follows:

(4)

which is accurate to 5% for ⟨r⟩ smaller than 0.05, assuming Gaussian statistics.

We turn next to the task of relating these measures of roughness to parameters of analytical PDFs used to characterize roughness. One form is the uniform random tilt (URT) PDF, which many investigators have used to compute the radiative consequences of ice crystal roughness [e.g., Macke et al., 1996]. URT assumes a random assortment of zenith tilt angles distributed uniformly between zero and a maximum value, θ_max. Integration of equation 2 with this assumption leads to a mean roughness given by the following:

(5)

which inverts (for small ⟨r⟩) to

(6)

where t_M is the conventionally defined distortion, or roughness parameter, for URT [Macke et al., 1996]. Equation 6 enables us to parameterize the URT PDF based from the observed mean roughness measure ⟨r⟩.

A second analytical form is the Gaussian PDF (also Cox-Munk or Gram Charlier) [Cox and Munk, 1954], which forms the basis of the widely used radiative scattering database of Yang et al. [2013]. The term originates in an assumption of Gaussian statistics of Z² in each spatial dimension (i.e., separately in and ). Expressed in terms of the roughness variable r, the Gaussian distribution is written as follows:

(7)

The Gaussian PDF is thus parameterized by σ_CM (the “scale” parameter) [Yang and Liou, 1998]. In a perfect Gaussian distribution, the mean and standard deviation of Z² are identical; hence, one has a choice in parameterizing it from observed data from either < Z² > or σ². Here we choose the following:

(8)

A third analytical form is the Weibull PDF [Weibull, 1951; Dodson, 2006], widely used in interpreting the effects of ice crystal roughness appearing in cloud optical scattering experiments [Shcherbakov et al., 2006a; Baran and Labonnote, 2007]. Expressed in terms of the roughness variable r, this PDF is written as follows:

(9)

where Z² is given by equation 3.

The Weibull PDF is thus seen to have a shape parameter, η_w ≤ 1, in addition to a scale parameter, σ_w. The two Weibull parameters are related to the observed roughness measures according to [Johnson et al., 2002; Dodson, 2006]

(10)

and

(11)

Approximate inversion of these equations is accomplished by Taylor expansion of the ratio obtained from equations 10 and 11, followed by an iterative inversion to solve for η_w. The result is as follows:

(12)

which is accurate to better than 2% for η_w > 0.5. σ_W² is then obtained by a rearrangement of equation 10:

(13)

The GNBF algorithm yields surface information from crystals in the SEM as they lie, with facets oriented nearly orthogonal to the viewer, but not perfectly so. Therefore, the retrieved surfaces are rotated to coincide with a plane obtained by bilinear least squares fit to all points in a given ensemble of points. Ensembles typically span 50 × 50 pixels or more, which is judged to be large compared the roughness scale. Values of r and Z² are determined on these rotated surfaces, following which ⟨r⟩, < Z² >_, and σ² are obtained as described above. For visual analysis, the resulting r values are binned in intervals of ~0.01 and normalized. We have found that visual presentation of the resulting probability functions is more useful on a semilog plot, which calls for graphing rρ instead of ρ.

As discussed elsewhere, t_M, σ_CM, and σ_W are homologous roughness parameters in the sense that all indicate the degree of roughening or distortion [see, e.g., Neshyba et al. [2013, Figure 2] and give rise to similar light scattering effects (see the comparison in Geogdzhayev and van Diedenhoven [2016]). For this reason, we compare observed PDFs to all three analytical forms, recognizing that the Gaussian distribution is a special case of the Weibull distribution. URT and Gaussian parameters (t_M and σ_CM) are obtained using equations 6 and 8, respectively. To parameterize the Weibull parameter, σ_W, we use equation 13. To do so, however, requires prior knowledge of the Weibull shape parameter η_w. In principle one could use equation 12 for this purpose, but we have found this to be problematic: values of η_w obtained in this way are typically smaller than visual inspection of graphs of rρ would indicate. Instead, we base η_w on a qualitative, visual best fit criterion. We return to this point in section 12.

3.1 Characterization of SEM Response to Ice Surface Topography

To determine the three-dimensional structure of an object such as an ice crystal from SEM images, it is necessary to know how the local surface topography of a material relates to the backscattered electron signal recorded at detectors A–D. Based on the light scattering model presented in Blinn [1977], we predicted that this response would depend on projections and , where is a surface normal vector, points from the surface to detector I, and is the beam vector (see Figure 2). We therefore examined the dependence of backscattered intensity on these projections. For the crystal shown in the inset to Figure 3, for example, we identified prismatic facets a and b and drew projected vectors , , and (the latter corresponding to the crystallographic c axis). These vectors were then used to compute a true (three-dimensional) surface normal vector, , for each facet, by a procedure described in Appendix A. Projections and were computed for each facet/detector combination, and the backscattered intensities were recorded. This process was repeated over a series of micrographs taken in 5° increments from 0° to 15° stage tilt angles. Examination of the resulting data set showed that nearly all the variability in backscattered intensity depended on the difference between projections, . Therefore, we define a backscattered intensity response variable

(14)

and graph the resulting locus of points in Figure 3. The figure suggests a linear dependence:

(15)

where I specifies a detector (A–D). Parameters m_I and b_I are therefore empirical parameters determined for any given crystal. From a physical standpoint, b_I may be thought of as a background brightness, and m_I as a sensitivity. F_I, like c_I, is given in BIU, defined in section 2.

While the foregoing establishes the form of the backscattered intensity response function, as a practical matter we must parameterize the function for each scenario in the SEM viewing window. This is because parameters m_I and b_I vary somewhat from crystal to crystal, due to the presence of nearby crystals that influence the path of backscattered electrons as they travel from crystal to detector. It is cumbersome, however, to use the stage-rotation method described above for each new scenario. Instead, we chose crystals that exhibited three smooth faceted surfaces of known orientation and used backscattered intensities from a single stage orientation for calibration. For example, for the crystal displayed in Figure 4a, we drew projected vectors , , and and calculated surface normal vectors and of the corresponding prismatic facets. In addition, the normal vector to an adjacent pyramidal facet, designated , was obtained by rotating by 28° along . Three backscattered intensities, obtained by averaging brightness values from rectangular segments on the corresponding facets, are also computed. This procedure yields three values of backscattered intensity as a function of the response variable, s_I, for each detector, from which parameters m_I and b_I may be analyzed by a best fit least squares criterion, as shown in Figure 4b. Parameterizations for this and other crystals are tabulated in Table S1 in the supporting information.

3.2 Formulation of GNBF Inversion for Retrieving Surface Heights From SEM Micrographs

3.2.1 GNBF Inversion Equation

With a parameterized response function in hand, the next step is to formulate an algorithm to retrieve surface heights from SEM images. We seek an algorithm that yields a global solution while minimizing the effects of noise. The algorithm applied here is Gauss-Newton in a Bayesian framework (GNBF inversion), which is designed to optimize such properties. For problems that are moderately nonlinear, but which are linear in a small region about the solution, iteration can be used effectively to converge on the solution. The iterative GNBF inversion equation takes the following form [Rodgers, 2000]:

(16)

We refer the reader to Rodgers for a derivation and full description of this equation and present here only a brief overview with specifics relevant to this work.

In equation 16, the retrieved variable is the vector of surface heights, Z, which is retrieved as the change from an a priori value Z_a. The change from the a priori is the change needed to bring the calculated value (c_n) into agreement with the observed backscatter intensity (c_obs). Here c_n is a function of Z and is calculated from the forward model presented in section 5. However, agreement is only necessary to within the uncertainty in c_obs and c_n, set by the variance S_e (assuming Gaussian statistics). Furthermore, in GNBF inversion, the optimal value is found that also takes into account the variance in the a priori, S_a. Finding the optimal value subject to expected deviations in both c_obs (due to error) and Z_a (due to variability in the surface heights) helps constrain retrieved values to be realistic by essentially weighting the solution toward the a priori. For our purposes, a large variance in Z_a can be allowed such that the retrieval is mainly determined by the observation vector, c_obs. This vector has units of backscatter intensity, while Z has units of height. Thus, it is necessary to transform between them. This is accomplished using a variation matrix equation involving the kernel K, as will be shown below.

Three major reformulations are needed to apply GNBF inversion to the problem at hand. First, SEM images consist of four two-dimensional matrices of backscattered electron intensity (because there are four detectors). These data need to be multiplexed into a single, one-dimensional vector, c_obs. The x-y coordinates of these images must also be multiplexed into a single (shorter) one-dimensional vector. Third, the matrix formulation for the kernel K needs to be developed in a way that relates these observations directly to surface heights. These vector and matrix formulations are addressed below for continuity, but the reader uninterested in these details may skip directly to section 8.

3.2.2 Vector and Matrix Formulation

For clarity, we develop the required formulas here in the context of a 3 × 3 grid, such as displayed in Figure 2b, in which backscattered intensities are understood to originate from four triangular “pixels,” each received by four detectors. (An analogous one-dimensional development is given in the supporting information.) Based on this picture, surface heights are specified by an 8 × 1 matrix as follows:

(17)

Normal surface vector components (in x and y directions) are specified by 4 × 1 matrices:

(18)

which are combined into a single 8 × 1 matrix:

(19)

Observed backscattered intensities at a given detector I are given by the 4 × 1 matrix:

(20)

which are combined (using four detectors) into the 16 × 1 matrix:

(21)

Next we define a 4 × 4 diagonal matrix that contains the dependence of the backscatter response function, equation 15, on x direction gradients:

(22)

and y direction gradients

(23)

and combine them into a 4 × 8 matrix:

(24)

Variations in the observed intensity at all four detectors can now be expressed as a function of variations in the surface normal x and y components according to

(25)

Surface normal components can be obtained from surface heights according to

(26)

where M_xy is defined by

(27)

in which M_x and M_y are gradient operator matrices in the x and y directions; these are given explicitly for the one-dimensional case in the supporting information. We next shift the variation operator (δ) to the right, giving the variation matrix equation

(28)

in which we have defined

(29)

The variation matrix equation, equation 28, expresses changes in backscattered electron intensities (δc) resulting from an arbitrary variation in surface heights (δZ) as a matrix problem. It bears noting that the shifting of the variation operator leading to equation 28 is a key part of the development, as doing so ensures contiguity of the surface. It is, moreover, in the form required to proceed with GNBF inversion algorithm, with the kernel K defined by equation 29. In the example developed here (Figure 2), the number of pixel intensities (16) exceeds the number of surface heights (8). Inversion of the variation matrix equation in this case is therefore overdetermined: the number of known observations (length of c) is greater than the number of unknown heights (length of Z). Applications of GNBF are always set up in this way so as to ensure stability of the inversion.

3.2.3 Application of GNBF Inversion

To apply the GNBF algorithm, we define δZ appearing in equation 28 as

(30)

where Z_n is a previously obtained (or initial) vector of surface heights, and Z_n+1 is the result of a desired subsequent iteration. Using Z_n, we calculate c_n = F_I (Z_n) for each detector, and express the variation in c as

(31)

We construct S_a as a diagonal matrix whose elements equal the square of the estimated standard deviation in the heights, Z; we typically specify this standard deviation as ~10 µm in our retrievals. Similarly, S_e is constructed as a diagonal matrix whose elements equal the square of the estimated uncertainty in the observed backscattered intensity. We typically specify this uncertainty as ~2%. A sensitivity analysis studying the effect of varying S_a and S_e is described below, in section 12. We use a priori values Z_a = 0 and an initial solution Z_n=0 = 0. Because GNBF is applied iteratively, it is not necessary for the forward model, F_I (Z_n), to be linear, but rather only that it be weakly nonlinear and characterized by an error contour surface with a single minimum. We find that only three iterations are needed for convergence in most cases.

In practice, application of the GNBF inversion algorithm is limited by the size of the kernel, K. The number of elements in K increases as the square of the number of pixels, which itself scales as the square of the length of a side of a roughly square subset (or “panel”) of an SEM image. On a laptop computer, we find that GNBF inversion is limited to panels up to about 50 × 50 pixels. With the help of a graphical processing unit, we can increase this to panels of about 100 × 100 pixels. Analysis of larger subsets of a given SEM image is done by patching together GNBF-derived panels in a composite reconstruction. Discontinuities in these reconstructions, where panels are adjacent to one another, are therefore often evident.

GNBF inversion is validated by comparing retrieved surface angles of a smooth crystal to known crystal facet orientations. For the surface shown in Figure 5, the retrieved angle between prismatic facets Pr1 and Pr2 is 52°, a 13% error from the presumed angle of 60°.

3.3 Calculation of Roughness Statistics

In previous work [Pfalzgraff et al., 2010], the authors described a distinction between ice crystal roughness associated with ablation versus growth. Here we describe our efforts to quantify this distinction. Crystals were grown and calibrated (i.e., values of m_I and b_I were determined) at a temperature of −36°C and a chamber pressure of 50 Pa. Growth scenarios presented below were obtained by monitoring crystals over a period of a few minutes as they continued to grow at temperatures below −33°C. Ablation scenarios were similarly obtained, but by raising the temperature of the Peltier cooling element to −32°C, which is just above the equilibrium temperature, or higher. Surfaces were then retrieved using GNBF inversion using values of m_I and b_I obtained for that crystal and characterized in terms of roughness according to the methods described in section 2.

Figure 6 shows growth roughness on crystal 2016-06-30_ice4, case 4.1, after several minutes at a temperature of −33°C. An SEM image from detector A is shown in Figure 6a, in which azimuthally anisotropic roughness can be seen as vertical trenches in the image. A horizontal intersection between two prismatic facets occurs in this field of view, but it is scarcely visible by the A detector. It is better seen by the B and D detectors, because of their orientation below and above this intersection. Regions to be reconstructed are indicated by the boxed segments. The A and B detector grids shown in Figure 6b illustrate the different perspectives provided by each backscatter detector. The A detector highlights the trench-like roughness feature, whereas the B detector highlights the facet intersection. Figure 6c shows the surface heights retrieved using GNBF inversion, in which both the facet intersection and the roughness are evident. Roughness measures of the upper portion of this area (on the prismatic facet labeled “Pr1”) are Z² = 0.022, r = 0.010, and σ = 0.15. Figure 6d shows the distribution of observed roughness, along with a URT distribution parameterized according to equation 6. It is evident that URT provides a poor representation of the observed distribution. Also shown in Figure 6d is a Gaussian PDF, parameterized by equation 8, which clearly matches observations better than URT. As for the Weibull PDF, our experience has shown that it is difficult to evaluate its possible advantage over the Gaussian PDF from such plots: the lines with varying values of the Weibull shape parameter tend to lie too close together. Better discrimination is possible by visual inspection in a log space, as shown in Figure 6e. Based on this graph, we assign a value of η ≈ 0.95. Because this value is close to 1, we infer that the Weibull PDF contributes little beyond a simple Gaussian. For the same reason, as indicated by equation 8, σ_CM ≈ 0.15.

We next examine a second case of growth roughness, on crystal 2016-06-30_ice5, case 5.1, also at −33°C. Figure 7a shows an SEM image of the facet surface, with the region to be reconstructed indicated by boxes. Figure 7b shows an expanded view of the retrieved region, paired with the result of a forward model calculation based on the retrieved surface. The reconstructed surface itself is shown in Figure 7c, with the vertical scale exaggerated to highlight roughness features. Figure 7d compares experimental and Weibull PDFs along with observed values. We find roughness measures of r = 0.019 and σ = 0.22. A slight shoulder peak is evident at r = 0.065, which corresponds to a tilt angle of approximately 20°. A Gaussian distribution (η = 1) appears to fit the distribution better to the left of this shoulder, while a shape parameter of η = 0.8 does better to the right. Overall, we assign a value of η = 0.9. Figure 7e depicts the reconstructed surface with equal vertical and horizontal scales.

We present in Figure 8 our analysis of crystal 2016-06-30_ice1, roughness case 1.2.1, imaged under growth conditions at −36°C. Figure 8a suggests that a pyramidal facet located in the upper portion of the image is smooth, while the prismatic facet exhibits significant growth roughness. The reconstructed region includes the roughest part of the prismatic facet, with results shown in Figure 8b. As expected, retrieved roughness features are clearly azimuthally anisotropic, although these features are less ordered than in the previous example. Figure 8c shows the corresponding roughness PDF (prismatic facet only), characterized by r = 0.033 and σ = 0.30. The shape of this distribution is best fit by η = 0.85, with a shoulder once again evident. Figure 8d presents the three-dimensional surface with equal vertical and horizontal scales.

Figure 9 displays results for the same crystal as in Figure 8 but focusing on a different region of the prismatic facet. In terms of the symmetry of the roughening, much of the same conclusions hold for this roughening, although the depth of the roughening appears greater here (Figure 9b). As Figure 9c demonstrates, this region is distinguished by a marked bimodality in the shape of the roughness distribution appears: a pronounced shoulder peak is present at r = 0.15 equivalent to a tilt angle of approximately 30°.

Figure 10 explores roughness on a crystal that was first ablated by slowly reducing temperature to −28.5°C then regrown at −36°C. The retrieval region, highlighted in Figure 10a, is located where the pyramidal and rounded facet were evident, while the crystal was growing, but after ablation and regrowth it cannot be assigned to any particular facet category. Figure 10b shows the reconstructed surface, in which the roughness appears azimuthally isotropic. Figure 10c compares the observed PDF to Weibull functions. This case, and similar ones presented in Figures S1 and S2 in the supporting information, indicates that isotropic roughness is best described by a Gaussian distribution.

Figure 11 shows a series of surfaces captured at approximately 1°C intervals between −33°C and −28.5°C. As the temperature increases into the ablation regime, roughness features grow deeper and wider. Between −33°C and −31°C, this causes an increase in r from 0.011 to 0.031. At −29.5°C, we observe a slight reduction in roughness, to r = 0.030. Analysis of the lower facet of this crystal over approximately the same temperature range is presented in Figure S3 in the supporting information, with similar results.

4.1 Sensitivity of Retrieval to GNBF Parameters

As formulated here, S_e and S_a are specified as diagonal matrices with elements corresponding to the square of the standard deviation in BIU (S_e) and the square of the uncertainty in the a priori heights (S_a). A sensitivity study was performed to estimate the influence of the values of these elements on retrieved roughness statistics. From an SEM image, surface heights were retrieved using three different S_e matrices, based on standard deviations of 0.03, 2.2, and 7.1 BIU. Other retrievals were performed using three different S_a matrices, based on standard deviations in heights of 0.70, 3.2, and 17 µm. Across all retrievals, the maximum variation in retrieved surface heights was on the order of only 0.2 µm, and roughness statistical measures (e.g., r) agreed to two significant figures. Observations of other cases showed similar robustness to variations of S_a and S_e.

We also investigated the sensitivity of roughness to the grid size of retrieved segments. Two retrievals of an identical region on the crystal 2016-06-30_ice4 (shown in Figure 6) were performed, one using a single 90 × 90 grid and the other using nine 30 × 30 grids. The resulting mean roughness values differed on the order of 5%.

A key motivation for choosing a Bayesian formalism for inverse retrieval is that the algorithm finds the optimal solution within the solution region characterized by the uncertainties. This framework avoids the extreme sensitivity to noise that is “a common feature of exact solutions to retrieval problems” [Rodgers, 2000]. This is particularly important for acquiring roughness statistics because small-scale noise will increase r. Qualitative inspection of results (e.g., Figure 6b) shows that forward modeled images do indeed exhibit reduced small-scale pixel-scale variation in backscatter intensity compared to observations, as desired.

4.2 Variation in Response Function Parameters

As described in section 5, we calibrated the backscatter response parameters m_I and b_I for each detector and crystal and used those parameters to retrieve surface heights of those crystals after roughening. We believe that these parameters vary from detector to detector because of inherent differences in detector sensitivity. For example, detector A generally records brighter backscattered intensities than the other detectors. It is unclear why, however, the backscatter response parameters should vary from scenario to scenario (i.e., from crystal to crystal). We speculate that detector sensitivity depends on the proximity of other ice crystals, which may create a lensing effect due to local variations in water vapor concentration. To minimize this possibility, we selected relatively isolated ice crystals for our analysis.

4.3 Trends From Roughness Statistics

Roughness statistics are summarized in Table 1. The naming convention for cases is as follows: the first number refers to crystal identity, second refers to the particular roughening scenario, and the third differentiates between different analyses of the same image. Regarding the degree of roughening, we see that values of r reach as high as 0.045 (σ = 0.39). Regarding the shape parameter, we see that ablation roughening is best described by η = 0.8, azimuthally isotropic roughening by η = 0.95, and azimuthally anisotropic roughening in the growth regime by η in the range 0.75 to 0.95. As ice crystal scattering matrices also vary with varying shape parameters [Geogdzhayev and van Diedenhoven, 2016], these results suggest that in situ measurements of scattering properties and remote sensing results may contain more facet-specific information than has been previously appreciated, as well as information that permits distinction between growth and ablation conditions.

In section 3, we noted that use of equation 12 to compute η_w typically yields smaller values compared to what one would infer from visual inspection of graphical representations. A cursory comparison of the two methods indicates that the differences can be substantial. For example, Table 1 shows that crystal 2016-06-30_ice4, case 4.1 (which appears in Figure 6), was visually judged at η_w = 0.95, whereas equation 12 yields η_w = 0.88, a 7% reduction. The reason for the bias, we suspect, is that in visual estimation, one tends to emphasize values near the peak of a given PDF, whereas this does not occur when using equation 12. More work on this will be required to resolve this matter.

4.4 Relationship to Previous Results

In situ studies of natural ice crystals at South Pole Station [Shcherbakov et al., 2006a, 2006b] result in retrieved values of η_W between 0.73 and 0.77, generally lower than values obtained here. Because our study focuses primarily on prismatic facets, the disagreement may be due to effects of crystal facets not studied here, such as basal or rounded facets. We think it more likely, however, that when faced with strongly non-Weibull PDFs (such as the bimodal distribution appearing in Figure 9c), in situ retrieval algorithms compensate by assigning excessively low values of η_W. In terms of scale parameters, our results span a similar range to that reported from the above in situ studies (σ_W up to ~0.25) but are small compared to the results of remote sensing studies (σ_W up to ~0.5) [Cole et al., 2014; van Diedenhoven et al., 2014b; Baran, 2015; Hioki et al., 2016]. But our results are obtained under a rather narrow range of conditions, −39°C to −29°C and near-equilibrium humidity, whereas remote sensing studies span a much larger range. Moreover, field studies suggest that roughness scale is sensitive to cloud top temperature [van Diedenhoven et al., 2014b] and possibly humidity [Baran, 2015].

A distinct advantage of the present method over the method based on prismatic facet intersections described previously [Neshyba et al., 2013] is that it is less restrictive: one does not need to find cases in which roughness appears at these intersections nor does one need to assume that such roughening is representative of facet interiors. Indeed, Magee et al. [2014] found that it was rare to find well-resolved roughness that intersected facet edges and therefore could not obtain quantitative data on much of the roughness they observed. A second advantage is that the present method, in retrieving heights as a function of two horizontal dimensions, provides far more data and therefore greatly increases confidence in the statistical properties of roughening.

Magee et al. [2014] also present evidence that roughness occurs on scales as small as the submicron level, which is below the imaging resolution attainable by the SEM used for our study. Further study is necessary to elucidate the relative importance of mesoscale versus submicron-scale roughness in relation to optical scattering.

The method of retrieval via GNBF inversion developed in this paper can be applied to other materials for which quantitative data concerning surface structure is of interest. However, several conditions must be met by the material in question: it must be homogeneous, such that all variation in backscatter intensity is due to surface tilt, and it must be continuous so that gradients may be calculated at all points. The first condition may be circumvented by coating techniques such as sputtering with gold palladium, if the desired features are large enough that they will not be obscured by the coating.