Spectral Albedo of Dusty Martian H2O Snow and Ice

Recent evidence of exposed H2O ice on Mars suggests that this ice was deposited as dusty (<∼1% dust) snow. This dusty snow is thought to have been deposited and subsequently buried over the last few million years. On Earth, freshly fallen snow metamorphoses with time into firn and, if deep enough, into glacier ice. While spectral measurements of martian ice have been made, no model of the spectral albedo of dusty martian firn or glacier ice exists at present. Accounting for dust and snow metamorphism is important because both factors reduce the albedo of snow and ice by large amounts. However, the dust content and physical properties of martian H2O ice are poorly constrained. Here, we present a model of the spectral albedo of H2O snow and ice on Mars, which is based on validated terrestrial models. We find that small amounts (<1%) of martian dust can lower the albedo of H2O ice at visible wavelengths from ∼1.0 to ∼0.1. Additionally, our model indicates that dusty (>0.01% dust) firn and glacier ice have a lower albedo than pure dust, making them difficult to distinguish in visible or near‐infrared images commonly used to detect H2O ice on Mars. Observations of excess ice at the Phoenix landing site are matched by 350‐μm snow grains with 0.015% dust, indicating that the snow has not yet metamorphosed into glacier ice. Our model results can be used to characterize orbital observations of martian H2O ice and refine climate‐model predictions of ice stability.

determinant of albedo. Coarser snow grains reduce the number of scattering events per unit pathlength, thereby increasing the likelihood of photon absorption, and reducing the albedo. Snow that has metamorphosed into glacier ice absorbs up to five times more solar radiation than snow, which can lead to subsurface melting within the ice, as is seen on Earth (Brandt & Warren, 1993;Liston & Winther, 2005).
Recent observations have indicated that dusty H 2 O ice is currently being exposed at martian mid-latitudes (e.g., Byrne et al., 2009;Dundas et al., 2018;). The precise nature (grain size and dust content) of this exposed ice is currently uncertain, and estimates of the time scale of martian snow metamorphism range from decades/centuries (Clow, 1987;Kieffer, 1990) to millions of years (Bramson et al., 2017). While the spectral albedo of dusty martian H 2 O snow has been calculated for some cases (Clow, 1987;Cull et al., 2010;Gyalay et al., 2019;Kieffer, 1990;Singh et al., 2018), no model of dusty martian firn or glacier ice (i.e., H 2 O ice with small air content) exists at present. Modeling the spectral albedo of dusty snow, firn, and glacier ice will allow for detailed characterization of the observed exposures of dusty martian ice, and help improve global climate models (GCMs) that are currently unable to completely replicate the spatial distribution of ice deposits on Mars (Haberle et al., 2017;Madeleine et al., 2014;Naar et al., 2020).
In this study, we model the spectral albedo of dusty martian H 2 O snow and ice across the solar spectrum (0.3-3 μm), and compare our results with in-situ measurements of martian ice at the Phoenix landing site Smith et al., 2009).

Methods
We incorporate the spectral absorption properties of martian dust (Wolff et al., 2009) into models that have been used for mixtures of dust with terrestrial snow and glacier ice. Refractive index data for water-ice are taken from Warren and Brandt (2008). The densities of pure H 2 O ice and martian dust are assumed to be 917 and 1,300 kg/m 3 , respectively (Arvidson et al., 2004;Moore & Jakosky, 1989;Moore et al., 1999).
To calculate the albedos of dusty martian snow, we use Mie theory and the delta-Eddington method (Dang et al., 2015;Joseph et al., 1976;. For dusty martian firn or glacier ice, we use a "specular delta-Eddington" model that was initially developed by Mullen and Warren (1988) for modeling the albedo of terrestrial lake ice. Subsequently, it was modified for firn or glacier ice, and validated against spectral measurements of terrestrial firn and glacier ice, some of which contained volcanic ash (Dadic et al., 2013). This model incorporates an optional specular flat surface overlying a scattering layer whose properties are calculated using the delta-Eddington method. The overlying specular surface occurs in glacier ice, and to some extent in firn, but not in snow, leading to the authors using a "specularity parameter" s, which varies from 0 for an optically rough, scattering medium such as snow, to 1 for a flat, smooth surface such as ice (Dadic et al., 2013). Details are provided in Appendix A.
We model a variety of snow and ice grain sizes (Table 1) because martian ice has been observed in-situ only at the Phoenix landing site, where no direct measurements of grain size, dust content or density were made . In order to distinguish between snow and ice, all snow densities greater than 550 kg/m 3 were modeled as firn/glacier ice rather than snow because the primary scattering mechanism within dense ice is caused by the bubbles present within the ice rather than the snow grains (Warren, 2019). The "grain sizes" listed in Table 1 for glacier ice are derived from the Specific Surface Area (SSA; the area of air-ice interfaces per unit mass of ice), following Bohren (1983), who showed that bubbly ice can be modeled as very-coarse-grained snow. Unless specified, all calculations are performed assuming a solar zenith angle of 49.5° (the effective zenith angle for diffuse incidence), and assuming that the snow/ice is thick enough to obscure any underlying material (>∼1 m; termed semi-infinite). The grain radius for dust is set to 1.8 μm for all calculations shown in Section 3, based on orbital and in-situ observations of martian dust that indicate KHULLER ET AL.   (Wolff et al., 2009) and H 2 O ice (Warren & Brandt, 2008 Wolff et al., 2009). Model results are compared with measured albedo of observed in-situ excess ice (ice exceeding the soil pore space) at the Phoenix landing site (68. 22°N, 234.25°E;Blaney et al., 2009;Smith et al., 2009). Note that ice within soil pore spaces was also observed at the Phoenix landing site Mellon et al., 2009;Smith et al., 2009), but in this work we are focusing on albedos of excess ice.

Grain Size Effects
Figure 2 illustrates how larger snow grains increase the likelihood of photon absorption (Dang et al., 2015), which leads to a reduction in albedo. There is a noticeable transition between snow albedos and firn/glacier ice albedos. In particular, the absorption features at 1.5 and 2 μm caused by overtones and combinations of fundamental vibrational modes become damped. At these two wavelengths the albedo is very small; the glacier ice has higher albedo than coarse-grained snow because of specular reflection.

Martian Dust Effects
Martian dust is 3-9 orders of magnitude more strongly absorbing than H 2 O ice at wavelengths 0.2-1.0 μm ( Figure 1). Thus, small amounts of dust significantly reduce the single-scattering albedo  of H 2 O snow and ice.  is a dimensionless parameter that ranges from 0 to 1:  = 0 for pure absorption and  = 1 for pure scattering.
The size of a dust grain (set to 1.8 μm in Figures 3 and 4) affects both its absorption and its scattering. A small dust grain has more surface area per unit mass than a large grain, so its scattering-to-mass ratio is higher. But the absorption-to-mass ratio is also higher in smaller grains, because the interior mass of a large dust grain is not exposed to radiation and is thus "wasted" for absorption (because incident photons are quickly absorbed in the outer shell of the grain and do not reach the interior). Which of these two effects dominates depends on the wavelength, and the size of the ice and dust grains, as will be discussed below. Figure 3 illustrates how the addition of small amounts of dust (10 −4 %-0.1%, i.e., 1-1,000 ppm by weight) drastically alters the single-scattering coalbedo (1 -). Addition of dust to snow or ice dramatically increases (1 -) for λ < 1.2 μm, where the imaginary index of refraction for martian dust is several orders of magnitude larger than that of H 2 O ice, causing greater absorption. Also note the general increase in (1 -) with ice grain size for the pure cases (black curves in Figure 3), because larger ice grains are more absorptive.
KHULLER ET AL. Note. The specific surface area is the area of air-ice interfaces per unit mass of ice. It is reciprocally related to grain size, as SSA = 3/(r ice × 917 kg/m 3 ). The values of SSA for glacier ice are those measured by Dadic et al. (2013). The specularity parameter for firn should be intermediate between values for snow and ice, but varying s within wide limits did not improve the agreement between modeled and observed albedo, so Dadic et al. just set it to s = 1. SSA, specific surface area.

Table 1 Model Parameters Used for Different Types of Snow and Ice, Classified Based on Their Density and Approximate Grain Size
Figure 2. Albedo variation of clean snow, firn, and glacier ice albedo with grain size. The specific surface area (SSA) used for each corresponding grain size is as follows: 100 m (SSA  32.7 m 2 /kg), 500 m (SSA  6.5 m 2 /kg), 2.5 mm (SSA  1.3 m 2 /kg) and 14 mm (SSA  0.23 m 2 /kg).
At longer wavelengths, the imaginary index of ice is more similar to that of dust, and the greater volume fraction of ice in the mixture causes the effect of even 0.1% dust to be almost negligible for snow (Figures 3a and 3b). For firn and glacier ice, the addition of dust >0.01% causes the coalbedo (1 -) to decrease at longer wavelengths, indicating an increase in scattering within the ice-dust mixture. The coalbedo is approximately proportional to rn i ; that is, the product of grain radius r and imaginary index n i (Equation 10 of Bohren & Barkstrom, 1974). Although the imaginary index of dust is 1-2 orders of magnitude larger than that of ice at these wavelengths, the size of dust particles is 3-4 orders of magnitude smaller than the ice grains, allowing the addition of dust to reduce the coalbedo. Figure 4 shows the effects of small amounts of martian dust on snow and glacier ice albedo. In fine-grained snow, dust reduces the albedo for λ < ∼1.4 μm ( Figure 4a). As the ice grain radius increases, dust causes greater reductions of albedo (Figures 4b-4d), because radiation penetrates more deeply in coarser grained snow and ice, where it encounters more absorbing material before reemerging at the surface. However, adding dust to coarser grained snow and ice causes an increase in albedo for wavelengths larger than 1.2, 1.0, and 0.9 μm in Figures 4b-4d, respectively. At visible wavelengths, successively adding more dust to firn or ice sequentially decreases the albedo, but eventually the absorption by dust saturates, and the scattering by dust then becomes dominant, so the firn or ice albedo begins to increase toward the spectrum of pure dust. This reversal happens at ∼0.1% dust (Figures 4c and 4d).
10.1029/2021JE006910 4 of 12 Figure 3. Effect of dust amounts on the single-scattering coalbedo for 100 μm, 500 μm, 2.5 mm, and 14 mm snow or ice grain sizes. By "pure" glacier ice we mean glacier ice that contains air bubbles but no dust.

Comparison With Phoenix Ice Data
At present, the only available in-situ observations of excess martian ice are from the Phoenix landing site, taken by the Phoenix Surface Stereo Imager (Lemmon et al., 2008). Orbital measurements of martian ice over the entire solar spectrum are available from CRISM and Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité (OMEGA; Bibring et al., 2005), but atmospheric effects make quantitative analysis of dust content and ice grain size challenging. Thus, while the Phoenix measurements cover only part of the solar spectrum (∼0.45-1 μm), they can be compared with our model results for dusty snow and ice.
What was measured at the Phoenix landing site was the radiance reflected by the target (ice or soil), ratioed to the radiance reflected by a calibration surface (a "perfectly diffusive reflective white surface"), and corrected for the incidence angle of the sunlight (Drube et al., 2010;Zamani et al., 2009). The Phoenix team's conversion from measured radiance to the albedo plotted in Figure 5 thus implicitly makes the default assumption that the soil or ice target is an isotropic reflector, since its actual bidirectional reflectance distribution function (BRDF) is unknown.
The soil albedo observed at the Phoenix landing site is well matched by the model result for 1.65 μm martian dust, suggesting that small amounts of dust may be present on the surfaces near the ice ( Figure 5). The KHULLER ET AL.
10.1029/2021JE006910 5 of 12 amount of dust needed to completely obscure underlying soil and ice is 0.6 mm (equivalent to 10 e-folding depths), based on the absorption coefficient of martian dust at these wavelengths. Scattering by dust would make an even thinner layer opaque, so 0.6 mm is an upper limit.
From Figure 4, it is apparent that the grain size of the Phoenix ice must be less than 2.5 mm, because the addition of small amounts of dust does not increase firn or glacier ice albedo at wavelengths near 1 μm enough to match the Phoenix ice albedo of ∼0.6 ( Figure 5). Although the solar incidence angle for the Phoenix observations (∼45°) is slightly smaller than the angle used in Figure 4 (49.5°), this difference in angle causes a negligible effect at these wavelengths .
We varied the snow grain size to produce an agreement between modeled and observed albedo near wavelengths of 1 μm, where dust amounts up to 0.1% have minimal effect. This procedure led to us obtaining a snow grain radius of 350 μm (SSA  9.3 m 2 /kg). However, at shorter wavelengths, all grain sizes for clean snow produce albedos that exceed the Phoenix observations. To match the Phoenix observations, the inclusion of 0.015% dust (1.3 μm grain size) was required ( Figure 5). This result constrains the output from previous modeling work by Cull et al. (2010), who reported that this relatively friable ice contained up to 1% dust. Unfortunately, Cull et al. (2010) were unable to derive grain-size estimates because the bidirectional reflectance model they used would have required observationally determining the scattering properties (asymmetry parameter) of the ice to model the ice-dust mixture, which was not done at the Phoenix site. More recently, Gyalay et al. (2019) used a third, separate technique to estimate that the ice content was up to 93% by volume, or 87% by weight, based on their assumed densities for soil and ice. That result was obtained for a modeled ice grain size of ∼1 mm, with 60 μm soil grains (Gyalay et al., 2019). However, their estimated ice content is likely to be too low, because dust contents greater than 1% make dusty snow and ice indistinguishable from pure dust at visible wavelengths (see Figure S2 in .

Detecting H 2 O Ice Using HiRISE Images
Data with high spectral resolution but relatively coarse spatial resolution (∼18 m/pixel) are available from the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM; Murchie et al., 2007). For finer spatial resolution, recent attempts to identify H 2 O ice on Mars (e.g., Byrne et al., 2009;Dundas et al., 2018;Dundas et al., 2021) have used relative color variations within images formed from three bands in the visible and near-infrared (VNIR), with pixel size ∼25 cm, from the High Resolution Imaging Science Experiment (HiRISE; McEwen et al., 2007). The mapping of wavelength to false color is 0.874 μm (red), 0.694 μm (green), and 0.536 μm (blue). In stretched false-color images of likely ice-hosting features on Mars, materials considered to be ice are those that appear relatively "blue" (i.e., having relatively higher values in HiRISE's 0.536 μm filter).
However, our modeling results indicate that snow metamorphism and the addition of small amounts of dust drastically alter the albedo of H 2 O ice at the HiRISE filter wavelengths. Thus, it is possible for ice to appear relatively yellow/white in stretched, false-color HiRISE images even though it contains less than 1% dust . Dusty ice can also resemble nearby lithic material at visible wavelengths, when present as firn or glacier ice with 0.01%-0.1% dust (Figures 4c and 4d), potentially like some of the recently discovered scarps thought to be ice-rich (e.g., Figure 2 of Dundas et al., 2021). Greater amounts of dust (>1%) will cause snow and ice to be indistinguishable at these wavelengths.

Origin of Ice at the Phoenix Landing Site
A dusty-snow origin for the excess ice at the Phoenix landing site is questioned primarily because of the presence of a few pebbles overlying the ice Sizemore et al., 2015). While a few models have been proposed to explain the presence of excess ice in the near-surface (Fisher, 2005;Sizemore et al., 2015), those models are difficult to reconcile with decameter-thick ice deposits being discovered throughout the mid-latitudes of Mars (Byrne et al., 2009;Dundas et al., 2018;. These recent discoveries indicate that it is more likely that the majority of shallow ice on Mars was deposited as dusty snow, which was subsequently buried through lag buildup and aeolian activity. Snow burial likely caused dust, pebbles and boulders to be incorporated within and on top of the ice, as is seen on Earth (Boulton, 1978). Alternatively, lithics within and on top of the ice might also have been emplaced by impact ejecta and/or meteorites.
The optical properties of the Phoenix ice are also consistent with formation by dusty snowfall. Ice formed through vapor diffusion into soil pore spaces would imply a far greater soil:ice ratio (>∼50% soil), which is not consistent with the light-toned ice's high visible albedo. Additionally, ice in the form of lenses (Fisher, 2005;Sizemore et al., 2015) is typically dense, and therefore will not have enough scatterers (bubbles) to explain the Phoenix ice's relatively high albedo. Thus, the Phoenix ice likely formed by snowfall with small amounts of dust in the snow, which has metamorphosed over time into 350 μm-sized snow grains.
A 350-μm grain size for the ice exposed at the Phoenix landing site is smaller than the 700-800 μm grain size derived for H 2 O ice at the northern residual polar cap . However, the H 2 O ice at the cap experiences seasonal deposition and sublimation of H 2 O and CO 2 frost, which can mobilize dust and alter local thermal and atmospheric conditions (Kieffer, 1990;Langevin et al., 2005). These frost processes make it difficult to constrain the age and rate of metamorphism of the cap H 2 O ice. In contrast, the Phoenix ice was buried under a few centimeters of soil and dust. Radar detections suggest that the Phoenix ice might extend to depths of 9-66 m (Putzig et al., 2014). Assuming an average depth (or equivalently, subsurface ice thickness) of 40 m, this amount of ice could have been deposited as dusty snow over the last 1 Myr (Madeleine et al., 2014).
If the Phoenix ice originated as snow long ago, how did it evolve over this time? The cross-sectional area of snow grains (or equivalently, the snow grain size squared) grows linearly with time (Gow, 1969;Linow et al., 2012;Stephenson, 1967). Thus, if the age of the Phoenix ice is 1 Myr and the initial snow grain size was 50 μm, the rate of snow metamorphism for the Phoenix ice is 0.12 μm 2 /yr, which is an order of magnitude smaller than the rate of snow metamorphism we derive using estimates from Bramson et al. (2017). Bramson et al. (2017) obtained a value of 86 Myr for 550 kg/m 3 firn to metamorphose into 917 kg/m 3 density glacier ice. Based on Table 1, these density values correspond to grain sizes of 500 μm and 16 mm. Assuming these grain sizes implies that the rate of snow metamorphism is 16 0 5 86 10 3 2 2 6 2 mm mm / years m Note that these predicted rates of snow metamorphism can vary greatly with local atmospheric and thermal conditions. For example, ice that might have experienced melting ) is likely to have experienced far greater rates of metamorphism, because even slight melting results in a rapid increase in grain size (Warren, 2019).

Summary
We have modeled the effects of dust and snow metamorphism on the albedo of martian H 2 O ice. Both effects generally cause a reduction in albedo, thereby increasing the energy absorbed by the ice. For example, dust can reduce the albedo of fine-grained snow (100 μm) from its pure-snow value of 0.9-1.0, down to 0.4-0.7, for wavelengths short of 1 μm. Similarly, the addition of small amounts of dust (10 −4 %-0.1%) reduces the visible albedo of bubbly H 2 O ice to as low as ∼0.1, depending on the wavelength and the grain size of the ice. In some cases, the albedo of dusty firn and glacier ice can be lower than that of pure dust alone. At long wavelengths, 0.1% dust can increase the albedo of coarse-grained snow, firn and glacier ice (Figure 4). Accounting for these radiative effects is crucial to our understanding and detection of exposed H 2 O ice on Mars. 0.01%-0.1% dust can render firn and glacier ice indistinguishable from pure dust at HiRISE filter KHULLER ET AL.
10.1029/2021JE006910 7 of 12 wavelengths. Thus, dusty snow and ice might not necessarily appear relatively blue in stretched, false-color images (as is commonly assumed for martian H 2 O ice).
Our model matches measurements of the excess ice observed at the Phoenix landing site, using 350 μm snow with 0.015% dust. This result suggests that the snow has not yet metamorphosed into glacier ice. If the results of our modeling are incorporated into GCMs, simulations of the distribution of martian H 2 O ice may be improved. In addition, results can be compared with currently ongoing studies of atmospherically corrected orbital measurements of H 2 O ice using hyperspectral CRISM data (e.g., Pascuzzo et al., 2019).

A1. Albedo of Pure Snow and Dust
For snow and dust, we use a Mie scattering code (Mätzler, 2002) to calculate single-scattering quantities separately, based on terrestrial models for snow (Dang et al., 2015;. Scattering by snow and dust grains is assumed to be similar to spheres in their far fields, with radii ice r and dust r , respectively, in meters. Each particle has volume V r  4 3 3 /  and cross-sectional area 2 A r   (i.e., ice dust ice , , , V V A and dust A ). Based on the formulas for scattering by spheres, the complex refractive index   m  and the dimensionless size parameter x are required as inputs to the Mie scattering code: where  is the wavelength, and r is the average grain radius of the material. We then calculate the following quantities for ice and dust, separately: where ext  is the extinction cross section, ext Q is the dimensionless extinction efficiency,  is the single-scattering albedo, scat  is the scattering cross section, g is the asymmetry factor and  is the scattering angle (see van de Hulst (1957) for a review of these terms). The absorption and scattering efficiencies abs Q and sca Q are also calculated for each material. Note that we average over a small range of grain radii (r r  /10) to eliminate the ripple generated by Mie theory for a monodispersion .
Using the delta-Eddington approximation for a semi-infinite material, we obtain the following transformations for g and : The spectral albedo, A of pure snow and dust is given by (for the semi-infinite case):

A2. Albedo of Pure Firn and Glacier Ice
The spectral albedo of firn and glacier ice accounts for specular reflection at the surface. Based on terrestrial models for firn and glacier ice, the albedo A is given by: where s is the "specularity parameter," which specifies the degree to which the upper surface behaves specularly; s ranges from 0 to 1. 1 R is the external reflection coefficient (air to ice), whereas 2 R is the internal reflection coefficient (ice to air). inc A and trans A are the delta-Eddington albedos calculated for incident and transmitted light at angles and t  , respectively. Scattering is assumed to take place only due to air bubbles within the ice. Snell's law gives: where rat m is the ratio of the real components of the refractive indices of the transmitted media to the incident media (i.e., air: 1 and ice: ice n ).
This is the same as Equation 5 of Dadic et al. (2013). Adding sca  and abs  gives the extinction coefficient The asymmetry parameter g is taken from Mullen and Warren (1988). The single-scattering albedo  and the optical depth  are then given by     sca ext / and    z/ ext , where z is the thickness of the ice. z ≥ 20 m is assumed to be optically semi-infinite, so z = 20 m is assumed for all cases.

A3. Albedo of Dusty Snow, Firn, and Glacier Ice
For snow and dust, using the absorption and scattering efficiencies abs Q and sca Q and the areas of the particles, the absorption and scattering cross-sections can be calculated: The spectral albedo of dusty snow and glacier ice can then be found using Equation A7 and A8, respectively, with the mixture optical properties derived in Equations A28-A30.