The Microwave Snow Grain Size: A New Concept to Predict Satellite Observations Over Snow‐Covered Regions

Abstract

Satellite observations of snow‐covered regions in the microwave range have the potential to retrieve essential climate variables such as snow height. This requires a precise understanding of how microwave scattering is linked to snow microstructural properties (density, grain size, grain shape and arrangement). This link has so far relied on empirical adjustments of the theories, precluding the development of robust retrieval algorithms. Here we solve this problem by introducing a new microstructural parameter able to consistently predict scattering. This “microwave grain size” is demonstrated to be proportional to the measurable optical grain size and to a new factor describing the chord length dispersion in the microstructure, a geometrical property known as polydispersity. By assuming that the polydispersity depends on the snow grain type only, we retrieve its value for rounded and faceted grains by optimization of microwave satellite observations in 18 Antarctic sites, and for depth hoar in 86 Canadian sites using ground‐based observations. The value for the convex grains (0.6) compares favorably to the polydispersity calculated from 3D micro‐computed tomography images for alpine grains, while values for depth hoar show wider variations (1.2–1.9) and are larger in Canada than in the Alps. Nevertheless, using one value for each grain type, the microwave observations in Antarctica and in Canada can be simulated from in‐situ measurements with good accuracy with a fully physical model. These findings improve snow scattering modeling, enabling future more accurate uses of satellite observations in snow hydrological and meteorological applications.

Introduction

Snow is a random heterogeneous medium composed of ice, air and possibly water and impurities. All its physical properties depend not only on the properties of these constituent materials but also on their geometrical arrangement at the micrometer scale, the so called microstructure (Torquato, 2002). This applies in particular to the electromagnetic properties that control the propagation of waves in snow, such as the scattering and absorption coefficients. Scattering in snow is caused by the dielectric contrast between air and ice, and its amplitude highly depends on the length scales of the microstructure. The “snow grain size” is an intuitive property commonly estimated in the field (Fierz et al., 2009). However, it is loosely defined from a geometrical point of view because snow crystals often have very complex shapes, leading to imprecise and subjective measurements. Moreover this single metric is insufficient to fully describe all the length scales. Finding a rigorous mathematical representation of the microstructure and prescribing its length scales from actual measurements in the field or from snow evolution model outputs is the biggest problem to be solved for modeling interactions between snow and electromagnetic waves. This step is crucial to ultimately predict satellite observables from snow measurable quantities (Figure 1) and conversely for a more reliable use of remote sensing to retrieve snow information for hydrological, meteorological and climate applications (Helmert et al., 2018; Hirahara et al., 2020). A paramount application is the retrieval of the snow height and snow water equivalent (snow mass on the ground), a major endeavor for snow hydrology (Derksen et al., 2019; Lievens et al., 2019; Pulliainen et al., 2020; Rott et al., 2008).

Figure 1. Two modeling chains for predicting satellite observables (brightness temperature, backscattering coefficient) from observed or modeled snowpack properties, (a) using empirical approach to adjust the microstructure parameters (top pathway) and (b) using the new physical pathway (bottom) presented in this study. The gray zone highlights where understanding is lacking, orange where understanding is progressing, green where established solutions provide sufficient accuracy.

In the visible and infrared spectral range, this problem has been solved a long time ago by introducing the optical diameter d_opt. This stems from early modeling work where snow was described as a random collection of identical non-overlapping spherical particles of ice (hard spheres (HS)) (Warren & Wiscombe, 1980). The sphere diameter, called the optical diameter or optical grain size (Wiscombe & Warren, 1980), proved to be a fruitful concept to predict scattering and absorption even when considering more complex geometrical shapes (Grenfell & Warren, 1999). For any shape, d_opt can be defined from the surface area of the ice/air interface S and ice volume V using d_opt = 6S/V. The relevance of the optical diameter comes from the fact that any medium exhibits a similar scattering behavior to HS if the particles have the same S/V ratio (Grenfell & Warren, 1999), a property called S/V equivalence. This equivalence is not strict, as the scattering behavior has a small, residual dependence on particle shape (Picard, Arnaud, et al., 2009). Hence, all modern snow optical radiative transfer models use the optical diameter to predict the albedo (Domine & Shepson, 2002). Measurements of the optical diameter and of the related metric called “specific surface area” defined as SSA = 6/ρ_iced_opt, has become considerably easier in the recent decades, with a variety of available techniques based on adsorption of methane (Legagneux et al., 2002), high-resolution 3D images obtained by micro-computed tomography (Kerbrat et al., 2008), or optical reflectance (Arnaud et al., 2011; Gallet et al., 2009; Matzl & Schneebeli, 2006; Painter et al., 2007). The latter is the most convenient technique in the field, and has enabled the collection of large data sets of SSA (Vargel et al., 2020).

In the microwave range, the problem stated above remains largely unsolved since a quantity equivalent to the optical diameter is hitherto missing. A first possible approach to represent snow is considering HS with diameter d_opt, as in the visible and infrared ranges. Results show however that this approach underestimates the scattering amplitude in the microwave range (Brucker, Picard, et al., 2011). Artificially increasing the sphere diameters by a constant empirical factor was found to be an effective solution for local-scale studies (Brucker, Picard, et al., 2011; Picard et al., 2014) but requires specific regional adjustments (Roy et al., 2013; Vargel et al., 2020). This demonstrates the lack of geometrical insight which prevents generalization and application in robust satellite retrieval algorithms at large scales. A popular extension of the HS representation is the sticky hard sphere (SHS) model (Macelloni et al., 2001; Picard et al., 2014; Tsang et al., 1985) that aggregates the particles into clusters under the effect of an ad hoc attractive force (Figures 2a and 2b). Because the clusters span a wider size than their constituent particles, the scattering amplitude is enhanced with respect to randomly positioned spheres. This highlights the prominent role of the inter-particle arrangement in the microwave range as opposed to the visible and infrared ranges where only the individual particle size and shape matter. The strength of the attractive force can be adjusted by a parameter called stickiness (Tsang et al., 1985). Despite being rigorously defined in terms of the pair potential between the particles, it is impossible to directly relate this parameter to actual measurable microstructural properties.

Figure 2. Comparisons of different microstructures composed of non-overlapping disks (in 2D for illustration purpose, (a, b, and c)) and (d) the resulting normalized auto-covariance function. (d) The shape of the auto-covariance informs about (a) the size of the particles, as it decreases more sharply from the origin small disks (a and b) than for (c) large ones, and (b) has a longer or more positive tail for (b) sticky disks than for (a) more randomly positioned disks. Both increasing size and stickiness enhances scattering.

A second approach to describe snow is considering a two-phase porous random medium, which is more general than assuming a particle collection. Such a medium is completely defined by the indicator function I(r₀), which takes the value 1 when ice is present at the position r₀ in 3D space and 0 otherwise (Figure 3a). This binary representation has been successfully used to interpret data recorded with the small-angle scattering (SAS) of neutron or X rays techniques allowing the investigation of the structure of many materials (Schmidt, 1991; Svergun & Koch, 2003). The main reason for this success is that the scattering amplitude is closely related to the auto-covariance function γ(r) of the indicator function (Porod, 1951; Torquato, 2002). The auto-covariance measures how rapidly the indicator translated by some increasing distance r, I(r₀ + r), becomes different from the original indicator I(r₀) (Figures 2d and 3b). r is called range or lag. The auto-covariance captures a combination of both size and arrangement of the medium structures. A fast decrease of the covariance at the origin (r ≈ 0) indicates a medium with small structures (Figures 2a and 2d), while a slow decrease indicates large structures (Figures 2c and 2d). Because of the importance of this decrease rate, the Porod length was introduced early in the development of the SAS technique (Porod, 1951). It is defined as the inverse of the decrease rate at the origin l_p = −(γ′(0))⁻¹. An important consequence established by (Debye et al., 1957) is that l_p is mathematically and unequivocally related to the optical diameter d_opt and density ρ for any porous medium with a smooth interface, $${l}_{\mathrm{p}}=\frac{2}{3}\left(1-\rho /{\rho }_{\mathit{ice}}\right){d}_{\text{opt}}$$ (Mätzler, 2002). An equivalent expression is obtained as a function of SSA, l_p = 4(1 − ρ/ρ_ice)/SSAρ_ice. This general result is highly relevant for snow since SSA and density are measurable quantities. Objective values of l_p can thereby be obtained easily. However, this is insufficient to predict microwave scattering because the long range behavior of the auto-covariance is important but unconstrained by l_p. For instance the clustering of spheres, as obtained with the SHS model, results in a covariance function with a more positive tail (Figures 2b and 2d) than when the spheres are fully randomly positioned (Figures 2a and 2d), resulting in stronger scattering for the same sphere size.

Figure 3. 3D Microstructure of a snow sample composed of alpine rounded grains. (a) segmented micro-CT image representing ice (blue) and air (void). (b) Normalized auto-covariance of the sample.

For an isotropic medium with randomly distributed and varying shapes, so-called Debye medium (Debye et al., 1957), the auto-covariance function has a decreasing exponential form γ(r) = exp(−r/l_c). l_c is called the correlation length and is equal to the Porod length in this particular case. For this reason, these two lengths have been often confusingly called “correlation length” in the snow literature (Mätzler, 2002; Royer et al., 2017). Given that l_p, or likewise l_c, can be derived from measurements, this microstructure representation has become popular to compute snow scattering (Mätzler & Wiesmann, 1999). Unfortunately, values of l_c derived from measurements in this way tend to overestimate scattering in many cases. Empirical scaling factors α have been therefore introduced to adjust l_c = αl_p (Mätzler, 2002), which we refer to here as the scaled exponential microstructure (sEXP hereinafter). However this microstructure lacks any physical justification because the scaling factor has to be adjusted from region to region (Vargel et al., 2020), in full similarity to the empirical adjustments required for the SHS model outlined above. Other forms of auto-covariance function have been introduced for snow, in particular some based on cut-leveled Gaussian Random Fields (Ding et al., 2010; Sandells et al., 2021). However, they also introduce parameters that are difficult to relate to snow physical properties.

In summary, the snow-microwave community faces a two-fold problem: first, it is not yet established which representation of snow microstructure (sEXP, SHS, or another) is optimal for electromagnetic scattering calculations, and second, the parameters (diameter, stickiness, correlation length) to run the scattering models are often not measurable and need to be adjusted empirically. They are also usually specific to a single microstructure representation, which causes incompatibility and confusion on how to relate them to one another and to measurable quantities.

We propose an alternative approach to solve this problem without relying on empirical adjustments (Figure 1, bottom pathway). Following the concept of the optical diameter, we introduce a new metric called the microwave grain size (l_MW) which optimally predicts scattering (Section 2). In this paper, we then show that:

1. The microwave grain size is a unifying parameter because common microstructure representations for snow (sEXP, SHS and others) can be reformulated with this new parameter, the density and the Porod length. Moreover, for a given microwave grain size value, the scattering amplitude at low frequencies becomes almost independent of the microstructure representation.

2. The microwave grain size is the product of the Porod length and a new introduced factor K that describes how widely the length scales vary in the microstructure. The so-called polydispersity K carries information on the shape of the particles and their relative arrangement. Furthermore, we find here that K is fairly constant (≈15% variations) for a wide range of convex-grained snows, thus providing a physical way to estimate l_MW from l_p and thus from easily measurable quantities in the field (SSA and density).

3. Taking K as depending on the type of snow grains (among classes in a universally accepted classification, (Fierz et al., 2009)) only is an efficient way to predict microwave observation from satellites which is comprehensively demonstrated for 104 sites in Antarctica and in Canada (Section 4).

This work hence provides a better understanding of the snow microstructure and a robust way to predict microwave scattering from measurable or observable snow properties using traceable physical relationships. These new formulations are implemented in the open-source Snow Microwave Radiative Model (SMRT). These findings and this model will help to rigorously link snowpack evolution model outputs to microwave emission and backscatter models inputs, which open great perspectives to improve the retrieval of crucial variables, such as the snow height or snow water equivalent, using remote sensing.

Background The Microwave Grain Size

A natural definition of the microwave grain size follows from the Born approximation (Born, 1926; Ding et al., 2010; Mätzler, 1998), used in several scientific domains (Gille, 2000; Porod, 1951; Teubner & Strey, 1987), and relating the scattering amplitude σ_s to the auto-covariance function γ(r) of a porous isotropic medium: [Image Omitted. See PDF]here k = 2πf/c is the wavenumber (f is the wave frequency and c is the speed of light) and ɛ is an electromagnetic term depending on the density ρ and the wavenumber. It only shows small variations in the case of dry snow in the frequency range 1–100 GHz (Löwe & Picard, 2015) and is not a source of uncertainties. The microstructure information is carried by $$\tilde{\gamma }(k)$$, the 3D Fourier transform of the isotropic auto-covariance function γ(r). We introduce the microwave grain size l_MW by noting that $$\tilde{\gamma }(k)$$ has the dimension of a cubic power of a length. Taking the static limit (k = 0), a simple definition follows: [Image Omitted. See PDF]As a direct consequence of this definition, the scattering amplitude is exactly proportional to the cubic power of l_MW in the static regime, and Equation 2 remains a good approximation in the low frequency limit (1–85 GHz for most snows), as long as the k⁴ term dominates the frequency variations over the electromagnetic and microstructure terms. This implies that knowing l_MW solves the problem of the scattering calculation. Nevertheless, no method exists to obtain or measure l_MW for snow yet. In the following we show that this parameter has very relevant properties and we devise a method for estimating it from measurable quantities.

The Unifying Role of the Microwave Grain Size

We first show how the microwave grain size is related to the specific parameters of some commonly-used microstructure representations. For this, we use the definition (Equation 2) with either the integration of the real space auto-covariance function γ(r) or the Fourier space expression at the origin $$\tilde{\gamma }(0)$$. The details of the calculations are reported in Text S1 and Table S1 in Supporting Information S1, and we briefly summarize the final results here. For the scaled exponential representation it is trivial to show by integration of the exponential function (2) that: [Image Omitted. See PDF]

The microwave grain size would thus coincide with the widely used “exponential” correlation length (Mätzler, 2002) if the auto-covariance were strictly exponential. It also appears proportional to the Porod length. Similarly for the SHS microstructure, the microwave grain size can be related to the sphere radius, density and stickiness. Since l_p is a function of radius and density, the microwave grain size can also be related to the Porod length as follows: [Image Omitted. See PDF]where the function K_SHS only depends on the stickiness τ and the snow density ρ (Table S1 in Supporting Information S1). This expression shows a clear and natural separation between the size (carried by l_p) and the effect of packing of the spheres (carried by K_SHS). Another widely used auto-covariance function was proposed by Teunber and Strey (TS representation hereinafter) for microemulsions of oil in water (Teubner & Strey, 1987). It has been previously introduced for snow (Löwe et al., 2011; Picard et al., 2018; Sandells et al., 2021). Here again a simple relationship can be obtained for TS: [Image Omitted. See PDF]where we use the dimensionless parameter q = l_p/d_TS to define the TS representation instead of the repeat distance d_TS as suggested by (Ruland, 2010). In contrast with these examples, not all the microstructure representations have closed analytical form yet. This is the case of the Gaussian Random Field microstructures as defined in (Ding et al., 2010) or (Sandells et al., 2021).

These three examples highlight that l_MW can be computed for different analytical forms of the auto-covariance function and related to the specific parameters of these forms. It is not guaranteed that an analytical expression of l_MW always exists for any microstructure representation, since the integration of γ(r) may not be carried out in closed form. However, when they exist, these relationships make the different microstructure representations comparable. For instance in a study using the SHS representation, l_MW can be calculated from the radius, density and stickiness, and can then be compared to another study using the scaled exponential representation where l_MW is simply related to the scaling coefficient and the Porod length. This provides a way to re-evaluate past studies.

Noting the central role of the microwave grain size, we went a step further and inverted these relationships to obtain all specific parameters of the considered microstructure representations (sEXP, SHS and TS) as a function of the triplet microwave grain size l_MW, Porod length l_p and snow density ρ only (the equations are reported in Text S2 and Table S2 in Supporting Information S1). The fact that such a common set of variables exists is an important and new result because it provides a unified way to parametrize different microstructure representations. Furthermore, given the definition of the microwave grain size, it is guaranteed that different microstructure representations predict the same scattering amplitude in the low frequency limit when the same microwave grain size value is used as input. Only at higher frequencies some differences between the microstructure representations may appear for a given microwave grain size, but the re-parameterization in terms of common triplet remains nevertheless effective and relevant. This result renders the choice of the best snow representation a secondary problem, and conversely implies that measuring the microwave grain size or deriving its value from measurable quantities become the primary task to be solved in order to predict snow scattering in the microwave range.

The Microwave Grain Size From Measurable Quantities

To obtain the microwave grain size from measurable quantities, it is necessary to reveal its fundamental link to geometrical characteristics of the microstructure. We established here a relationship between the microwave grain size and the chord length distribution (CLD), independent of the particular choice of the auto-covariance functional form.

Chords are line segments intersecting an infinite line with the two phases of a porous medium (Figure 4a) (Torquato, 2002). The CLD of each phase is a statistical characteristics of the microstructure (Figure 4b). Intuitively, the average chord length μ₁ of the ice (and air) gives information on the size of the ice grains (and of the air pores respectively). It is related to l_p for any microstructure by: μ₁ = l_p/(1 − ϕ) where ϕ is the fractional volume of the ice phase, that is, the ratio between snow density ρ and pure ice density ρ_ice = 917 kg m⁻³ (Ruland, 2010). It can be therefore unequivocally estimated from the measurable SSA and density. The higher order moments (μ_i, i > 1) of the CLD carry information on how dispersed the chord lengths are, that is, all higher order moments have a small value only if all the chords have similar lengths, and a high value for complex shapes or a large range of sizes. This property is called polydispersity. The second moment μ₂ was previously used to parametrize a generalized version of the TS microstructure model (Ruland, 2010). An important conclusion of this study was that the scattering amplitude increases not only as a function of size (the first order moment) but also with polydispersity (the higher order moments). We adapt and generalize here this idea to a wide class of microstructures. To this end we establish a general relationship linking the microwave grain size and the chord length moments. This is achieved in two steps (details of the calculation are given in the Text S3 in Supporting Information S1), first by relating the microwave grain size to the second derivative of the Laplace transform of the auto-covariance function $$\hat{\gamma }(s)$$: [Image Omitted. See PDF]and second by using an approximation for the Laplace transform of the CLD established in a previous study (Roberts & Torquato, 1999). This approximation assumes that the pore chord distribution is exponential, a class of microstructures known as the Boolean model (Bilodeau et al., 2007). In such a model, the solid phase is built up by randomly positioning a finite set of primary shapes (e.g., spheres, cubes, polyhedra, etc., or any combination of them) in space with possible overlap. In such a model, the pore CLD is exponential if the primary shapes are all convex (Bourgeois & Lyman, 1997). We verified how close to an exponential is the pore CLD for 167 snow samples collected in the Alps from in-lab snow growth experiments (Figure S1 in Supporting Information S1) and concluded that the Boolean model applies well to snow. Note that even with convex primary shapes, the resulting microstructure has concave parts, as it is common in snow (depth hoar, grain boundaries), because overlaps are allowed in the Boolean model.

Figure 4. (a) A 2D slice of the 3D Microstructure shown in Figure 3a with examples of chords in the air (orange segments) and ice (green segments). (b) The distribution of the lengths of all the chords in the 3D samples for the air and ice.

The relationship obtained after applying these two steps links the microwave grain size to the first four moments of the ice CLD: [Image Omitted. See PDF] [Image Omitted. See PDF]here we only consider snows with ϕ < 0.5 (density less than 468 kg m⁻³) because for ϕ > 0.5 it is recommended to swap air and ice (air primary shapes in an ice background) (Dierking et al., 2012). This relationship highlights the proportionality between the microwave grain size and the Porod length through the factor K. We call this latter factor the “microwave polydispersity” because it only involves ratios between the higher order moments and the first order moment of the CLD, and thus measure the chord polydispersity. As opposed to (Ruland, 2010) we demonstrate here that the second order moment is insufficient to fully characterize the polydispersity as relevant to microwave scattering, the first four moments are all required.

Before studying this equation in its general form, it is instructive to consider the case ϕ = 0, a medium with isolated grains and very low density even though it does not apply to snow. The microwave polydispersity of such a sparse medium writes: [Image Omitted. See PDF]

It only depends on the first and fourth moments of the ice CLD which can be related to the volume V and surface area S of the particles using the Cauchy formula (Mazzolo et al., 2003), leading to: [Image Omitted. See PDF]This equation offers a practical means to compute the microwave polydispersity for any geometrical particle with known surface area and volume (when the medium is sparse). Moreover, it gives an intuitive understanding of the polydispersity by noting that the ratio $$\frac{S}{{V}^{2/3}}$$ is related to the isoperimetric shape factor f₁ = 6V/π^1/2S^3/2, a common measure of sphericity of particles (Redenbach et al., 2012). f₁ indeed takes its highest possible value for spheres and decreases with the particle elongation. The microwave polydispersity K_sparse is proportional to $${f}_{1}^{-2/3}$$, implying that spheres are the least efficient scatterers, and the scattering amplitude increases with elongation for a given Porod length. This result may explain why representing snow as non-overlapping ice spheres usually underestimates scattering and that large empirical scaling factors had to be used in the past to reconcile model simulations and observations(Brucker, Picard, et al., 2011; Picard et al., 2014; Roy et al., 2013). To conclude for sparse media, the microwave grain size can be interpreted as the product of an elongation indicator (K_sparse) and the particle size (l_p).

In the case of dense media such as snow, the polydispersity given by Equation 8 involves two additional terms in ϕ and ϕ², with a more complex combination of CLD moments. Furthermore, the second and third moments cannot be related to S and V only. Despite this complexity, the formulation provides several hints. First, it confirms the idea of (Ruland, 2010) about the influence of the chord polydispersity on scattering. Second, it shows that the polydispersity tends to decrease with increasing density (the first order term in ϕ is negative and (1 − ϕ)^−2/3 decreases with increasing ϕ) at least for moderate densities (ϕ² ≈ 0). This implies that the microwave polydispersity recovers a well known and important effect in dense packings, where the scattering amplitude of packed particles is lower than the sum of individual particle scattering (Tsang & Kong, 2001). And last, Equation 8 allows us to estimate the polydispersity value from the CLD, which itself can be obtained from micro-CT imaging of real snow. This equation hence provides a means to obtain the polydispersity and then the microwave grain size from measurable quantities.

Materials and Methods Micro-CT Data and Chord Length Distribution

The data set used here to compute CLD was first presented in (Löwe et al., 2013). It comprises 167 snow samples scanned with X-ray tomography (micro-CT), producing 3D images at a resolution ranging from 5.1 to 10.7 μm. The samples are in fact of two categories: 37 of them are individual samples collected in the Alps while the remaining was obtained by sampling at different times from 6 in-lab snow maturation experiments. These experiments differ from each other by the imposed thermal gradient conditions, from isothermal to 100 K m⁻¹. All the samples were assigned to a snow type (depth hoar, rounded grains, faceted crystals, decomposing and fragmented precipitation particles, melt forms and precipitation particles) according to the international classification of seasonal snow on the ground (Fierz et al., 2009). The data set is therefore quite heterogeneous and is not representative of any snow on Earth, but is adequate to illustrate the effect of snow types on polydispersity. The CLD of ice and air was extracted from each 3D image after binarization, by drawing lines in the vertical and two perpendicular horizontal directions as presented in (Krol & Löwe, 2016).

Snow In-Situ Measurements

In-situ measurements were collected in Antarctica and Canada to compute microwave grain size and perform the microwave simulations. In Antarctica, snow properties were measured at 18 sites (Table S3 in Supporting Information S1) over a large range of latitudes during three scientific traverses, namely Vanish (2011–2012), ASUMA (2016–2017) and EAIIST (2019–2020). Additional measurements were taken in 2011 at Dome C (Picard et al., 2014). A relatively similar protocol was applied at every site. A borehole was drilled up to a depth of typically 8 m (4.1–17.9 m). The extracted core was sliced in ≈10 cm long pieces. Snow density was obtained by measuring the diameter, height and mass of each cylindrical slice. If a slice was not cylindrical, the height was recorded, and the density was set to that of the nearest cylindrical slice. The SSA profile was measured by short-wave infrared reflectometry using the Posssum and Asssap instruments (Arnaud et al., 2011; Libois et al., 2015). On ASUMA and EAIIST, Asssap was used to record the SSA profile along each extracted snow core of 50–100 cm length. The profiles were then assembled and the small gaps between each core were filled by linear interpolation. The profile resolution is about 1 cm. At Dome C and at sites S2b and S4 on Vanish, Asssap was used to take a single record for every 10 cm slice in a cold chamber in France (Picard et al., 2014). At point S2 on Vanish, Posssum (Arnaud et al., 2011) was directly used in the borehole to record the full profile at 1 cm resolution. This profile is however short, only 4.9 m. Both instruments, Posssum and Asssap are based on the same principle and have been inter-calibrated many times. Their accuracy was estimated to 15% against independent SSA measurements (Arnaud et al., 2011). Since the measurements from a single borehole were used for each simulation that was then compared to satellite observations representative of a 12.5 km (or 25 km) wide pixel, it is expected that the intra-pixel spatial variability is a large source of uncertainties. This prevents performing a very precise site-by-site comparison between simulations and observations. The annual mean temperature at each site was measured with a Pt100 sensor at 10 or 20 m depths, after 24hr stabilization. The complete temperature profile was not recorded because of its changing nature. For this reason, the simulations are conducted with an uniform temperature equal to the temperature measured at 10 or 20 m depth, and they are compared with the annual average brightness temperature. The density and SSA profiles are also considered independent of time. This approximation is valid because due to the cold conditions and the low accumulation, the main rapid changes only occur in the topmost ≈20 cm of the snowpack, the remaining being stable over years on the Antarctic Plateau. This is applied to the low frequencies (10, 19 and 37 GHz) where the snowpack portion contributing to emitted signal is larger than about one m depth, and this is mathematically justified by the quasi-linearity of the temperature dependence in the heat equation and the radiative transfer equation in snow (Picard, Brucker, et al., 2009). In contrast at 89 GHz, because the radiation is emitted by the topmost ≈20 cm of the snowpack, and the snow properties were measured in summer, the simulations use an uniform temperature equal to the mean December-January 2 m air temperature extracted from the ERA5 reanalysis. The results are compared with the average brightness temperature over the same months.

In Canada, 86 sites (Table S4 in Supporting Information S1) were sampled over a large latitudinal range. The density, SSA and temperature profiles were measured in snowpits down to the ground as detailed in (Vargel et al., 2020). The density was measured with a density cutter and a scale. The SSA was measured on samples extracted from the pits using the IRIS instrument (Montpetit et al., 2012) based on short-wave infrared reflectometry as Asssap and Posssum. The profiles of temperature and the soil temperature were recorded for each pit.

Microwave Simulations

The SMRT (Picard et al., 2018) is used to conduct the simulations of microwave thermal emission. The model represents the snowpack as a stack of horizontal layers specified with the in-situ properties as follows. The layering is directly derived from the density profile. In Antarctica, because some profiles are too short (e.g., S2 on Vanish) with respect to the microwave penetration depth at the lowest frequency (10 GHz), the modeled snowpack is extended down to 30 m depth by repeating the lower meter of the measured profile. The SSA which is usually sampled with a higher resolution than density is averaged for each density layer. The Porod length l_p is then deduced from density and SSA. In Canada, with a snow height rarely exceeding 1.5 m, the soil is a significant contributor to the microwave signal at low frequencies (10 and 19 GHz) and certainly plays a small role at 37 GHz as well. Unfortunately the soil characteristics relevant to microwave simulations (soil permittivity, surface roughness, …) are in general difficult to measure, and were not available here. This problem was solved by optimization of the soil parameters by (Vargel et al., 2020) using the observations at low frequencies. We have taken here the soil parameters of that study without any further adjustment.

To explore the role of microstructure representation, SHS, sEXP and TS are considered for most simulations. The original version of TS is limited to K < 1 (Teubner & Strey, 1987) but has been extended as proposed by (Ruland, 2010). The latter is implemented in SMRT as detailed in the Text S2 in Supporting Information S1.

The other settings of SMRT are common to previous studies (Picard et al., 2018; Vargel et al., 2020). In short, the Improve Born Approximation is used to compute the scattering and absorption coefficients in each layer and the Discrete Ordinate method (DORT) solves the radiative transfer equation for the whole snowpack account for multiple scattering between the layers. The outputs for each site are the brightness temperature at four frequencies and at horizontal and vertical polarizations.

Microwave Observations

Microwave observations were compared to the model simulations in order to retrieve the polydispersity and assess the simulation performance. In Antarctica, the microwave brightness temperature observations at 10 19, 37 and 89 GHz were recorded by the Advanced Microwave Scanning Radiometer 2 (AMSR2) sensor onboard Japan's Global Change Observation Mission 1st - Water “SHIZUKU” (GCOM-W1) satellite. We extracted the observations at the nearest pixel of each site (Table S3 in Supporting Information S1) from the National Snow and Ice center AMSR-E/AMSR2 Unified Level 3 daily product, version 2. The product has a resolution of 25 km at 10 GHz and 12.5 km at the higher frequencies. The observations were averaged over the period 2013–2019. The typical brightness temperature accuracy is ±1.5 K.

In Canada, the observations at 86 sites (Table S4 in Supporting Information S1) were obtained with ground-based radiometers operating at the same frequencies as AMSR2 (Vargel et al., 2020) though not all the frequencies and polarizations were observed at all sites due to instrumental failure or availability. The accuracy is typically 2 K. These data mainly differ from satellite data by the small field of view of the sensor which is at the meter scale, and is coincident with the snow properties measurements.

Results Polydispersity of Snow Samples

Figure 5 shows the microwave polydispersity K for the 167 samples taken in the field in the Alps or from in-lab snow growth experiments and imaged by micro-CT. The graph distinguishes two categories of snows as a function of the grain shape, with convex grains on the one hand and depth hoar on the other hand. Convex grains include rounded and faceted grains (typical of alpine dry snow) and melt forms (occurring during melt). There are grouped together because their respective mean polydispersity is 0.72 ± 0.084 (1σ, n = 53) for rounded grains, 0.71 ± 0.073 (n = 33) for faceted grains, 0.68 ± 0.028 (n = 5) for melt forms, showing no significant differences (pair-wise Welch's t-test, p > 0.05). Depth hoar, also known as cups because of their hollow shape, features higher values 0.85 ± 0.081 (n = 62) than the other grains, with a significantly different mean (p ≪ 0.05). Meanwhile, we note that our values compare well with values (0.8–1.2) obtained for the Boreal Finnish depth hoar in a recent investigation (Leinss et al., 2020) where the empirical scaling factor of the exponential function was determined using micro-CT images (according to our Equation 3 the polydispersity K is equal to this empirical factor α).

Figure 5. Polydispersity K of convex grains (blue) and depth hoar (orange), (a) calculated for 167 snow samples using micro-computed images (blue and orange vertical bars) (b) obtained in this study from theory (vertical dotted lines) or microwave retrieval (blue vertical and orange horizontal dashed lines), (c) and derived from (Leinss et al., 2020) which use micro-computed images (horizontal solid line). The vertical lines are used for values determined with a reasonable accuracy, while the horizontal lines are used when a wider range of values is determined.

These first results obtained with micro-CT images show that K spans a relatively narrow range 0.71 ± 0.078 (1σ, n = 91) for rounded grains, faceted grains and melt forms, if compared to the ≈10-fold potential range of variation of l_p. This result suggests that when micro-CT measurements are not available, running microwave simulations with a constant value of K for these convex grain shapes could be sufficient. The next section tests this hypothesis.

Retrieval of the Polydispersity From Microwave Observations and In-Situ Data

We use the microwave grain size l_MW to predict snowpack microwave emission first in Antarctica and second in Canada. The in-situ measurements provide the profiles of SSA and density, from which l_p can be deduced without approximation. Since no coincident micro-CT measurements were taken, we deduce l_MW by assuming that the polydispersity K is constant (but unknown) for the rounded and faceted grains, the prevailing snow grain types on the Antarctic plateau. These measurements and derivatives, along with an assumption on the microstructure representation, are sufficient to fully prescribe the microstructure in every snow layer. We performed the SMRT simulations for three different microstructure representations (sEXP, SHS, and TS) parameterized with the unfiying triplet l_MW = Kl_p, l_p, and ρ. For each microstructure representation, the optimal K value was determined by minimizing the root mean square error (RMSE) calculated between the simulated and observed brightness temperatures at 19 and 37 GHz and at vertical polarization (Figure 6a). We then test the simulations with the optimal K on a wider set of frequencies (10, 19, 37, 89 GHz) and at both vertical and horizontal polarizations (Figure 7 and S2 in Supporting Information S1).

Figure 6. Root mean square error and bias between simulations and observations calculated at 19 and 37 GHz, at vertical polarization, as a function of the polydispersity value (a) applied to all Antarctic sites over the whole profiles and (b) applied to all Canadian sites for the depth hoar layer.

Figure 7. Observed (cross) and simulated (circle) brightness temperatures at four frequencies and vertical (blue) and horizontal (yellow) polarizations at 18 sites in Antarctica (sorted from the inner plateau to the coast, Table S3 in Supporting Information S1) using sticky hard spheres and the optimal polydispersity of 0.64.

The RMSE calculated at two frequencies and vertical polarization features a clear minimum, as a function of K, of 5.8 K, 5.7 and 6.2 K for sEXP, SHS and TS respectively (Figure 6a). When the simulations with the optimal K are run at the four frequencies and two polarizations, the average RMSE is 11.4 K, 11.3 K, 11.79 K for sEXP, SHS and TS respectively (Figure 7). Both results show small differences in performance between the microstructure representations. This reflects past findings where different microstructure representations have been used with equal success (Royer et al., 2017; Vargel et al., 2020). This is an expected outcome of the microwave grain size definition as discussed above. Split per frequency, the RMSE is the lowest at 37 GHz, and increases at 19, 10 and 89 GHz (Figure 4a). We attribute these variations mainly to the in-situ measurement uncertainties, and the difference of scale between the in-situ and the satellite measurements. At 10 GHz and 19 GHz, the microwave emanates from the surface to about 15–20 m and 5–10 m depth respectively, whereas the measurements were taken up to only ≈8 m on average (Table S3 in Supporting Information S1). Even though we extended the simulated snowpack downward (Section 3.2), this is approximate and may explain part of the uncertainties in the results at the two lowest frequencies. Conversely, at 89 GHz, the microwaves emanate from the top 20 cm of the snowpack. This zone was sampled for all the cores but with a vertical resolution of 10 cm that is too coarse for this high frequency. The frequency 37 GHz is optimal given our experimental sampling, with waves mainly coming from the uppermost first meter, where accurate and detailed measurements were taken at all sites. Regarding polarization, the performance is better in vertical polarization (blue in Figure 6) than in horizontal polarization (yellow), which is a classical result (Durand et al., 2008; Wójcik et al., 2008) explained by the insensitivity of the vertical polarization to the snowpack density layering. Overall the modeling errors are of the same order as in other studies where optimizations were applied (Macelloni et al., 2001; Picard et al., 2014), and the model shows excellent skills to reproduce the latitudinal gradient showing brightness temperatures increasing from the plateau to the coast. We conclude that assuming a unique constant K for rounded and faceted grains is suitable to predict the microwave signal in Antarctica, given the uncertainties in the in-situ measurements and the difference of scale between the satellite and in-situ observations.

The optimal polydispersity value is 0.63, 0.64, 0.60 for sEXP, SHS and TS respectively (Figure 4a). These values fall in the lower range of K obtained from micro-CT on the alpine rounded and faceted grains (Figure 5). This result is remarkable because the two estimates are fully independent, providing for the first time a link between the microwave-optimized scaling factor and its microstructural origin. Furthermore, the three optimal K values for the different microstructures are close to one another (within 7%) which comes from the unifying character of the microwave grain size. While this Antarctic data set provides a first confirmation that a constant K is suitable for microwave simulations, the variety of grain types is limited, only rounded and faceted grains are present on the Antarctic Plateau.

We further test our hypothesis on the Canadian data set where highly metamorphized snow is omnipresent as depth hoar. The typical eastern Canadian Arctic snowpack consists of an upper part of rounded or faceted grains overlying a bottom part of depth hoar. For the upper part, we make and test the hypothesis that the optimal K value obtained in Antarctica also applies in the Canadian environments. In contrast, to account for the particular scattering efficiency of the depth hoar in the lower part, we consider a specific value for depth hoar (K_DH). This value is obtained by optimization as done previously, by minimizing the difference between simulations and observations at 19 and 37 GHz in vertical polarization.

The RMSE calculated at two frequencies and vertical polarization shows a minimum, with values of 22.7 K, 20.7 and 21.6 K K respectively for sEXP, SHS and TS (Figure 6a). These three values are of the same order. However for SHS, the minimum is not marked and the bias never reaches 0 K, mainly due to a systematic overestimation of the simulated brightness temperature at 37 GHz. The simulations also become numerically unstable (diagonalization error in the DORT solver in SMRT, (Picard et al., 2018)) for large polydispersity at 89 GHz preventing the exploration of polydispersity values above 2.3. SHS appears to be unable to cope with high polydispersity and to produce strong enough scattering for a given sphere size. This shows the limit of the sphere model even with highly clustered particles. We conclude that SHS is unsuitable for depth hoar, in line with past studies (Löwe & Picard, 2015; Vargel et al., 2020). The results in brightness temperature (Figure 8 and S2 in Supporting Information S1) again show the good skills of the model. The simulations at 89 GHz, showing virtually no bias (e.g., 1.2 K for TS, p-value of 0.7), are very insightful because only the upper layer contributes at this high frequency and these simulations are therefore independent of the K_DH optimizisation. This confirms that the polydispersity estimated in Antarctica and used without adjustment in the upper layer here applies well to the rounded and faceted grains in Canada. In contrast, the lower frequencies are sensitive to the depth hoar layer, and do depend on the optimized K_DH value. The optimal value only weakly depends on the microstructure representation choice, 1.25 and 1.5 for sEXP and TS respectively. However, the determination is relatively imprecise as shown by the wide minima in Figure 6b. This is particularly true for TS, the RMSE changes by less than 1 K over the range 1.2–1.9. If instead of the RMSE minimum, we consider a null bias as an optimization criterion, we would obtain optimal polydispersity of 1.4 and 1.7 for sEXP and TS respectively. Despite these uncertainties, we conclude that the optimal K_DH is certainly above 1, which is significantly higher than the polydispersity of rounded and faceted grains. It is also significantly higher than the alpine depth hoar polydispersity estimated from micro-CT. This suggests that scattering by the depth hoar in Canada is much stronger than that in the Alps (for a given l_p). Observations indicate that the structure of hoar is indeed different between these regions (Domine et al., 2016; Satyawali & Schneebeli, 2010). The eastern Canadian Arctic depth hoar is often of centimeter size and is more developed due to the very strong vertical temperature gradient prevailing during the entire winter season. The high polydispersity could thus be explained by the large ratio between the micrometer scales (the steps and the thin walls of the depth hoar crystals) and the centimeter size of the crystals or even the long range organization between the crystals as in columnar depth hoar (Figure 9). In the Alps, depth hoar is often tinier and less structured because the thermal gradients are weaker and operate over a shorter period (mostly the beginning of the snow season) which justifies a smaller polydispersity.

Figure 8. Observed (cross) and simulated (circle) brightness temperatures at four frequencies and vertical (blue) and horizontal (yellow) polarizations at 86 sites in Canada using the Teubner and Strey microstructure. Sites are listed in Table S4 in Supporting Information S1.

Figure 9. Macrophotography of rounded grains and depth hoar. The white bar indicates the 1 mm scale.

Discussion and Conclusion

This study establishes a fully tractable chain of physical links to conduct simulations of microwave scattering from measurable snow physical properties (Figure 1). For each snow layer, density and SSA or optical diameter d_opt provide the Porod length l_p which is then converted to the microwave grain size l_MW by multiplication with the microwave polydispersity K. We showed that assigning a constant value to K depending on the traditional grain shape leads to satisfactory simulations. An optimal value of ≈0.6 for rounded, faceted and melt forms and 1.2–1.9 depending on the microstructure representation for depth hoar in the eastern Canadian Arctic was obtained. The confidence in these values is relatively high for the former group, composed of convex grains, because we obtained a similar estimate with two independent methods (CLD direct calculation and microwave retrieval). However, this comparison was not performed at the same site, because of the lack of coincident micro-CT and microwave observations. For depth hoar, the value is more uncertain, but it is certainly much larger than that of the rounded and faceted crystals, which can be understood by the morphological differences (Figure 9) and the possibly wider range of structure in depth hoar (e.g., soft depth hoar, indurated depth hoar from wind slabs, indurated depth hoar from melt freeze layers, columnar depth hoar) (Domine et al., 2018).

An immediate application of this new tractable chain is to perform microwave simulations with the outputs from state-of-the-art snowpack models such as CROCUS (Vionnet et al., 2012) and SNOWPACK (Lehning et al., 1999). These models predict snow evolution from timeseries of meteorological conditions. As their outputs include all the variables required by our chain (density, SSA and traditional grain shape), it becomes obsolete to rely on empirical coefficients (Brucker, Royer, et al., 2011). This achievement should increase the interest in microwave satellite observations to assess or constrain the snowpack models in the future. Our findings open new perspectives in large scale simulations of microwave signatures (Pulliainen et al., 2020) and in data assimilation of microwave observations in snow hydrological models (Durand & Margulis, 2006).

In the future, instead of relying on the traditional grain shape to infer the polydispersity value, direct and more precise values could be obtained. For in-situ surveys, polydispersity can be obtained from snow samples imaged by micro-CT, although this involves significant work. From a modeling point of view, CROCUS and SNOWPACK already have a “sphericity” prognostic variable to represent grain shape. Unfortunately the “sphericity” definition established three decades ago (Brun et al., 1989) is not compatible with the isoperimetric shape factor which we demonstrated to be equivalent to the microwave polydispersity in sparse media. More work is needed to relate these quantities. An even more advanced and promising avenue is the future snow evolution models that are expected to describe metamorphism laws more closely to the microstructure (Leinss et al., 2020). A model able to predict the evolution of the auto-covariance function or of the CLD from fundamental thermodynamic principles would indeed enable seamless predictions of the microwave polydispersity.

However, there are still some important unsolved issues. From a theoretical point of view, a better understanding of the peculiar geometrical features of the microstructure controlling the CLD is needed. Although our results for isolated convex grains are simple and intuitive (the polydispersity K is a measure of grain sphericity), the situation for dense media seems more complex. The equation (Equation 8) established to estimate the polydispersity as a function of the chord length moments gives a practical way to compute K from micro-CT, but does not reveal exactly which geometrical features of dense media control the polydispersity. The long-range order in the medium, characterizing how grains are arranged relatively to each other, is known to influence K (Chen et al., 1990; Ruland, 2010) but investigations on the order in snow microstructure is lacking. It will also be important to determine whether the polydispersity can be assumed to be constant for depth hoar crystals grown in different conditions. The range of depth hoar polydispersity estimated in the present study is about 1.2–1.9 (50% variation) from microwave and even larger when including the calculation from micro-CT. It is certainly the largest source of l_MW uncertainties considering that l_p can be derived from measurements of SSA, with 15% uncertainties, and density with 10% uncertainties. This large range is not a surprise according to our field experience. Depth hoar is certainly the snow type with the largest visual variations in crystal size, shape and order across the world. Columnar depth hoar features the largest and most organized crystals (Figure 9) and is expected to yield high polydispersity, while some depth hoar found at the bottom of the alpine snowpack is often small and random. Further investigation on depth hoar with micro-CT is required. This study also assumes an isotropic medium from the very beginning although snow geometrical properties are known to be different in the vertical and horizontal directions (Krol & Löwe, 2016). A possible approach may be to consider a different microwave grain size for each Cartesian direction (Leinss et al., 2020).

Introducing the microwave grain size l_MW and the polydispersity K provides a way to relate the different microstructure models but it does not solve the problem of choosing the most adequate microstructure representations for snow. All the representations reach approximately the same RMSE after optimization of the polydispersity, except SHS in the case of depth hoar. In light of the results, the scaled exponential may seem attractive because of its simplicity and efficiency, and the scaling factor introduced empirically in the past by fitting exponential curves to measured auto-covariance functions (Krol & Löwe, 2016; Mätzler, 2002) or by microwave optimization (Royer et al., 2017). It appears to correspond to the polydispersity K, but it is not strictly similar as fitting an exponential curve to the auto-covariance function may differ from integrating this function with Equation 6. The scaling factor of 0.75 (on average) established by (Mätzler, 2002) is close to our estimates of K for alpine snows, and the increasing trend from fresh snow and decomposing particles to faceted grains and to depth hoar (Table S5 in Supporting Information S1) has similarities with our findings (Figure 5). Despite these great advantages, the sEXP does not respect the required mathematical properties of an auto-covariance function at the origin (Torquato, 2002). The impact of this inconsistency on microwave scattering is negligible because the long-range behavior of the auto-covariance function matters most (Equation 6). However, it is important in the optical range (Krol & Löwe, 2016) making sEXP unsuitable for a unified treatment of snow in the optical and microwave ranges. The SHS representation has a fully valid auto-covariance function and yields the best performance in Antarctica, but clearly fails to represent depth hoar. Teubner– Strey seems adequate for any type of grains, even though performance in Antarctica is slightly reduced compared to SHS. To discriminate the representation performances, future work should investigate the snow microwave response at higher frequencies (e.g., 150 GHz available on the Microwave Humidity Sounder) where the microstructure details play a more prominent role. This requires higher resolution measurements of snow properties than what collected so far.

The theory developed in this study is of interest beyond snow and microwaves. It clarifies how the microstructure of a porous medium controls wave scattering, when the wavelength is larger than the grain size, independently of the constituent materials and of the wave nature. The microwave grain size and the polydispersity K as defined here are new and general metrics useful to investigate a variety of media and waves. Conversely, we expect that the tied theoretical links will help to transfer new knowledge from the materials science to snow scattering and ultimately contribute to more efficient remote sensing applications.

Acknowledgments

The European Space Agency is funding the development of the Snow Microwave Radiative Model model (Microsnow project). The traverse data were obtained through the French Agence Nationale de la Recherche (EAIIST grant ANR-16-CE01-0011, MONISNOW grant 1-JS56-005-01, ASUMA grant ANR-14-CE01-0001 ASUMA, Vanish grant ANR-07-VULN-013), the Institut Polaire Français Paul-Emile Victor (IPEV), the National Antarctic Research Program (PNRA, grant EAIIST PNRA16-00049-B). EAIIST was also supported by the BNP-Paribas Foundation through its Climate Initiative program. Technical supports during ASUMA and EAIIST were provided the French national ice core drilling program (F2G) and the EQUIPEX CLIMCOR (ANR-11-EQPX-0009-CLIMCOR). The Canadian field campaigns were supported by Natural Sciences and Engineering Research Council of Canada and Polar Knowledge Canada.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

All microwave and snow data needed to evaluate the conclusions are available in the repository https://doi.org/10.18709/perscido.2022.05.ds367 (Picard, Löwe, et al., 2022). The Snow Microwave Radiative Model model code is available from https://github.com/smrt-model/smrt (latest version) and the specific version used in this study with unified microstructure is available from https://doi.org/10.5281/zenodo.6518996 (Picard, Sandells, & Löwe, 2022). The code to produce the main result figures (Figures 5–8) is made available from https://doi.org/10.5281/zenodo.6519037 (Picard, Sandells, & Löwe, 2022).