SNICAR-AD v3: A Community Tool for Modeling Spectral Snow Albedo

Abstract. The Snow, Ice, and Aerosol Radiative (SNICAR) model has been used in various capacities over the last 15 years to model the spectral albedo of snow with light-absorbing constituents (LAC). Recent studies have extended the model to include an adding-doubling two-stream solver and representations of non-spherical ice particles, carbon dioxide snow, snow algae, and new types of mineral dust, volcanic ash, and brown carbon. New options also exist for ice refractive indices and solar zenith angle-dependent surface spectral irradiances used to derive broadband albedo. The model spectral range was also extended deeper into the ultraviolet for studies of extraterrestrial and high-altitude cryospheric surfaces. Until now, however, these improvements and capabilities have not been merged into a unified code base. Here, we document the formulation and evaluation of the publicly-available SNICAR-ADv3 source code, web-based model, and accompanying library of particle optical properties. The use of non-spherical ice grains, which scatter less strongly into the forward direction, reduce the simulated albedo perturbations from LAC by ~9–31 %, depending on which of the three available non-spherical shapes are applied. The model compares very well against measurements of snow albedo from seven studies, though key properties affecting snow albedo are not fully constrained with measurements, including ice effective grain size of the top sub-millimeter of the snowpack, mixing state of LAC with respect to ice grains, and site-specific LAC optical properties. The new default ice refractive indices produce extremely high pure snow albedo (> 0.99) in the blue and ultraviolet part of the spectrum, with such values measured so far only in Antarctica. More work is needed particularly in the representation of snow algae, including experimental verification of how different pigment expressions and algal cell concentrations affect snow albedo. Representations and measurements of the influence of liquid water on spectral snow albedo are also needed.


In 2011, a web-based 470-band implementation of SNICAR was launched for informal educational purposes. This web model was executed more than 185,000 times between June 2011 and June 2020 by users from dozens of countries and became more widely used than anticipated. During this period there have been numerous improvements to the model (e.g., Flanner et al., 2014;Singh and Flanner, 2016;Cook et al., 2017;He et al., 2018;Dang et al., 2019), but 60 they have not been assimilated into a unified code base. Furthermore, model users have requested improved documentation of the web-based model and accompanying code. Together, this motivated the creation and release of SNICARv3 in June 2020 and SNICAR-ADv3 in January 2021, along with this study. Dang et al. (2019) compared several two-stream models of snow albedo, including SNICAR, with 16-stream solutions from the Discrete Ordinates Radiative Transfer (DISORT) model. They found that the delta-Eddington adding-doubling two-stream snow algae. Sect. 4 describes modeled spectral and broadband albedo sensitivities to different environmental features, model options, and LAC concentrations. Finally, Sect. 5 presents evaluations of the model against field measurements of spectral snow albedo from different environments.  Multi-layer two-stream models utilize these fundamental quantities to solve for the upward and downward radiative fluxes at each layer interface within each spectral band, from which spectral albedo at the model top, radiative absorption within each layer, and spectral transmittance through the column are derived. Ice grains scatter strongly in the forward direction because they are generally much larger than the wavelengths of interacting light. Strong forward peaks in the scattering phase functions necessitate transformations or analytic manipulations of the two-stream input variables to accurately represent fluxes (e.g., 105 Joseph et al., 1976;Wiscombe and Warren, 1980;Bohren and Huffman, 1983). The delta scalings of τ , ω, and g (Joseph et al., 1976, Appendix A) are applied to account for forward scattering in all versions of SNICAR.
While previous versions of SNICAR adopted a tri-diagonal matrix two-stream solver (Toon et al., 1989), SNICAR-ADv3 utilizes an adding-doubling solution (Briegleb, 1992;Briegleb and Light, 2007), which has several advantages (Dang et al., 2019). First, the adding-doubling framework enables internal Fresnel layers to be incorporated into the model, paving the way 110 for unified treatment of snow and ice. Second, compared with 16-stream albedo solutions, the delta-Eddington adding-doubling approximation provides more accurate albedo estimates, especially under diffuse conditions (Dang et al., 2019). Third, the adding-doubling solution is stable under all conditions except when SZA= 90 • (i.e., when the sun is exactly on the horizon), whereas the tri-diagonal matrix formulation can encounter rare singularities across a broader range of SZA, dependent on the column optical properties and the approximation applied. Fourth, the adding-doubling solver is more computationally efficient 115 than the formulation from Toon et al. (1989). For these reasons, we have transitioned SNICAR to the adding-doubling solver, i.e., SNICAR-AD (Dang et al., 2019). SNICAR-ADv3 combines delta scalings with the Eddington two-stream approximation (Briegleb and Light, 2007), whereas previous versions of SNICAR combined the delta scalings with either the Hemisphericmean or Eddington two-stream approximations as formulated by Toon et al. (1989). Another difference between these solutions is that fluxes under diffuse illumination are solved for via angular integration of direct-beam incidence at eight Gaussian points 120 in SNICAR-AD (e.g., Briegleb and Light, 2007), whereas diffuse incident flux is input as a distinct upper-boundary term in the solution from Toon et al. (1989). The two-stream equations applied in SNICAR-ADv3 are listed in Appendix A. Because of the mathematical singularity that occurs at cos(SZA) = 0, the allowable SZA range is limited to 0 − 89 • .

Spectral grid
The new model simulates albedo across a wavelength (λ) range of 0.2 − 5.0 µm at 10 nm resolution. There are 480 bands 125 in total, centered at: [205 nm, 215 nm, 225 nm, . . . , 4995 nm]. The model was extended from a lower wavelength bound of 300 nm to 200 nm to handle more of the UV spectrum. Though there is insufficient surface irradiance at wavelengths of 200−300 nm to appreciably affect broadband albedo of typical Earth surfaces, surface albedo in this spectral range is important for studies of Mars' polar ice caps (France et al., 2010;Singh and Flanner, 2016;Singh et al., 2018), high-altitude regions of Earth, and photochemical reactions on Earth and other planets (e.g., Singh, 2020). All files in the accompanying library of 130 constituent optical properties for SNICAR-ADv3 are defined on this spectral grid, though many of these definitions require uncertain extrapolation to the shortest and longest wavelengths.

Bulk layer properties
The spectrally-dependent bulk layer properties that enter the two-stream solver depend on spectrally-resolved optical properties and mixing ratios of all constituents present within the layer. The layer extinction optical depth of constituent n depends on its 135 mass extinction cross-section (k e,n , units of m 2 kg −1 ) and layer mass burden (L n , units of kg m −2 ) as: The layer burden of each constituent depends on its mass mixing ratio (q n , units of kg n kg −1 snow ) as: where ∆z is the snow layer thickness and ρ s is the density of snow in the layer. Here, the mass of snow and ρ s include the 140 mass of ice and all other constituents present in the snow. Note that in previous versions of SNICAR, mass mixing ratio was defined as the constituent mass per unit mass of ice. Differences between these two definitions are negligible for the trace mixing ratios of impurities (10 −9 to 10 −4 ) that are typically specified in the model. The definition of q applied here is a true mixing ratio so that q n = 1 and enables users to simulate the albedo of non-ice substrates by specifying q = 1 instead of q = ∞. It also conforms with a definition of snow density consistent with that typically measured in the field, e.g., using a mass 145 measurement that includes all constituents. Ice grains are therefore treated as any other constituent in this representation, but instead of querying for a mass mixing ratio of ice, the mass burden of ice is derived from the user-specified snow density by assuming all non-LAC mass of snow is ice: where ice is identified as constituent n = 1 of N total constituents. 150 The bulk layer optical properties are then formulated using the traditional technique that assumes particles scatter independently in one another's far fields, which is generally valid when inter-particle spacing is substantially longer than the wavelength of interacting light. This assumption is evaluated in detail by Wiscombe and Warren (1980, Section 7), who concluded that although the assumption is not strictly valid for snow, its contribution to bias in snow albedo modeling is likely quite small and with unclear sign, at least when snow density is less than ∼ 450 kg m −3 . With τ n defined for each constituent, including ice, 155 the layer extinction optical depth is the sum from that of all N constituents: The bulk layer single-scatter albedo and scattering asymmetry parameter are calculated with optical depth weighting and scattering optical depth weighting, respectively, of each constituent: The optical properties for ice and other constituents included with SNICAR-ADv3 are described in Sect. 3.

Broadband albedo and surface irradiance
The code and web-based model provide spectral albedo, α(λ), calculated in each of the 480 spectral bands, and solar broadband 165 albedo (α), which is weighted by surface spectral irradiance, F ↓ (λ), as follows: As a practical matter, surface irradiance band fractions are applied and provided in the web-based model output. These fractions are normalized to sum to 1.0 over the SNICAR spectral range such that broadband albedo is simply: where α i is the albedo of band i and f ↓ i is the fraction of surface irradiance within band i.  Table 3. The fine-grain (re = 100 µm) base case is applied here.
Several options for surface spectral irradiance are included in the SNICAR-ADv3 package, derived using different atmospheric profiles listed in Table 1. Surface spectral irradiances are calculated separately from SNICAR with the DISORT-based Shortwave Narrowband (SWNB2) model (Zender et al., 1997;Zender, 1999), using standard atmospheric vertical profiles of water vapor, ozone, and other gases, and a lower-boundary spectral albedo typical of a snowpack with an effective grain size 175 of 100 µm. To obtain f ↓ , we interpolate the output of SWNB2, whose native spectral resolution ranges from 0.3 to 25 nm, to the SNICAR spectral grid. High-altitude options are provided for Summit, Greenland and a generic "high mountain" environment, which are derived by truncating the sub-Arctic summer and mid-latitude profiles to have surface pressures of 796 hPa and 556 hPa, respectively. We also provide a top-of-atmosphere irradiance option, which is useful for modeling surfaces of Mars (e.g., Singh and Flanner, 2016) and other bodies in our solar system with thin atmospheres, and for exploring impacts 180 of atmospheric attenuation on broadband albedo. The top-of-atmosphere spectral irradiances used to drive SWNB2 are from Matthes et al. (2017), averaged over three solar cycles. The percentages of total solar irradiance in this dataset residing outside the SNICAR spectral range, at λ < 0.2 µm and λ > 5.0 µm, are only, respectively, 0.013% and 0.066%. Surface spectral irradiances and α are sensitive to cloud and atmospheric conditions (e.g., Wiscombe and Warren, 1980), and can also be somewhat sensitive to the radiative transfer model and top-of-atmosphere irradiance data used (Bair et al., 2019).

185
Spectral irradiances associated with clear-sky or cloudy atmospheres are selected internally in the model to match the userspecification of direct or diffuse incident light, respectively. The cloudy irradiances are modeled with a liquid cloud of optical thickness 10 at λ = 500 nm, located at a pressure of 800 hPa or in the bottom-most atmospheric layer of profiles with surface pressure less than 800 hPa. The cloud droplet properties vary spectrally based on Mie calculations of water spheres with effective radius of 10 µm. Clear-sky irradiances are calculated for the full range of SZA (0 − 89 degrees at 1 degree resolution)  Ice grain optical properties are derived from Mie calculations of ice spheres and adjustments for scattering by non-spherical particles (Fu, 2007;He et al., 2017). Users specify one of three datasets of ice refractive indices, leading to retrieval of distinct Mie properties from the accompanying library of optical properties in netCDF format. The complex refractive indices (M = m r + m i i) from the original data are linearly interpolated to the SNICAR spectral grid, except as described below. The H 2 O ice refractive index options in SNICAR-ADv3 are: 200 1. Warren (1984) / Perovich andGovoni (1991): m r data from Warren (1984) are applied across the solar spectrum. m i data from Perovich and Govoni (1991) are applied over λ = 250 − 400 nm and data from Warren (1984) are applied at all other wavelengths. The data points from Warren (1984) at 210 and 250 nm are adjusted slightly to achieve a smooth transition in interpolated data between the two datasets over λ = 210 − 250 nm. 2. Warren and Brandt (2008): m r and m i data compiled by Warren and Brandt (2008) are applied across the solar spectrum.

205
Note that the ice absorption data over λ = 600 − 1400 nm in this compilation originate from Grenfell and Perovich (1981).
3. Picard et al. (2016) / Warren andBrandt (2008): m r data from Warren and Brandt (2008) are applied across the solar spectrum. The 320-600 nm m i data from Picard et al. (2016) are adopted, and data from Warren and Brandt (2008) are applied at longer wavelengths. The 320 nm value reported by Picard et al. (2016) is also applied constantly over 210 200-320 nm, thus retaining the same feature of constancy over this spectral range that is exhibited in the Warren and Brandt (2008) dataset.
Values of m i , which govern ice absorptivity, differ by two orders of magnitude or more in the short blue and UV part of the spectrum (Fig. 2a) across these three datasets, but are very similar in the NIR (defined here as λ = 0.7 − 5.0 µm), where ice is absorptive. UV and blue ice absorptivity is extraordinarily low and notoriously challenging to measure because very 215 long optical path lengths are needed. Dataset #1 is largely superseded by Dataset #2, the more recent and comprehensive compilation from Warren and Brandt (2008), but is retained because the original version of SNICAR applied these data, and they arguably still represent an upper bound of ice absorptivity. Warren and Brandt (2008) argue that pure ice m i is no larger There is much less uncertainty in m r , and agreement is excellent between Warren (1984) and Warren and Brandt (2008) in the solar spectrum (Fig. 2c). (Note that these and subsequent spectral figures use log-spacing to emphasize the shorter wavelengths which contribute disproportionately to broadband albedo.) Mie properties are calculated using the solver from Bohren and Huffman (1983) for a large range of ice grain effective 225 radii (r e , or surface-area weighted mean radius), and using each of the three sets of H 2 O ice refractive indices. The optics library includes individual files for r e values ranging from 30-1500 µm, discretized by 1 µm. Properties are first calculated for 5000 individual spheres with log-spaced radii ranging from 0.05 to 10000 µm, then weighted according to a log-normal size distribution for each of the specified r e values. The relative number of particles with radius r, n(r), in the log-normal distribution is: where σ g is the geometric standard deviation of the size distribution, fixed at 1.5 for all of our ice particle calculations, and r n is the number-median radius of the size distribution, related analytically to r e as: The resolved effective radius of the discretized distribution is: 235 r e = r 3 n(r) r 2 n(r) Comparisons of the user-specified and resolved effective radii for our calculations show near-perfect agreement, indicating that the size distributions are adequately resolved with these model parameters. We assume a bulk H 2 O ice density of 917 kg m −3 for deriving mass-normalized optical properties.
Finally, users can specify one of four ice particle shapes, identified by He et al. (2017)  Optical properties used for spheres are simply the Mie properties of the user-specified r e . For the other shapes, the userspecified r e represents the radius of the sphere with equivalent specific surface area (SSA, or surface-area to volume ratio) as that of the non-spherical particle (He et al., 2017). Because equal-SSA spheres provide good proxies for the extinction cross-section and single-scatter albedo of non-spherical particles (e.g., Grenfell et al., 2005), Mie-generated k e and ω values are also retrieved for the non-spherical particles. The scattering asymmetry parameter (g), however, is adjusted using the 250 parameterization described by He et al. (2017) and He et al. (2018), structured on the scheme developed by Fu (2007).
Corrections to g are parameterized as a function of wavelength, user-specified effective radius, particle shape, and particle aspect ratio (the ratio of grain width to length). The asymmetry parameter for hexagonal plates (g hex ) follows Fu (2007): where g is parameterized as a function of aspect ratio using the wavelength-dependent coefficients listed in Tables 1 and 2 of   255 Fu (2007). For other shapes, g is derived from g hex with: where C is a correction factor that depends on r s , the sphere radius with the same orientation-averaged projected area-tovolume ratio of the non-spherical particle, and on the particle shape factor (SF ), defined as the ratio of r s of the nonspherical grain to that of an equal-volume sphere: where a 0−2 are wavelength-dependent coefficients listed in Table 3 of He et al. (2017). Note that for convex shapes (here, spheroids and hexagonal plates) r s = r e and for concave shapes r s > r e in general, and for Koch snowflakes in particular r s = r e /0.544 (He et al., 2018).
These empirical parameterizations are based on detailed geometric optics ray-tracing calculations (Fu, 2007) that include 265 surface wave interactions and diffraction (Liou et al., 2011(Liou et al., , 2014Liou and Yang, 2016). The parameterized values of g at 8 wavelengths are then interpolated to the SNICAR spectral grid using a piecewise cubic Hermite interpolating polynomial (Matlab function 'pchip'). Within the code, users can alter the aspect ratio (within the range of 0.1 − 20) and shape factor, but the web-based model adopts the default aspect ratios listed above and associated shape factors, which are presented in Table 1 of He et al. (2017). The parameterized value of g can exceed 1.0 in rare instances with large spheroids, and consequently g is 270 capped at 0.99. Implications of particle shape for simulated albedo are described in Sect. 4.1.

CO 2 ice
To facilitate the simulation of carbon dioxide ice surfaces like those found on Mars, we also include optical properties for CO 2 ice grains, as applied in SNICAR by Singh and Flanner (2016) and Singh et al. (2018). We apply concatenated CO 2 ice refractive indices over wavelengths of 0.2-1.8 µm and 1.8-5.0 µm from Hansen (2005) and Hansen (1997), respectively. UV spectrum and m i values are somewhat uncertain. CO 2 ice m i is relatively featureless at wavelengths shorter than 1 µm and reaches a minimum near 250 nm (Hansen, 2005). Unlike H 2 O, however, CO 2 ice m i is punctuated by numerous sharp absorption bands throughout the NIR, though averaged over the NIR CO 2 ice is less absorptive than H 2 O ice.

280
Optics files are included in the library for CO 2 r e ranging from 5 − 1500 µm at 1 µm resolution, calculated using the same approach and parameters for log-normal size distributions of H 2 O. The library also includes optics files for r e values up to 20,000 µm at 500 µm resolution. This larger range is provided to accommodate sensitivity studies that probe uncertainty in CO 2 snow morphology on Mars (e.g., Singh and Flanner, 2016). We assume a bulk CO 2 ice density of 1500 kg m −3 . The same parametric correction of g for non-spherical ice particles applied to H 2 O grains (He et al., 2017) is enabled for CO 2 ice. The 285 adjustment is specific to H 2 O and hence will be biased somewhat when applied to CO 2 , but perhaps not by much because of their similarity in m r (Fig. 2b), which for a given shape is the main determinant of the scattering phase function, especially in weakly absorbing portions of the spectrum (e.g., Räisänen et al., 2015). We also note that some of the available shape options, notably hexagonal habits, are unrealistic for CO 2 snow. Hence, the non-spherical parameterizations for CO 2 ice should be used with caution.

Light-absorbing constituents
The primary application of SNICAR has been to study the effects of LAC on snow albedo. In this section we describe the optical properties of various types of LAC that are included in the accompanying library of Mie properties. The variation in observed LAC properties is immense and far greater than the subset of properties included here. The library presented here constitutes a sample of properties that have been used in previous works and which are broadly representative of constituent 295 properties. We remind readers, though, that in many cases properties specific to one's application or measurements may need to be generated and applied.
A key property that governs the potency of LAC in perturbing snow albedo is the mass-absorption cross-section (k a ): Although the LAC scattering properties are also needed in our representation (Eqs. 4-6), they generally have negligible 300 bearing on snow albedo at typical LAC mixing ratios (∼ 10 −10 to 10 −4 ), as snow scattering at λ < 1.0 µm is dominated by ice grains, and impurity scattering is too small to increase snow reflectance in the absorptive NIR bands of ice. Hence we only report k a values in this section, and remind readers that scattering properties are also included in the optics files. All of the properties described below are generated with Mie calculations using the standard and coated-sphere solvers of Bohren and Huffman (1983). Table 2 lists the resolved size distribution parameters and 500 nm k a values for all LAC species included in 305 the library. In addition to the resolved effective (r e ) and number-median (r n ) radii, we also list the resolved mass mean (or volume mean) radius, r m : In the log-normal size distribution, r m is related analytically to r n as:

Black carbon
Black carbon (BC) is generally defined as the strongly-absorbing component of carbonaceous aerosols, and consists largely of elemental carbon. We apply the BC optical properties described by Flanner et al. (2012). These are derived from the spectrallyresolved refractive index parameterization provided by Chang and Charalampopoulos (1990), but adjusted with linear offsets to achieve m r = 1.95 and m i = 0.79 at λ = 550 nm, as recommended by Bond and Bergstrom (2006) in their comprehensive 315 review. The parameterization we apply over 0.2 − 5.0 µm, with λ in units of micrometers, is therefore: a The size parameters and densities of coated particles refer to the whole particle, whereas ka is defined with respect to core mass only. Core particle sizes are identical to their uncoated counterparts. b The ranges of dust densities and ka values span the four types of dust described in Sect. 3.3.3 We assume log-normal size distributions with r n = 40 nm and σ g = 1.8, intermediate values from a compilation of measurements of fresh soot (Bond et al., 2006b). The particle density is assumed to be 1270 kg m −3 , which achieves k a = 7.5 m 2 g −1 320 at 550 nm, conforming with the central recommendation of Bond and Bergstrom (2006). BC in snowpack has been found to have a larger size distribution than atmospheric BC in some studies (Schwarz et al., 2013), but not in others (Sinha et al., 2018). The two types of BC that can be specified are uncoated (representative of fresh, externally-mixed soot) and sulfatecoated (representative of aged, internally-mixed soot). The refractive sulfate coating, with properties from Hess et al. (1998), has an outer radius 2.15× larger than the uncoated BC and produces an absorption enhancement, per unit mass of BC, of 1.5 325 (Bond et al., 2006a). Various new and more sophisticated treatments of internally-mixed BC in ice have been explored (e.g., Liou et al., 2011;Flanner et al., 2012;He et al., 2017). All of these studies indicate that BC absorbs more solar energy, per unit mass, when it is embedded inside weakly absorbing matrices like ice, with absorption enhancement factors typically ranging from 1 − 2, but with substantial variability associated with the morphology and size distributions of the inclusions, coatings, and matrices (e.g., Cappa et al., 2019). The sulfate-coated BC species therefore serves as a first-order proxy for BC that is 330 internally-mixed in any weakly-absorbing, refractive agent (ice, sulfate, organics, etc). More accurate techniques, such as the representation of multiple inclusions of BC residing in non-spherical ice particles (He et al., 2017), should be applied when The "Base Case" algae consists of dry cell mass fractions of 1.5% chlorophyll-a, 0.5% chlorophyll-b, 5% photoprotective carotenoids, and 0% photoreceptive carotenoids. Single-pigment scenarios with 1% mass fraction are shown in the other curves. The mean cell radius in all cases is 10 µm. Note that ka is normalized to the mass of the entire cell, which is 78% water by volume. The large ka values for snow algae at λ > 2.5 µm (exhibited in all curves) are due to absorption by water.
particle/inclusion habits and mixing geometry are known, but it is quite rare that these microphysical details are measured or simulated. The spectral profiles of k a for the two species of BC included in the library are shown in Fig. 3a.

Brown carbon 335
So-called "organic carbon" is co-emitted with BC in the combustion process, and generally consists of weakly-absorbing aerosol. Although most organic carbon has a sufficiently high ω to not appreciably influence snow albedo, the "brown carbon" subset absorbs strongly in the UV/blue and weakly in the red (Andreae and Gelencsér, 2006), and is therefore of interest for snow albedo impacts (e.g., Lin et al., 2014). The optical properties of organic and brown carbon exhibit enormous diversity due to differences in solubility, size distribution, and molecular composition that are associated with combustion characteristics 340 and the composition of the parent material (e.g., Sun et al., 2007;Laskin et al., 2015). Organic and brown carbon aerosols also originate from primary emission of biogenic species and secondary aerosol formation from biogenic gas emissions, adding even more to the diversity of the optical properties of this class of aerosol.
Here, we provide properties for two versions of brown carbon with identical size distributions and sulfate coatings as the BC species, but derived from brown carbon imaginary refractive indices measured by Kirchstetter et al. (2004). These properties 345 are more absorptive than most other derivations of organic and brown carbon optical properties (e.g., Hoffer et al., 2006;Sun et al., 2007;Chen and Bond, 2010;Lack et al., 2012;Laskin et al., 2015), and thus may serve as an approximate upper bound on albedo impacts from organic/brown carbon. The native data from Kirchstetter et al. (2004) extend from 350 − 700 nm, and are extrapolated to 200 nm using a cubic polynomial (Matlab function 'pchip') that preserves the tendency towards increasing absorptivity with shorter wavelengths. Data at λ > 700 nm are linearly tapered down to a value of 10 −5 at λ = 5 µm.

350
The real component of the refractive index is held constant at 1.53. A chart showing how these properties compare with other measurements of brown carbon is shown in Figure 12 of Laskin et al. (2015). The k a values from our determination are presented in Fig. 3b, which shows UV absorptivity that is nearly as large as that of BC, but a much sharper decline in absorptivity with wavelength. As with BC, the absorption enhancement associated with the sulfate coating is ∼ 1.5 in the visible, but it could actually range from 1 (i.e., no enhancement) to 2 or more.

Mineral dust
Snow darkening resulting from the deposition of aeolean dust has been documented in many cases (e.g., Painter et al., 2007).
Per unit mass, dust is much less absorptive than BC, but is often present in mixing ratios that are orders of magnitude greater.
Dust particles in the atmosphere and snow are typically larger, and span a wider range of sizes, than those of BC, motivating the creation of more options in SNICAR for dust particle sizes. Soil dust also exhibits substantial diversity in mineral composition, 360 leading to variability in absorptivity. Minerals with oxidized iron, in particular, tend to strongly absorb light and thus iron content is one of the key determinants of dust absorptivity.
In SNICAR-ADv3 users can specify dust mixing ratios in each of the following five particle size bins: 1. 0.05 < r < 0.5 µm 2. 0.5 < r < 1.25 µm 365 3. 1.25 < r < 2.5 µm 4. 2.5 < r < 5.0 µm 5. 5.0 < r < 50 µm The smallest four of these bins match those that have been used in global aerosol transport studies (Zender et al., 2003;Mahowald et al., 2006;Scanza et al., 2015). The largest bin was added because measurements show that very large dust 370 particles are often present in snowpack (e.g., , especially in patchy snow of arid environments. The dust optical properties for each of these bins are derived from Mie calculations assuming a log-normal size distribution with analytic r e = 1.38 µm and σ g = 2.0 (and a mass median diameter of 3.5 µm), matching the parameters used by Scanza et al. (2015) and others, and truncated to each of the size bounds listed above. We perform Mie calculations on 1000 log-spaced particle radii within each size bin. The resolved r e , r n , and r m values for each size bin are listed in Table 2. 375 We provide optical properties for four compositions of dust that have been used in previous dust-on-snow studies, representative of properties from different geographies and derived using diverse techniques and data: 1. Saharan dust: For this species we apply mineral fractions from the "central hematite" scenario presented by Balkanski et al. (2007), which produced good agreement with AERONET measurements in regions strongly affected by Saharan dust. These mineral volume fractions are: 31.5% illite, 24% kaolinite, 23% montmorillonite, 14% quartz, 6% calcite, and 380 1.5% hematite. We apply the Bruggeman mixing approximation (Bruggeman, 1935) to derive a particle-mean dielectric constant ( = M 2 ) from that of each of the i minerals ( i ) and mineral volume fractions (V i ). Specifically, we use an iterative secant method  to find such that: We then use the particle-average refractive indices in Mie calculations with the size distribution parameters listed above.

385
The spectrally-resolved mineral refractive indices we use are essentially those tabulated in the supplement of Scanza et al. (2015), which originate from several sources (Egan et al., 1975;Egan and Hilgeman, 1979;Querry, 1987;Long et al., 1993;Rothman et al., 1998, A.H.M.J. Triaud, personal communication, 2005. We also assume the same mineral from Mie inversions and measurements of the spectral albedo of optically-thick dust samples of known particle size distribution . This technique finds the m i value at each wavelength that optimizes the comparison between measured and modeled dust spectral albedo. The real component (m r ) was assumed to be 1.525 for the inversions. Due to an apparent inconsistency between the m i values reported by  and their measured dust 395 particle size distributions, we we-ran the inversion to produce revised m i values, assuming a lognormal size distribution with the same resolved effective radius (r e = 2.3 µm) of the measurements. The revised m i values are larger than those from the original inversion. We extrapolated m i to wavelengths shorter than 350 nm with Matlab function 'pchip'. At λ > 2.5 µm, m i is assumed to equal 0.00165, the 2.5 µm value derived from our inversion. We assume a dust density of 2600 kg m −3 and σ g = 1.32, consistent with . This dataset has the advantage of deriving from optical refractive indices of individual minerals, as occurs in the modeling of dielectric properties for the other types of dust we provide. This technique may, however, be susceptible to errors in albedo measurements, unaccounted light transmission, the assumption of spherical particles, and uncertainty in m r .
3. Greenland dust: Our approach for calculating Greenland dust properties follows that used for the Saharan dust, but with 405 mineral composition from Polashenski et al. (2015). These mineral fractions were estimated from linear mixing model calculations using measurements of aluminum, iron, non-sea salt calcium, water insoluble potassium, and arsenic in snow samples collected across Greenland. Three scenarios of hematite fraction were derived by Polashenski et al. (2015),  (Singh and Flanner, 2016;Singh et al., 2018). The refractive indices used to generate our Martian dust are derived from spectral measurements from the Mars Reconnaissance Orbiter (Wolff et al., 2009. We apply the shortest-wavelength (λ = 263 nm) value from this dataset constantly across the short-spectrum (λ = 200 − 263 nm) portion of SNICAR. We also assume a particle density of 2000 kg m −3 in deriving the optical properties, consistent with Singh and Flanner (2016).

420
Spectral mass-absorption cross-sections (k a ) for these dust types and size bins are shown in Fig. 3c-f. Per unit mass, small dust particles are more effective absorbers than larger ones. All dust species are more absorptive in the blue than red, and smaller particles exhibit a stronger spectral dependency in absorption. The jagged feature near 350 − 400 nm in the Saharan and Greenland dust originates from the hematite refractive indices. Mie non-linearities are also apparent. For example, the finest Greenland dust is the most absorptive of all types, which is not surprising given its relatively large hematite content, 425 but the mid-size Martian dust is more absorptive than its counterparts. Overall, Colorado dust has lower visible m i than the other dust species included in our library, but is an important component of our collection because of the diversity in dust absorptivity that it indicates, the importance of southwest Colorado for dust-on-snow studies (e.g., Skiles and Painter, 2019), and the unique way in which the properties were determined. Singh and Flanner (2016) Fig. 3f for reference). These properties are also included in the accompanying library.

Volcanic ash
Deposition of volcanic ash on snow, while highly episodic, can substantially lower albedo (e.g., Conway et al., 1996;Young et al., 2014;Gelman Constantin et al., 2020). We apply ash refractive indices from the "central forcing" scenario presented by Flanner et al. (2014) (their Figure 2). These are derived from a variety of measurements following the 2010 Eyjafjallajökull eruption, including chemical and optical measurements of ash in the atmosphere and snow Bukowiecki et al., 2011), and inversions of sun photometer, ground-based lidar, and particle soot absorption photometer measurements (Hervo et al., 2012;Toledano et al., 2012;Weinzierl et al., 2012). These central scenario m i values are very similar to those of ash collected from the 1980 Mount St. Helens eruption (Patterson, 1981). The 500 nm m i value we apply is 0.0044. T -Matrix calculations (Mishchenko and Travis, 1998) and found that, although the scattering properties differ substantially between equal-volume spheres and non-sphere ash particles, the k a values differed by at most 16%, leading to the conclusion that particle shape contributed only second-order uncertainty to the estimation of k a , the key optical property for snow impurity albedo studies. The volcanic ash k a values are depicted in Fig. 3g, which shows lower absorptivity than all types of dust presented except for Colorado dust.

Snow algae
Algae blooms occur in seasonal snowpack and in the ablation zones of ice sheets and glaciers (e.g., Takeuchi et al., 2001;Painter et al., 2001;Cook et al., 2017). Snow and glacier algae exhibit enormous genetic and phenotypic diversity (e.g., Lutz et al., 2014), and are thus difficult to represent using traditional physical models. Our goal here is to provide a firstorder representation that accounts for light-absorption by some of the key pigments present in algae, utilizing the approach 455 introduced by Cook et al. (2017). Given the focus of this study on snow albedo modeling and the pending development of a separate glacier albedo model (Whicker et al., 2021), we focus exclusively on snow algae, which are often reddish or greenish and can be detected remotely because of their unique spectral characteristics (e.g., Painter et al., 2001;Wang et al., 2018). Our approach should be considered an experimental work-in-progress that will evolve with improved observational constraints on the nature and determinants of snow algal absorption, as well as improved modeling techniques for representing heterogeneous 460 and irregularly-shaped media like algal cells (Cook et al., 2020).
We utilize the spectral mass absorption cross-sections of various algal light-absorbing pigments from Dauchet et al. (2015) to derive pigment m i values and assume a constant m r = 1.5 for all pigments (Pottier et al., 2005;Cook et al., 2017). We apply the Bruggeman approximation described earlier to mix these properties (and those of water) in various proportions and then conduct Mie calculations using the cell-average refractive indices. We assume a cell water volume fraction of 0.78 0. The numbers in brackets indicate the discretized dry cell mass fractions for which we created optical properties. Chlorophylla pigments are present in all photosynthetic snow algae, and the accessory chlorophyll-b pigment is also common (e.g., Remias et al., 2005). The carotenoid properties refer to generic classes of pigments, with the photoprotective carotenoids particularly relevant, as algae often produce these pigments to protect themselves from UV radiation and excessive solar heating (e.g., Remias et al., 2005). The dry cellular matter unaccounted for with the pigments is assumed to have constant M = 1.5 + 0i 475 (i.e., to be non-absorbing) (Pottier et al., 2005;Cook et al., 2017). All dry matter is assumed to have a density of 1400 kg m −3 (Dauchet et al., 2015), yielding a mean cell density of ρ alg = 1088 kg m −3 . The m i properties from Dauchet et al. (2015) are only defined down to λ = 350 nm, and the use of extrapolating functions is suspect due to the sharp absorption features of pigments. We therefore take a somewhat conservative approach of assuming that m i at λ = 200 nm (our lower model boundary) is equal to half of the published values at λ = 350 nm, and interpolate linearly between 200 nm and 350 nm.

480
We assume Gaussian size distributions of algal cells, instead of log-normal distributions as applied for all other LAC. The standard deviation (σ) of the distribution is assumed to be 10% of the mean cell radius (r alg ), and for each mean radius we perform Mie calculations over 200 equally-spaced size bins ranging across r alg ± 4σ. The mean cell radii for which we generated properties are: 1, 2, 5, 10, 15, 20, 25, 30, 40, and 50 µm, with r alg = 10 µm set as the default. Remias et al. (2005) report a mean diameter of 14.9 ± 5.7 for quasi-spherical cells of C. nivalis, one of the most prevalent species of algae in 485 mountain and polar snowpack.
We generate a library of algae optical properties, building on the "BioSNICAR" library developed by Cook et al. (2017), with five dimensions: mean cell radius and dry cell mass fractions of each of the four pigments listed above. In total, the sampling and discretization over these five dimensions results in 125,440 unique combinations of Mie properties, and the algae component of the optics library therefore occupies the majority (8.2 GB) of its disk space. Users specify an algae cell number 490 concentration (N ) in units of cells per mL of meltwater, which is consistent with how the microbiology community measures and reports cell concentrations. Internally, however, SNICAR treats all LAC in terms of mass mixing ratios and mass-specific optical properties (e.g., Eqs. 1-3). We therefore convert the user-specified cell number concentration to a mass mixing ratio (q alg ) using Gaussian statistics of the size distributions applied in Mie calculations, i.e.: q alg = 10 −15 N ρ alg 4 3 π r 3 alg + 3r alg σ 2 (21) 495 where the factor of 10 −15 is needed with r in µm, N in # mL −1 , and ρ in kg m −3 .  Dry cell mass fraction of photosynthetic carotenoids 0% * Snowpack density and underlying ground albedo have no impact on simulated albedo when the snowpack is semi-infinite.
The base case cell radius, and that of all curves shown in Fig. 3h, is 10 µm, slightly larger than that measured by Remias et al. (2005) for C. nivalis, and slightly smaller than the medium scenario given by Cook et al. (2017). Sharp absorption features are apparent in Fig. 3h, including the chlorophyll-a absorption peaks near 430 and 670 nm that are exploited for remote sensing retrievals (Wang et al., 2018), and the more closely spaced chlorophyll-b peaks. These become partially obscured when other 505 pigments are present, as seen in the base case absorption profile, and can become even more obscured when other broad spectrum photoprotective pigments are present (Cook et al., 2020). Note that k a shown in Fig. 3h is that of the entire cell, which is 78% water by volume. The absorption peak near λ = 3 µm seen in all cases is due to cell water. Normalization of cell k a to just the dry cell mass, as is sometimes reported, would produce 3.53× larger values.
The techniques we utilize to represent algae are imperfect. First, the Bruggeman approximation assumes perfectly homo-510 geneous mixing, whereas pigments have discrete sizes, shapes, and positions within the cell. In cases where such features of the cell are known, techniques like the dynamic effective medium approximation (Chýlek and Srivastava, 1983), discrete dipole approximation (Draine and Flatau, 1994), or geometric-optics surface-wave approach (Liou et al., 2011) might be more appropriate. Second, Mie approximations assume spheres. While glacier algae are often filamentous and generally highly nonspherical (e.g. Yallop et al., 2012) and more appropriately represented with geometric optics or other techniques (Cook et al.,515 2020), some of the most common types of snow algae like C. nivalis actually are quasi-spherical (Remias et al., 2005;Lutz et al., 2014), suggesting that Mie approximations may be reasonable for snow algae. Third, other types of pigments are present in algae, particularly specific types of carotenoids whose spectral properties may differ from those applied here (e.g., Remias et al., 2012;Williamson et al., 2020). Our technique of generating an offline library of optical properties that comprehensively spans the parameter space precludes the use of many types of pigments. A more efficient technique may be to generate algae 520 optical properties concurrently with queries for albedo calculations, using Mie Theory or other parametric techniques. We leave it to future studies, however, for deeper dives into the representation of snow algae.

Model Sensitivities
In this section we briefly explore the sensitivity of modeled spectral and broadband albedo to some of the model parameters and design features described above. Due to the immensity of the parameter space, this analysis is not comprehensive. We define 525 a base state of model features, listed in Table 3, from which we perturb individual parameters. Because effective snow grain size (r e ) is a dominant source of variability in snow albedo across Earth's snowpacks, we explore many of the sensitivities under both a coarse-grain (aged snow) base case with r e = 1000 µm and the default fine-grain (fresh snow) base case with r e = 100 µm. The base parameters listed in Table 3 are also the default parameters used in the web-based model.

Pure snow 530
Surface spectral irradiance has no impact on spectral albedo, but influences broadband albedo (α) as shown in Eq. 7. Table 1 lists the base case α for all atmospheric profiles. Under cloudy skies a higher proportion of the surface irradiance lies in the visible spectrum, where snow is most reflective, and hence α is larger under cloudy conditions (e.g., Wiscombe and Warren, 1980). Other factors equal, clear-sky α also increases with increasing water vapor content due to near-infrared absorption in the atmosphere, as seen by comparing α under winter and summer atmospheric profiles (Table 1). 535 Fig. 4 depicts modeled spectral and broadband albedo sensitivities to ice refractive index dataset (panels a, b, and d), snow r e (panels c and d), SZA with fine and coarse grains (panels e and f), and grain shape with fine and coarse grains (panels g and h). Sensitivity to the choice of H 2 O ice refractive index is negligible in most of the spectrum. Even in the UV/blue spectrum α differs by < 0.01 with r e = 100 µm between the two most recent and most reliable datasets (Warren and Brandt, 2008;Picard et al., 2016), but these differences increase to 0.025 with r e = 1000 µm. Differences are even larger between the Warren (1984) 540 and Warren and Brandt (2008) datasets but the former is obsolete. Ice absorption is so weak in the Warren and Brandt (2008) dataset that UV/blue snow albedo is nearly 1.0 regardless of grain size (Figs. 4a and 4b). CO 2 snow is more reflective than H 2 O snow, except in the highly absorptive NIR bands of CO 2 (Fig. 4d, Singh and Flanner, 2016). Snow r e (or SSA) is well-known to be one of the primary controls on snow albedo (e.g., Warren and Wiscombe, 1980), as shown in Fig. 4c. The influence is index datasets (a,b,d), ice particle effective grain size (c,d), solar zenith angle (e,f), and ice particle shape (g,h). Except for the parameter being varied, all snowpack properties are those of the fine-grain (re = 100 µm) or coarse-grain (re = 1000 µm) base case defined in Table 3.
Solar broadband albedos for each curve are also printed on the figures, with order corresponding to that of the legend.
primarily in the NIR part of the spectrum, and α ranges from 0.848 − 0.747 for r e = 100 − 1000 µm in our base case. It is also 545 well-known that clear-sky albedo increases with increasing SZA, as manifested in Figs. 4e and 4f. The range of α across SZA of 0 − 85 • is 0.819 − 0.862 in the fine-grain base case (Fig. 4e), and is 0.707 − 0.758 in the coarse-grain base case (Fig. 4f).
At large SZA, however, the albedo dependency is weakened due to a red-shift in surface spectral irradiance that occurs with Rayleigh scattering (van Dalum et al., 2020), which by itself decreases α and thus counteracts the tendency towards larger α due merely to increasing incidence angle. We see noticeable divergence in α between simulations with and without SZA-550 dependent surface irradiance at SZA >∼ 65 • due to more complex spectral shifts associated with long path-length gaseous absorption and Rayleigh scattering.
Snow grain shape and morphology have gained renewed attention for their importance in snow radiative transfer (e.g., Kokhanovsky and Zege, 2004;Liou et al., 2011;Libois et al., 2013;Räisänen et al., 2015;Dang et al., 2016;He et al., 2017). Figs. 4g and 4h show that non-spherical ice grains are associated with higher albedo than equal-SSA spherical grains, 555 particularly in the NIR. This is due to differences in the scattering asymmetry parameter (g), which is actually the only optical parameter that changes with grain shape in our parameterization (He et al., 2017). Non-spherical ice grains scatter less strongly in the forward direction than spheres (e.g., Fu, 2007;Libois et al., 2013;Räisänen et al., 2015;Dang et al., 2016), which decreases the penetration depth of radiation in snow and increases the probability that photons will re-direct out of the top of the snowpack, thus increasing albedo when all other factors are equal. For the three non-spherical grain shapes and default 560 aspect ratios considered in SNICAR-ADv3, Koch snowflakes have the smallest g and largest α relative to spheres, followed by hexagonal plates and spheroids. The range in α across all four shapes is 0.827 − 0.854 in the fine-grain case (Fig. 4g) and is 0.712 − 0.756 in the coarse-grain case (Fig. 7h). Thus, as with most other factors, the influence of grain shape is greater in larger-grain snow. Because it is now clear that spherical ice grains have unrealistically large values of g, we assign a nonspherical grain shape, hexagonal plates, as the default grain shape in SNICAR-ADv3, though with the recognition that grain 565 shape is enormously diverse in real snowpack. Fig. 5 depicts the dependency of albedo on snow thickness, snow grain size, and grain shape. Light penetrates most deeply in coarse-grain, spherical, low-density snow, translating into a greater snow thickness needed to completely mask the underlying ground, or to achieve "semi-infinite" thickness. In the fine-grain base case (Fig. 5a), light penetration is shallow and α is high (0.762) even with a 1 cm thick snowpack overlying a relatively dark (α = 0.25) surface. In the coarse-grain base case ( Fig.   570 5b), however, we see stronger dependency of α on snowpack thickness. Even with a 20 cm thick snowpack, α is 0.026 less than its semi-infinite value. Finally, we see the importance of grain shape once again in Fig. 5c, which shows α ranging from 0.610 − 0.705 for different shapes with r e = 1000 µm and a snowpack thickness of 10 cm. These ranges narrow with smaller r e and thicker snowpack.

575
Here, we briefly describe some of the modeled albedo sensitivities to LAC in SNICAR-ADv3. Sensitivities to BC and the variant of brown carbon we apply are shown in Fig. 6. Part-per-billion mixing ratios of both species can substantially reduce snow albedo, and more so in coarse-grain than fine-grain snow. In the base case, the reduction in α caused by 10 ppb of uncoated BC increases from 0.0023 to 0.0070 with an increase in r e from 100 µm to 1000 µm. The ratio of coarse-grain α reduction to fine-grain α reduction, greater than 3 at small q BC , decreases slightly to 2.7 with q BC = 1000 ppb. The albedo spectra show 580 that BC impacts snow albedo rather uniformly across the visible spectrum, and with some impact up to λ ≈ 1.0 µm, whereas brown carbon exerts a much greater influence at UV/blue wavelengths and has almost no impact at λ > 0.7 µm. Broadband albedo reduction from brown carbon is 0.24 − 0.30× that from equal mixing ratios of BC over the parameter ranges shown in Fig. 6, though we again remind readers that the k a values we apply for brown carbon are higher than most estimates.
The sulfate-coated versions of BC and brown carbon, with larger k a , drive larger albedo reduction per unit mass of particle.

585
Across the r e and mixing ratio ranges shown in Fig. 6, the coated particles reduce α by 1.2 − 1.4× more than their uncoated counterparts (not shown). We suggest that these coated species may serve as reasonable proxies for internally-mixed particles in ice, as the albedo reduction enhancement was shown to be similar (1.3−1.6×) for BC residing inside snow grains compared with externally-mixed BC (He et al., 2018).
An important consequence of shallower light penetration in snow with non-spherical ice grains is that the albedo perturbation 590 from LAC is reduced (Dang et al., 2016), as fewer LAC are exposed to light. Figs. 6c and 6f show the broadband albedo reduction from BC and brown carbon in snow with different ice grain shapes. The use of hexagonal plates as the default grain shape in SNICAR-ADv3 instead of spheres reduces the simulated albedo impact from BC by ∼ 24%, with only slight variability associated with grain size or BC amount. The BC albedo perturbation in snow composed of spheroids and Koch snowflakes is ∼ 9% and ∼ 31% less, respectively, than in snow composed of spheres, and these ratios are very similar with 595 brown carbon (Fig. 6f). Because this variability is associated with the penetration depth of light in snow, it also highlights the importance of accurately representing the vertical distribution of LAC, which can be modeled with multi-layer configurations of SNICAR-ADv3.
Albedo changes induced by dust and volcanic ash are shown in Fig. 7. Mixing ratios up to 1000 ppm of fine (0.05 < r < 0.5 µm) Saharan dust reduce albedo by up to 0.15 in fine-grain snow (Fig. 7a) and 0.21 in coarse-grain snow (Fig. 7b). This type mixing ratio of 100 ppm, we see that the albedo reduction from Saharan dust ranges from 0.028 − 0.053 in the fine-grain base case, depending on the size bin of particles, with fine dust having the greatest mass-normalized influence (Fig. 7c). For a given set of snow conditions and particle size, dust type contributes to modest variability in albedo change (Fig. 7d), with Martian dust the most absorptive and Colorado dust the least. Other observed mixtures of dust would produce greater variability (Scanza 605 et al., 2015), however, and our representation of dust may evolve in the future to accommodate this diversity. The volcanic ash from Eyjafjallajökull has a similar impact on α as Saharan dust, though via a more uniform spectral influence that gives the snow more of a greyish/brownish appearance (Fig. 7e). As with dust, larger volcanic ash particles reduce snow albedo less per unit mass than smaller ones (Fig. 7f), and the disparity widens with increasing snow grain size (not shown).
Snow algae can introduce unique spectral characteristics to snow albedo (Fig. 8). With the base case snow and algae prop-610 erties (Table 3), algae number concentrations of 10 3 − 10 6 mL −1 reduce snow albedo by 0.004 − 0.184 (Fig. 8a), increasing to 0.014 − 0.336 in the coarse grain base case, which for algae is a much more realistic scenario because they generally only proliferate in snow that has been wet for an extended period of time. Observed cell concentrations in snow and glacier algal blooms are ∼ 10 3 − 10 5 mL −1 (Painter et al., 2001;Yallop et al., 2012;Lutz et al., 2014;Cook et al., 2020;Williamson et al., 2020), though measurements are often made from a thicker layer of snow than that in which most of the algae reside, implying 615 higher concentrations at the very top of the snow. Simulated albedo impacts at 0.2 < λ < 0.35 should be discounted due to our lack of pigment m i data in this spectral range (Sect. 3.3.5). Unsurprisingly, larger cells cause greater reduction in albedo when comparing fixed number concentrations (Fig. 8b). Normalizing to mass, however, as we do for other LAC, would show that smaller algae absorb more per unit mass (not shown). Fig. 8d shows snow spectra for idealized heavy-load cases where only one type of pigment is present in the cell (as in Fig. 3h), included simply to highlight that unique spectral features can arise 620 when the pigments are blended in different proportions. The two chlorophyll pigments create reflectance peaks at slightly different wavelengths near λ = 500 nm, for example, and some combinations of carotenoids and chlorophyll cause distinctively green or red snow, perhaps similar to the "watermelon snow" that is caused by C. nivalis and often observed by backcountry enthusiasts. Algae, even more so than other LAC, tend to be concentrated near the very top (e.g., < 1 cm) of the snow, suggesting that use of the multi-layer model with thin layers near the surface is necessary to properly resolve the influence of 625 algae on albedo (Cook et al., 2017(Cook et al., , 2020. We also remind readers that this should be considered an experimental technique to represent snow algae, and that more observational and model development work is needed on this topic, including the exploration of techniques such as the SNICAR-GO formulation (Cook et al., 2020;Williamson et al., 2020) for representing highly non-spherical glacier algae. below, we include model scenarios limited to the resolution and extent of measurements, as well as additional scenarios that apply speculative but reasonable assumptions. Unless otherwise stated, model base case parameters are applied.
We utilize data from the following studies: 1985-1986 and 1990-1991, respectively, across a spectral range of 0.31 − 2.5 µm. We apply the diffuse spectral albedo 640 and standard deviations listed in Table 6 of that study, which represent average data from the two sites and show a narrow range of variability. In the base model scenario, we apply the vertical profiles of grain size and snow density reported in Tables 3 and 5 of that study. The measurements of grain size were conducted visually with a reported accuracy of about 25 µm. We also assume a sulfate-coated BC mass mixing ratio of 0.3 ppb, as reported, which has negligible impact on albedo.

645
- Hudson et al. (2006): This study focuses on snow bidirectional reflectance, but also includes spectral albedo measurements under diffuse sky conditions from Dome C, Antarctica, across a spectral range of 0.35 − 2.5 µm. We use the data shown in Figure 6 of that study, which represent an average of five measured spectra. Vertical grain size measurements were not included, but the authors state that the surface snow was composed of grains with radii of 50-100 µm, with little temporal variability. We therefore apply single-layer model scenarios with r e of 50 and 100 µm. We also apply the 650 reported BC mixing ratio of 3 ppb. Surface grain sizes were retrieved by comparing measured spectral albedo across the 1.03 µm absorption feature with a 655 lookup table of modeled reflectances from spherical ice grains (Nolin and Dozier, 2000). For consistency, we therefore also apply model spheres when using the reported grain sizes, though we include additional scenarios with hexagonal plates. We apply the reported BC and dust mixing ratios measured from each site, assuming the BC to be sulfate-coated and the dust to be fine Saharan dust. The reported dust concentrations are sufficiently low to have little impact on our analysis.

660
- Hadley and Kirchstetter (2012): Controlled laboratory experiments were performed in this study with varying amounts of flame-generated soot incorporated into aqueous suspensions that became mixed with quasi-spherical ice grains. We thus apply spherical ice grains in our model comparison with this dataset. Spectral albedo measurements, limited to λ = 400 − 900 nm, were made on 5 cm thick snow samples with direct zenith illumination, and scaling factors were applied by the authors to derive semi-infinite snow albedo. BC mass mixing ratios up to 1680 ppb were generated in 665 snow with multiple ice effective grain sizes, and we present the measured spectra shown in Figure 1 of Hadley and Kirchstetter (2012). The k a of the soot was estimated to be 15 m 2 g −1 at λ = 532 nm, higher than our standard soot, so we apply our more absorptive sulfate-coated soot for this comparison, with implications discussed below.
- Brandt et al. (2011): This study also created artificial snow infused with commercial soot, but did so outdoors on a 75 m 2 field. The frozen water droplets were quasi-spherical with an estimated effective grain size of 60 µm. Spectral albedo 670 measurements over λ = 0.35 − 2.5 µm were conducted under diffuse (cloudy) illumination on snow with and without a substantial load (∼ 2500 ppb) of BC. The "Aquablack-162" soot applied in this study had a smaller and narrower size distribution than the standard BC used in SNICAR, so for this comparison we created a separate species of BC with r m = 65 nm, σ g = 1.3, and k a = 6.0 m 2 g −1 at λ = 550 nm, matching the reported specifications by Brandt et al. (2011). (We also include this version of BC in the optics library.) A second artificial snowpack with a different type of 675 soot was presented in this study but was not rigorously characterized and showed noisier albedo measurements, so we do not study it here.
- Aoki et al. (2000): These authors measured spectral albedo of natural snowpacks in Hokkaido, Japan, coincident with vertical profiles of snow grain size, density, and impurity concentrations. The bulk impurity concentrations were not partitioned by species, though mean particle absorptivity was derived. Mixing our standard BC and fine Saharan dust 680 optical properties to match this absorptivity suggests a mean BC percentage by mass of ∼ 0.5 − 1.5%, though with considerable uncertainty. We assume a 98.5%/1.5% partitioning of dust and coated BC in our default scenarios for comparison, but find in sensitivity studies that larger BC fractions improve the spectral agreement. We apply the 3-layer vertical profiles of snow grain size, density, and impurity loads shown in Figure  - : Field measurements of spectral snow albedo (0.35 < λ < 1.5 µm) and vertical profiles of snow grain size, dust, and black carbon concentrations were conducted in this study throughout the 2013 spring melt season at the Senator Beck study site in the San Juan Mountains of Colorado, an area that periodically experiences large dust deposition events. We select three profiles for comparison: One with relatively clean snow before much melting had 690 occurred (13 April), one with moderate dust loads after a deposition event (28 April), and one with very heavy surface dust burdens after an additional deposition event and substantial melt had occurred (3 May). The vertical profiles of impurities, snow grain size, and density imposed in the model included measurements in 10 layers at 3 cm resolution and an 11th deep layer. Snow grain size measured with contact spectroscopy represents a sphere-equivalent grain size, and thus we apply spheres for this comparison. We apply the Colorado dust optical properties and particle size distribution 695 unique to this environment (Section 3.3.3, , along with the direct or diffuse light conditions that coincided with the measurements. For the 13 and 28 April cases we show a second model simulation with a thin surface layer composed of smaller or larger snow grains to demonstrate that improved comparison is achievable. Furthermore, top-layer dust concentration is set to that of the second-from-top, dustier layer in this sensitivity simulation for 28 April. 5.1 Clean snow 700 Fig. 9 shows model-measurement comparisons of snow with extremely low LAC content. The diffuse sky comparison with Grenfell et al. (1994) also includes 2-layer models of spheres and hexagonal plates with a very thin 0.25 mm surface layer.
Because measured grain size was not vertically-resolved to better than 5 mm, Grenfell et al. (1994) demonstrated how an unresolved but plausible thin surface layer of fine-grained snow could substantially improve model-measurement agreement in the NIR. Our model sphere scenario is identical to that utilized by Grenfell et al. (1994) and is included for reference. Fig.   705 9a shows that even better agreement at λ > 2.0µm can be achieved with a hexagonal plate scenario. All scenarios slightly underpredict a single measurement at 1.9 µm. All scenarios also produce slightly higher albedo than measured at λ < 0.4 µm.
An alternative explanation for model low bias in the 1.4 − 1.9 µm region is that the NIR refractive index data used as model input are biased. Indeed, Carmagnola et al. (2013) and Dumont et al. (2021) show that use of alternative NIR refractive index data (Grundy and Schmitt, 1998) yields better agreement with snow reflectance measurements from Greenland and France,710 without the need to assume thin, fine-grained surface layers in the model.
The model comparison with Hudson et al. (2006) (Fig. 9b) applies hexagonal plates with grain size not well-constrained by independent measurements. In addition to the single-layer semi-infinite cases with r e of 50 and 100 µm, we also include a 3-layer model with r e = 75 µm (top 0.25 mm), r e = 100 µm (0.25 − 0.50 mm depth), and r e = 175 µm below. This scenario was chosen again to illustrate that excellent agreement can be achieved across the NIR when very thin model layers are applied 715 near the surface. Measurements from Hudson et al. (2006) at λ < 0.4 are slightly higher than those from Grenfell et al. (1994), resulting in better model-measurement agreement in the UV. Model curves are slightly too low at λ = 1.9−2.0 µm and slightly high at λ = 1.2 − 1.35 µm, though the latter could probably be remedied with a 4-or 5-layer model.
The clear-sky clean-snow case from Casey et al. (2017), depicted in Fig. 9c, shows strong agreement with the matching spherical-grain scenario from SNICAR-ADv3 and a hexagonal plate scenario with slightly larger effective grain size than the 720 sphere r e retrieved by Casey et al. (2017). The agreement, however, is imperfect in the visible and UV spectrum, where the measurements show lower albedo than the model. This discrepancy may result from uncertainty in the measurements, which were of directional reflectance and required an anisotropic correction factor to determine albedo. We expect truly pristine snow to have higher visible and UV albedo, akin to what was measured by Hudson et al. (2006), unless the blue/UV ice absorption coefficients applied from Picard et al. (2016) and Warren and Brandt (2008) are biased low. Another source of uncertainty 725 when evaluating albedo under clear-sky conditions is the fraction of diffuse light, especially at short wavelengths, that was present for the measurements. Generally lacking this information, we assume unidirectional irradiance for comparisons against the clear-sky cases, though this is an inaccurate assumption in the blue spectrum where Rayleigh scattering is important.  (Fig. 10a-c), which are of natural snow with BC mass mixing ratios of 2400, 3300, and 490 ppb, are generally excellent. The spectral shape and magnitude of modeled albedo in the perturbed part of the spectrum (λ < 1.0 µm) are near-perfect in panels (b) and (c), but do not perfectly match the measured albedo in Fig. 10a, which shows snow with a slightly more brownish hue than that produced by the standard SNICAR BC. Smaller size distributions of BC can be more absorptive in the blue than red (as with the Aquablack described below), creating a spectral response of the nature 735 shown, but BC size distributions were not reported. The hexagonal plate model curves in Figs. 10b and 10c include a very thin (0.25 mm) surface layer with grain size selected to optimize agreement. The single-layer sphere cases are identical to those specified by Casey et al. (2017).

Snow with Black Carbon as the Dominant Source of LAC
The comparison against Hadley and Kirchstetter (2012), for snow r e of 55, 65, and 110 µm, is shown in Fig. 10d-f. In general, this comparison is also quite promising throughout the visible spectrum and for the entire range (0 − 1680 ppb) of BC 740 mixing ratios explored. We note, however, that there is a discrepancy between the estimated k a of the soot applied in that study (15 m 2 g −1 at λ = 532 nm) and that of the coated BC applied in SNICAR (11.5 m 2 g −1 ). If a form of BC with k a = 15 m 2 g −1 were applied in SNICAR, the model would overestimate the BC impacts on snow albedo, with albedo reductions roughly 30% greater than shown. A key finding from the measurements of Hadley and Kirchstetter (2012) is that the albedo perturbation from BC increases with increasing snow grain size, to a degree that is in-line with the model results presented here and confirming 745 earlier theoretical arguments for this phenomenon . wavelengths. Agreement is reasonable in this case, but Fig. 10h shows that improved agreement results when more detailed knowledge of the impurity properties is incorporated. First, Brandt et al. (2011) note that the city water supply used to create the snow "contained a small amount of a brown absorber (perhaps rust or humic acid)". Hence instead of applying 12 ppb of BC for the clean snow case, we apply 100 ppb of brown carbon, which produces an improved spectral match with the clean snow case. Second, when we instead apply the Aquablack BC properties described above, the modeled albedo of the contam-755 inated snow takes on a slightly more brownish hue, leading to a near-perfect match with measurements. Finally, Brandt et al.
(2011) speculated that the hydrophillic soot they used may have resided inside the ice grains, though they could not confirm this with observations and found that modeling with an external mixture of BC, as done here, produced good albedo agreement.
Internally-mixed BC in the model would produce a larger albedo impact (He et al., 2017), worsening the albedo comparison. The snow produced by Hadley and Kirchstetter (2012) may also have contained internally-mixed BC because the soot was 760 incorporated into an aqueous suspension, but this is also unconfirmed.

Snow with Dust and Black Carbon
BC and dust constitute the dominant sources of LAC in most natural snowpack, though their relative proportions vary widely in space and time. Measurements by Aoki et al. (2000) of natural snow with moderate impurity loads are shown in Fig. 11a-c.
When we combine a spherical grain model, as used by Aoki et al. (2000), with the same r e profiles they applied, we find 765 similarly excellent agreement with measurements at λ > 1.5µm and slight model overprediction at shorter NIR wavelengths, similar to the comparison they presented with a different radiative transfer model. Assuming hexagonal plates with larger grain sizes, which are encouragingly more similar to the r e values derived through image processing by Aoki et al. (2000), also produces good agreement in the NIR, especially near the 1.3 µm ice absorption feature. Adding a 0.25 mm fine-grain layer produces excellent agreement throughout the NIR, depicted with the yellow curve in Fig. 11a. (For economy, we omit 770 this hypothetical layer in panels (b) and (c), though it also reduces model bias in those comparisons). Using the looselyconstrained BC/impurity mass fraction of 1.5% produces an excessive visible albedo for the scenes shown in Figs. 11a and 11b, but reasonable albedo for that shown in Fig. 11c. The scene in Fig. 11b is too brown and too dark to be replicated with our standard coated BC and dust, constrained with measured impurity loads, but by assigning 15% of the impurity mass to brown carbon we achieve an excellent match. This is merely illustrative, as Aoki et al. (2000) do not mention brown carbon, though 775 its presence can't be ruled out either. It is also quite possible that more absorptive (brown) dust was present than represented in our model. The observation of ∼ 20× higher impurity loads at the snow surface than at depth, as exhibited in Figs. 11a and 11c, demonstrates the importance of accurately resolving the vertical distribution of impurity load, especially after melt consolidation has occurred (Doherty et al., 2013).
The snow measured by  was contaminated primarily by dust, in widely varying amounts, and 780 secondarily by black carbon. Using the Colorado dust properties and measured vertical impurity distributions in our modeling produces good agreement in visible albedo for the 13-April case, insufficient darkening in the 28-April case, and excellent agreement in the extreme dust case of 3-May (Fig. 11d-f). The measured dust mixing ratio in the second-from-top snow layer for the 28-April case was three times larger than that of the top layer. A sensitivity study with a new 1 cm surface layer darkened with this larger dust concentration and a larger snow grain size produces good agreement (red curve in Fig. 11e). This 785 is merely illustrative, but possibly within the range of measurement uncertainty, considering the different footprints of spectral albedo and snow-pit excavations used for vertical profiles. Modeled albedo in the NIR is overpredicted in the two most dusty cases (28-April and 3-May).  also observed this phenomenon and discuss possible explanations including under-estimated snow grain size due to dust masking of the 1.03 µm ice absorption feature and complex interactions associated with internal dust-ice mixtures that are not represented in SNICAR. Our demonstration of improved NIR agreement with grain 790 size adjustments alone for the two April cases points to potential issues with measurements of grain size and shape, while our inability to resolve the NIR bias in the dusty 3-May case without wrecking the visible agreement suggests a need to further investigate the potential role of mixing state described by . Overall, however, we find the agreement between measurements and modeling across this large range of dust-induced darkening, with visible albedo ranging from 0.31 − 0.86, to be encouraging.

6 Conclusions
We have presented the formulation of a publicly-available model and accompanying library of optical properties that are used to simulate the spectral albedo of snow, dependent on many variables including the content of light-absorbing constituents (LAC).
Types of LAC included with SNICAR-ADv3 are black carbon, brown carbon, mineral dust, volcanic ash, and snow algae, with variants associated with coatings, particle size, mineralogy (for dust), and pigment content (for algae). SNICAR-ADv3 unifies 800 numerous improvements and capabilities that have been introduced to the model through separate studies, including application of the adding-doubling two-stream solver (Dang et al., 2019) and representations of non-spherical ice particles (He et al., 2017), carbon dioxide snow (Singh and Flanner, 2016), snow algae (Cook et al., 2017), different types of dust (Polashenski et al., 2015;Wolff et al., 2009) and volcanic ash , new ice refractive indices (Picard et al., 2016), and SZA-dependent surface spectral irradiances. For representing natural snow with unknown ice grain shape, we recommend 805 applying one of the three available non-spherical ice grain shapes (spheroids, hexagonal plates, or Koch snowflakes), which have smaller scattering asymmetry parameters than spheres and therefore reduce the radiative penetration depth and exposure of sub-surface LAC (Dang et al., 2016). The use of hexagonal plates instead of equal-SSA spheres reduces the simulated albedo impact of black carbon in single-layer simulations by ∼ 24%.
Compared with spectral albedo measurements of clean snow, SNICAR-ADv3 performs well across the solar spectrum, and 810 arguably better when non-spherical ice shapes are used, but in most cases a very thin (∼ 0.25 mm) surface layer composed of fine-grain snow must be introduced to match observations across the near-infrared spectrum. No observational studies, to our knowledge, have resolved snow grain size at this vertical resolution. Alternatively, uncertainty in ice refractive indices may be large enough to reconcile model-measurement discrepancies between wavelengths of 1.4 and 1.9 µm without the need to assume thin surface layers Dumont et al., 2021), suggesting a need for additional measurements of ice 815 properties in this part of the spectrum. Compared to spectral albedo measurements of snow laden with black carbon, dust, and (in one case) brown carbon, SNICAR-ADv3 also performs well when reasonable assumptions are made. A substantial source of uncertainty in some comparisons, however, is the optical properties of the LAC present in the snow being studied, as these vary substantially with source provenance, particle size, and atmospheric processing. LAC internally-mixed in weakly-absorbing media, including (e.g.,) sulfate and ice, are predicted to generally absorb ∼ 20 − 60% more solar energy per unit mass than 820 externally mixed LAC, but the mixing states of particles in snow are also rarely known. We see no systematic directional bias in the snow albedo reduction simulated by SNICAR-ADv3 when applying the standard LAC included in the optical property library, suggesting that these properties are reasonable to apply when site-specific LAC properties are unknown. A companion study (Whicker et al., 2021) extends SNICAR-ADv3 to represent glacier ice. Other target areas for model improvement include representation of liquid water and snow roughness , as well as experimental verification of simulated impacts 825 from snow algae and volcanic ash. It is likely that algae modeling techniques will need to be adapted. We hope that the webbased model, multi-layer source code, and accompanying optical property library of SNICAR-ADv3 prove useful for future research, education, and model development efforts.
Appendix A: The delta-Eddington adding-doubling solution We refer readers to Briegleb and Light (2007) for a detailed derivation of the delta-Eddington adding-doubling multi-layer 830 two-stream solution. Here, we present only the essential equations applied in SNICAR-ADv3.
Defining f = g 2 , the layer bulk optical properties (Eqs. 4-6) are transformed with delta-scalings to account for forward scattering as: The following intermediate variables are applied in the two-stream solution, with Eqs A7 and A8 dependent on the cosine of the incident zenith angle (µ 0 ): The layer reflectance to direct-beam radiation is: The layer transmittance to direct radiation is: The layer reflectance and transmittance to diffuse radiation are calculated by integrating the respective direct quantities over µ: where eight Gaussian angles are used for the numerical integration.
Properties of individual layers are then combined through the adding-doubling approach to produce bulk properties for mul-860 tiple layers. The reflectance and transmittance of two stacked layers, with layer 1 overlying layer 2 and direct-beam illumination on layer 1 are: T 12 (µ 0 ) = e −τ * 1 /µ0 T 2 (µ 0 ) + T 1 (µ 0 ) − e −τ * 1 /µ0 + e −τ * 1 /µ0 R 2 (µ 0 )R 1 T 2 1 − R 1 R 2 (A14) The bulk properties for these two combined layers under diffuse illumination of layer 1 are: The combined transmittances for illumination from below are identical whereas the reflectance for illumination from below, i.e. by upward scattered radiation that is assumed to be diffuse, is equivalent to Eq. A15 but with reversed layer subscripts.
The cumulative layer interface properties for downward and upward propagating radiation are solved through loops that integrate from the top-down and bottom-up of the column. These interface quantities are: -R up (µ 0 ): the reflectance of the entire column below the interface to downwelling direct radiation from above -R up : the reflectance of the column below the interface to diffuse radiation from above -R dn : the reflectance of the column above the interface to upwelling diffuse radiation from below -T dn (µ 0 ): the bulk transmittance between the model top and the interface with respect to direct radiation incident on the 880 model top -T dn : the transmittance between the model top and the interface to diffuse radiation With these quantities defined, the downwelling and upwelling fluxes at each layer interface (F ↓ dr , F ↑ dr ), normalized to the unit direct flux incident at the model top, are, respectively: The downwelling and upwelling fluxes at each layer interface (F ↓ df , F ↑ df ), normalized to the unit diffuse flux incident at the model top, are, respectively: Finally, the albedo is derived from fluxes at the top interface: Code and data availability. The web-based single-layer version of SNICAR-ADv3 can be run at: http://snow.engin.umich.edu. Matlab source code for the multi-layer model, and the accompanying netCDF library of optics files are available at: https://github.com/mflanner/ SNICARv3. A frozen branch of the code used in this manuscript is associated with the doi: 10.5281/zenodo.5176213. The scripts and data used to generate all plots are archived at: https://doi.org/10.5281/zenodo.5707933. A Python implementation of the multi-layer model,
Author contributions. MGF wrote the manuscript, assimilated the code updates, and created the web-based model. JBA performed model evaluations against spectral albedo measurements, JMC created the original representation of snow algae, CD and CSZ incorporated the adding-doubling solver, CH designed and incorporated the representation of non-spherical ice particles, XH provided top-of-atmosphere spectral irradiances, DS incorporated carbon dioxide snow and Martian dust optical properties, SMS provided Colorado dust optical properties, CAW merged the adding-doubling solver into our new code base, and CSZ and MGF created the original SNICAR model. All authors reviewed the manuscript.
Competing interests. The authors declare no competing interests.