Enhanced Snow Absorption and Albedo Reduction by Dust‐Snow Internal Mixing: Modeling and Parameterization

We extend a stochastic aerosol‐snow albedo model to explicitly simulate dust internally/externally mixed with snow grains of different shapes and for the first time quantify the combined effects of dust‐snow internal mixing and snow nonsphericity on snow optical properties and albedo. Dust‐snow internal/external mixing significantly enhances snow single‐scattering coalbedo and absorption at wavelengths of <1.0 μm, with stronger enhancements for internal mixing (relative to external mixing) and higher dust concentrations but very weak dependence on snow size and shape variabilities. Compared with pure snow, dust‐snow internal mixing reduces snow albedo substantially at wavelengths of <1.0 μm, with stronger reductions for higher dust concentrations, larger snow sizes, and spherical (relative to nonspherical) snow shapes. Compared to internal mixing, dust‐snow external mixing generally shows similar spectral patterns of albedo reductions and effects of snow size and shape. However, relative to external mixing, dust‐snow internal mixing enhances the magnitude of albedo reductions by 10%–30% (10%–230%) at the visible (near‐infrared) band. This relative enhancement is stronger as snow grains become larger or nonspherical, with comparable influences from snow size and shape. Moreover, for dust‐snow external and internal mixing, nonspherical snow grains have up to ~45% weaker albedo reductions than spherical grains, depending on snow size, dust concentration, and wavelength. The interactive effect of dust‐snow mixing state and snow shape highlights the importance of accounting for these two factors concurrently in snow modeling. For application to land/climate models, we develop parameterizations for dust effects on snow optical properties and albedo with high accuracy.


Introduction
Mineral dust, as one of the most common and abundant aerosols by mass, plays a critical role in the Earth and climate system (Ginoux et al., 2012;IPCC, 2013;Shao et al., 2011). It directly interacts with atmospheric radiation via light scattering and absorption, resulting in important direct radiative effect (Kok et al., 2017;Mahowald et al., 2014;Miller et al., 2004). It also affects cloud formation and precipitation by acting as cloud condensation nuclei (CCN) and ice nuclei (IN) (Creamean et al., 2013;DeMott et al., 2010;Huang et al., 2014). Moreover, mineral dust can further reduce snow albedo and accelerate snow/glacier melting upon deposition onto snow-/ice-covered regions (Lee et al., 2017;Qian et al., 2015;Skiles et al., 2018). Observations have shown substantial dust deposition and induced snow albedo reductions over the global cryosphere, including East Asian seasonal snowpack Wang et al., 2017), North American mountains Skiles & Painter, 2017), European Alps (Di Mauro et al., 2015;Dumont et al., 2017), the Tibetan Plateau Zhang et al., 2018), the Arctic (Doherty et al., 2010;Dumont et al., 2014), and the Antarctic (McConnell et al., 2007). For example, Skiles et al. (2012) showed daily mean radiative forcing of up to 214 W m -2 caused by dust deposition on snow based on 6 years' observations over Colorado mountains (United States), which leads to significant snow cover loss and accelerated snow melting. Zhang et al. (2018) found up to 10 W m -2 radiative forcing due to dust in snow over the Tibetan Plateau based on measured dust concentrations, which could reduce snow cover duration by up to 4 days. Thus, it is important to understand and accurately quantify dust-snow interactions and associated albedo reductions in order to better estimate dust-induced snow albedo forcing and subsequent hydroclimatic effects.
Previous studies have made great effort to physically model the effects of light-absorbing aerosols, including black carbon (BC) and dust, on snow albedo mainly by assuming aerosols mixed outside spherical snow grains (i.e., external mixing) (Aoki et al., 2011;Flanner et al., 2007;Tuzet et al., 2017;Warren & Wiscombe, 1980). On this basis, Dang et al. (2015) developed simplified parameterizations for snow albedo reduction caused by BC and dust externally mixed with snow spheres. Nevertheless, observations have shown that nonspherical snow grains are common in reality (e.g., Dominé et al., 2003;Erbe et al., 2003). A number of studies have recently investigated the impact of nonspherical snow grains in snow albedo modeling (e.g., Dang et al., 2016;He, Takano, Liou, Yang, et al., 2017;Kokhanovsky & Zege, 2004;Libois et al., 2013;Liou et al., 2014;Räisänen et al., 2017), indicating higher pure snow albedo but smaller aerosolinduced snow albedo reduction for nonspherical snow grains compared with spherical grains. Moreover, aerosols such as BC have been found to be mixed inside snow grains (i.e., internal mixing) in many cases (e.g., Aoki et al., 2000;Magono et al., 1979), which substantially enhance snow albedo reduction relative to BC-snow external mixing (Flanner et al., 2012;He et al., 2014). Further combining snow grain nonsphericity and BC-snow internal mixing reveals important interactive effects on snow albedo reduction (He, Liou, Takano, Yang, et al., 2018, He, Flanner, et al., 2018Saito et al., 2019). However, these studies only focused on BC-snow interactions, while how snow nonsphericity and dust-snow internal mixing influence snow albedo reduction has not been investigated. This highlights an urgent need to understand and quantify their roles in dust-snow interactions.
Direct snowpack observations, albeit rather limited, have shown evidence on the existence of dust-snow internal mixtures (e.g., Hörhold et al., 2012;Kumai, 1976;Spaulding et al., 2011). In general, aerosols tend to mix externally with snow grains through dry deposition and/or below-cloud scavenging, while aerosolsnow internal mixtures can form through in-cloud scavenging (Flanner et al., 2012). First, as efficient ice nuclei, dust particles can enter into ice crystals via several heterogeneous (i.e., contact, immersion, condensation, or deposition) ice nucleation processes (DeMott et al., 2010;Hoose & Möhler, 2012). Second, though initially insoluble, dust particles can be coated with soluble species (e.g., sulfate and nitrate) during atmospheric transport and become efficient cloud condensation nuclei to form liquid cloud droplets (Levin et al., 1996;Li & Shao, 2009). The dust-containing cloud droplets can further freeze homogeneously under very low temperature (lower than about -40°C) and/or go through the riming process to form dustcontaining ice crystals . Furthermore, ice crystal aggregation, accretion, and riming processes during ice growth and precipitation formation periods may lead to multiple dust inclusions in one ice crystal by merging multiple dust-containing cloud hydrometeors into a single one, while snowflakes formed from the growth of dust-containing ice crystals via the Wegener-Bergeron-Findeisen process may have only single dust inclusion. Besides, dust-snow mixing state could also be altered by snowpack processes (e.g., melting, sublimation, sintering, and refreezing) after snowfall and dust deposition (Flanner et al., 2012;Skiles & Painter, 2017), though observational constraints are rather limited currently. Therefore, it is imperative to account for dust-snow internal mixing in estimating the impact of dust deposition on snow albedo. particles like dust (Flanner et al., 2012). Recently, Liou et al. (2014) developed a stochastic aerosol-snow albedo model (SASAM) that explicitly resolves aerosols internally mixed with nonspherical snow grains based on a geometric-optics surface-wave (GOS) approach (Liou & Yang, 2016, and references therein). The model has been successfully used for BC-snow albedo studies , He, Liou, Takano, Yang, et al., 2018 and shows consistent results with the dynamic effective medium approximation method in the BC-snow case (He, Flanner, et al., 2018). This provides a powerful tool to investigate dustsnow internal mixing.
Therefore, in the present study, we extend the SASAM model capability to account for dust-snow internal mixing and investigate its impact. To our knowledge, this is the first study to explicitly quantify the combined effects of dust-snow internal mixing and snow grain nonsphericity on snow single-scattering properties and albedo (section 3). We also develop a set of new dust-snow parameterizations for land/climate modeling applications (section 3). We further discuss radiative implications for dust-snow internal mixing (section 4) and sensitivities to dust refractive index and size distribution (section 5). Finally, we describe how to apply our parameterizations and associated uncertainties (section 6) and conclude our findings (section 7).

Model Descriptions
We use a multilayer SASAM originally developed by Liou et al. (2014) and improved subsequently by He, Takano, Liou, Yang, et al., (2017), . It explicitly simulates snow grains of different shapes internally/externally mixed with multiple polydisperse aerosol particles and associated snow optical properties and albedo. The model has been primarily applied to BC particles in previous studies (e.g., He et al., 2014, He, Takano, Liou, Yang, et al., 2017, He, Liou, Takano, Yang, et al., 2018. In this study, we further extend the model capability to resolve dust-snow interactions, particularly nonspherical snow grains internally mixed with size-resolved dust particles. Specifically, prescribed snow grain structures are generated in a three-dimensional coordinate system, with dust particles randomly distributed inside/outside snow grains using a stochastic procedure . The number of dust particles mixed with snow grains is determined by dust and snow densities and sizes as well as dust concentration in snow. The spectral singlescattering properties of dust-snow mixtures are computed based on the geometric-optics surface-wave approach (He et al., 2016;Liou et al., 2011;Takano et al., 2013), which uses a Monte Carlo photon tracing technique to simulate light scattering and absorption processes for dust-snow mixtures. For an ensemble of randomly orientated snow grains, their single-scattering properties are further averaged over all directions. The computed single-scattering properties (i.e., single-scattering albedo, asymmetry factor, and extinction efficiency) of dust-snow mixtures are subsequently used to derive snow albedo through the adding/doubling radiative transfer scheme (Takano & Liou, 1989) in snowpack. The number and thickness of snow layers in the model are adjustable based on research needs. More details can be found in Liou et al. (2014).

Simulations and Parameterizations
In the present study, we consider dust internally/externally mixed with snow grains with volume-equivalent sphere radii (R v ) of 100, 500, and 1,000 μm for four typical shapes (Figure 1; see also He, Takano, Liou, Yang, et al., 2017 for details), including sphere, spheroid, hexagonal plate, and Koch snowflake, which generally represent the observed key snow shape characteristics (Dominé et al., 2003;Erbe et al., 2003). These observations suggested that fresh snow grains with relatively small sizes tend to be snowflakes or hexagonal shapes, while aged snow grains with relatively large sizes are close to spheroids or spheres. However, we note that observations of snow grain shapes across a wide size range under different snowpack conditions are still limited, highlighting the necessity of more extensive and systematic snow shape measurements to better constrain snow albedo modeling. We use the Warren and Brandt (2008) database for spectral snow refractive indices, with an ice density of 0.917 g cm -3 (Warren & Wiscombe, 1980). Based on the observed range of dust content in snow (e.g., Huang et al., 2011;Painter et al., 2012;Wang et al., 2017), we account for dust mass concentrations (C dust ) of 0, 1, 10, 50, 100, 500, and 1,000 ppm (or μg g -1 ) to demonstrate low to high dust pollution in snow, with a dust density of 2.5 g cm -3 (Zender et al., 2003). The spectral dust refractive indices are from Dang et al. (2015) based on a comprehensive compilation of available observations, except that the imaginary refractive indices at wavelengths of >1 μm are from d' Almeida et al. (1991). We also conduct sensitivity simulations using different dust imaginary refractive indices to quantify the associated impacts and uncertainties (see section 5.1). We assume a lognormal dust size distribution with a geometric mean diameter of 0.65 μm and a geometric standard deviation of 2.0 (equivalent to an effective radius of 1.1 μm), which represents dust from large-scale transport (Formenti et al., 2011;Maring et al., 2003) and may be smaller than that from local soil (Kok, 2011). The dust size effects on snow optical properties and albedo are further quantified through sensitivity simulations (see section 5.2).
We simulate spectral (0.2-5.0 μm with a 0.01-μm interval) albedo of clean and dust-containing snow based on the SASAM model (section 2.1), assuming a homogeneous semi-infinite snowpack with a solar zenith angle of 49.5°(i.e., insolation-weighted mean condition for sunlit Earth hemisphere). The effects of solar zenith angle on snow albedo have been quantified and parameterized in several previous studies (e.g., Marshall, 1989;Saito et al., 2019). We further integrate the spectral snow optical properties and albedo into broadband values through weighted averages using the incoming solar spectra under typical midlatitude wintertime clear sky and overcast cloud conditions following Flanner et al. (2007).
Based on the rigorous SASAM model simulations (hereinafter reference calculations), we develop a set of parameterizations for dust-induced changes in snow single-scattering properties and albedo at different wavelength bands for application to land and climate models. Four widely used band types are considered in this study, including the broadband type with two bands (visible and near-infrared (NIR)), the Community Land Model (CLM) band type (hereinafter CLM bands; Oleson et al., 2004), the Fu-Liou Radiative Transfer Model band type (hereinafter Fu-Liou bands;Fu & Liou, 1993), and the Rapid Radiative Transfer Model (RRTM) band type (hereinafter RRTM bands; Mlawer & Clough, 1997). We discuss uncertainties and limitations associated with model simulations and parameterizations in section 6.

Results and Discussions
3.1. Effects on Snow Single-Scattering Properties 3.1.1. Dust-Snow Internal Mixing Previous studies have well documented the single-scattering properties and albedo of clean snow and associated influences from snow grain and snowpack properties (e.g., Flanner et al., 2007;Liou et al., 2014;Warren, 1982;Warren & Wiscombe, 1980). Here, we focus on dust-induced changes in snow singlescattering properties (single-scattering albedo (ω), asymmetry factor, and extinction efficiency) and albedo. We find that dust-snow internal mixing has negligible effects on snow asymmetry factor and extinction efficiency, which is consistent with previous studies He, Takano, Liou, Yang, et al., 2017) and thus not discussed here. However, compared with pure snow, the snow single-scattering coalbedo (1-ω, defined as one minus single-scattering albedo) and hence snow absorption are substantially enhanced by Figure 1. Demonstration of dust-snow internal (i.e., dust mixed inside snow grains; top row) and external (i.e., dust mixed outside snow grains; bottom row) mixing for four typical snow grain shapes used in this study, including sphere, spheroid, hexagonal plate, and Koch snowflake.
dust-snow internal mixing mainly at wavelengths <~1.0 μm, with larger enhancements for higher dust concentrations ( Figure 2). Figure 3 quantitatively shows the dust-induced spectral snow single-scattering coalbedo enhancement (E 1-ω ; defined as the single-scattering coalbedo of dust-containing snow divided by that of pure snow). The enhancement is the strongest at wavelengths of 0.2-0.4 μm and decreases sharply at wavelengths of >0.4 μm, which reduces to 1.0 (i.e., no enhancement) at wavelengths longer than~1.0 μm, because of the strong dust absorption but weak snow absorption at shorter wavelengths. This is similar to the spectral feature of BC-snow internal mixing but has much smaller enhancement than that caused by BC with the same Figure 2. Spectral single-scattering coalbedo from reference calculations for clean (black lines) and dust-contaminated snow (colored lines) with a volume-equivalent sphere radius of 100 μm and four snow grain shapes, including sphere (first column), spheroid (second column), hexagonal plate (third column), and Koch snowflake (fourth column), for dust-snow internal (top row) and external (bottom row) mixing. Figure 3. Spectral snow single-scattering coalbedo enhancement (E 1-ω ; see text for definition) caused by dust-snow internal (left panel) and external (right panel) mixing relative to pure snow, based on reference calculations for different dust concentrations in snow (indicated by different colors). The lines and shaded areas are, respectively, the mean and one standard deviation (1σ) uncertainty range based on calculations using different snow grain shapes (i.e., sphere, spheroid, hexagonal plate, and Koch snowflake) and snow volume-equivalent sphere radii (i.e., 100, 500, and 1,000 μm). Note that the uncertainty (shaded area) is small, due to the small effects of variabilities in snow grain shape and size on E 1-ω (as opposed to the large effects of dust content variabilities on E 1-ω ). To highlight the effective dust effect and its spectral details, results at wavelengths of >1.7 μm are not shown due to negligible dust effects. See Figures S1-S2 for results of each calculation case. concentration in snow (He, Takano, Liou, Yang, et al., 2017), due to the weaker absorption ability of dust per mass. The dust-induced E 1-ω varies by several orders of magnitude with different dust mass concentrations in snow (C dust ), with up to 10 2 and 10 5 for C dust of 1 and 1,000 ppm, respectively (Figures 3 and S1). We further find very small effects of variabilities in snow grain size and shape on E 1-ω as opposed to the large effect of dust content variabilities on E 1-ω (Figure 3), which agrees with the findings in the BC-snow case (He, Takano, Liou, Yang, et al., 2017).
Moreover, the E 1-ω caused by dust-snow internal mixing shows a quantitative nonlinear relationship with C dust and can be parameterized through the following empirical formulation: where C dust has the unit of ppm and a i (i = 1 -3) is the wavelength-dependent parameterization coefficient (independent of snow shape and size). This parameterization is in very good agreement with the reference calculation based on the rigorous SASAM model simulations at different wavelength bands (Figures 4 and  S3), with the absolute normalized mean bias (NMB) of <0.06%, the normalized root-mean-square-error (RMSE) of ≤0.1, and the coefficient of determination (R 2 ) of >0.997. For direct application to land/climate models with explicit snowpack radiative transfer schemes, we provide the parameterization coefficients for both clear sky and overcast cloud conditions at the CLM, Fu-Liou, and RRTM wavelength bands in Table 1. Note that we only parameterize the E 1-ω at wavelength bands shorter than~1.2 μm, due to negligible dust effects on snow absorption (i.e., E 1-ω very close to 1.0) at longer wavelengths. Results at wavelength bands of >~1.2 μm are not shown due to negligible dust effects on E 1-ω . The red lines and shaded areas are, respectively, the mean and one standard deviation (1σ) uncertainty range based on reference calculations using different snow grain shapes (i.e., sphere, spheroid, hexagonal plate, and Koch snowflake) and volumeequivalent sphere radii (i.e., 100, 500, and 1,000 μm). Note that the red solid lines and blue dashed lines are almost overlapped in all cases because of very close values, and the uncertainty (shaded area) is rather small due to very small snow shape and size effects on E 1-ω as opposed to the large dust effect (see also

. Dust-Snow External Mixing
Compared to dust-snow internal mixing, we find very similar spectral impacts of dust-snow external mixing on snow single-scattering coalbedo ( Figure 2), with the major snow absorption enhancement at wavelengths <~1.0 μm. The enhancement (E 1-ω ) in spectral single-scattering coalbedo caused by dust-snow external mixing relative to pure snow has a similar pattern as shown in the case of dust-snow internal mixing but with a smaller magnitude by up to a factor of 2 (Figures 3 and S2). Following the case of dust-snow internal mixing, we also develop a parameterization for the E 1-ω caused by dust-snow external mixing based on the same formulation (equation (1)). We provide the parameterization coefficients (a i in equation (1)) for the CLM, Fu-Liou, and RRTM wavelength bands under clear sky and overcast cloud conditions in Table 1. The parameterization shows good accuracy at different wavelength bands ( Figures S4 and S5), with the NMB of <6.5%, the normalized RMSE of ≤0.3, and the R 2 of >0.98.

Effects on Snow Albedo 3.2.1. Dust-Snow Internal Mixing
Because of the enhanced snow absorption caused by dust-snow internal mixing (section 3.1.1), the dustpolluted snow albedo decreases markedly at wavelengths of <1.0 μm compared to pure snow, with stronger reductions for larger snow grain sizes and dust concentrations (Figures 5 and S6). The albedo reduction reaches the maximum at short wavelengths (0.2-0.4 μm) and decreases quickly at longer wavelengths to about 0 beyond~1.0 μm, with a faster decrease for a larger dust concentration and/or snow grain size ( Figure 6). This is a result of the dominant pure snow absorption at wavelengths of >1.0 μm (Warren, 1982). We find that the maximum spectral snow albedo reduction caused by dust-snow internal mixing is more than doubled for snow grain radii (R v ) of 1,000 μm relative to that for R v of 100 μm, with reduction of up to about 0.05 (0.2) and 0.2 (0.5) for R v of 100 (1,000) μm with C dust of 10 and 100 ppm, respectively. For fixed C dust and R v , spherical snow grains lead to the strongest albedo reduction at wavelengths of <1.0 μm, while the nonspherical counterparts have up to~35% weaker reduction depending on snow size, dust concentration, and wavelength, with the smallest reduction for Koch snowflakes ( Figure 6). This is because nonspherical snow grains tend to have smaller asymmetry factors and thus weaker forward scattering (Dang et al., 2016). . Spectral snow albedo from reference calculations for clean (black lines) and dust-contaminated snow (colored lines) with a volume-equivalent sphere radius (R v ) of 100 μm and four snow grain shapes, including sphere (first column), spheroid (second column), hexagonal plate (third column), and Koch snowflake (fourth column). Results for dust-snow internal (top row) and external (bottom row) mixing are shown. See Figure S6 for results with a R v of 1,000 μm. Figure 6. Spectral snow albedo reduction caused by dust-snow internal (top row) and external (bottom row) mixing from reference calculations for snow volumeequivalent sphere radii (R v ) of 100 μm (dashed lines) and 1,000 μm (solid lines) and four grain shapes, including sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange). Results for dust concentrations in snow (C dust ) of 10 ppm (first column), 100 ppm (second column), and 1,000 ppm (third column) are shown.
Further analysis shows that the snow albedo reduction caused by dust-snow internal mixing can be quantitatively and nonlinearly related to dust concentration and snow grain size and shape, following an empirical formulation derived from laboratory measurements (Hadley & Kirchstetter, 2012). This formulation is originally proposed for BC-induced snow albedo reduction and has been successfully used to develop BC-snow parameterizations (He, Liou, Takano, Yang, et al., 2018;Saito et al., 2019). We find that this formulation is also applicable to dust-induced snow albedo reduction (Δα; defined as pure snow albedo minus dustpolluted snow albedo) for semi-infinite snowpack as follows: where C dust has the unit of ppm, b i (i = 1 -3) is the parameterization coefficient depending on snow grain shape and wavelength, R 0 (100 μm) is the reference snow grain radius, and R e (unit: μm) is the snow grain effective radius (i.e., specific projected-area-equivalent sphere radius) defined by the snow grain volume (V snow ) and projected area (A snow ). This definition of snow grain effective radius is widely used in snow/ice crystal modeling and satellite retrievals (e.g., Jin et al., 2008; and can be quantitatively converted to other commonly used snow effective radii, including the volume-equivalent sphere radius, the projected-area-equivalent sphere radius, and the specific surfacearea-equivalent sphere radius. For a broader application of the present parameterization, the detailed conversion among these snow effective radii can be found in He, Takano, Liou, Yang, et al. (2017) and Saito et al. (2019).
This parameterization (equations (2-4)) well reproduces the results from the reference calculations at different wavelength bands under both clear sky and overcast cloud conditions (Figures 7 and S7-S10) with very high accuracy (Figures 8 and S14), including the NMB of <4%, the RMSE of ≤0.02, and the R 2 of >0.99. We provide the parameterization coefficients (b i in equations (2 and 3)) at the visible-NIR bands (Table 2), CLM bands (Table S1), Fu-Liou bands (Table S2), and RRTM bands (Table S3) for application to land/climate models with bulk snow albedo schemes. Note that we only parameterize the Δα at wavelengths shorter thañ 1.5 μm for the CLM, Fu-Liou, and RRTM bands, due to negligible dust effects on snow albedo reduction (i.e., Δα very close to 0) at longer wavelengths.

Dust-Snow External Mixing
The spectral pattern of snow albedo reduction caused by dust-snow external mixing is similar to that caused by dust-snow internal mixing but with a smaller magnitude (Figures 5 and 6). The absolute difference in albedo reductions between dust-snow external and internal mixing tends to be larger for larger snow grain sizes. For example, compared with internal mixing, the maximum albedo reduction due to dust-snow external mixing is smaller by up to about 0.05 and 0.1 for R v of 100 and 1,000 μm, respectively, with a C dust of 1,000 ppm ( Figure 6). Furthermore, the variation in albedo reductions across different snow shapes is larger for dust-snow external mixing than internal mixing (Figure 6), for example, with the variation of up to about 0.15 (0.10) for external (internal) mixing at visible wavelengths with a C dust of 100 ppm and a R v of 1,000 μm. Interestingly, we find negative snow albedo reductions at wavelengths longer than~1.0 μm for cases with both very large dust concentrations (C dust ≥ 1,000 ppm) and snow grain sizes (R v ≥ 1,000 μm), which indicate higher albedo at these NIR wavelengths for dust-snow external mixtures than pure snow. This is because of the stronger scattering and hence higher albedo at NIR wavelengths (longer than~1.0 μm) for pure dust compared with pure snow and is consistent with the findings in Dang et al. (2015). They showed that when dust concentrations become large, dust can affect snow albedo by not only absorption but also scattering at the NIR wavelengths. This results in higher albedo for dust-snow external mixtures than pure snow, which also depends on snow grain sizes (e.g., C dust ≥1,000 ppm for R v of 1,000 μm and C dust ≥10,000 ppm for R v of 100 μm). However, we note that this feature is not seen for dust-snow internal mixing ( Figure 6).
Following the case of dust-snow internal mixing, we also develop a parameterization for snow albedo reduction caused by dust-snow external mixing based on the same formulation (equations (2-4)). The  (2-4)) for different snow grain shapes (i.e., sphere, spheroid, hexagonal plate, and Koch snowflake) and volumeequivalent sphere radii (i.e., 100, 500, and 1,000 μm). Results at different wavelength bands are shown, including the visible (0.3-0.7 μm) NIR (0.7-5.0 μm) bands (first column), the Community Land Model (CLM) bands (second column), the Fu-Liou bands (third column), and the Rapid Radiative Transfer Model (RRTM) bands (fourth column). Also shown are the normalized mean bias (NMB), the root-mean-square-error (RMSE), the coefficient of determination (R 2 ), and the 1:1 ratio lines (dashed lines). See Figure S14 for results under the overcast cloud condition.  (Figures 7 and 8 and S11-S14), with the NMB of <7.3%, the RMSE of ≤0.025, and the R 2 of >0.99. We provide the parameterization coefficients (b i in equations (2 and 3)) for the visible-NIR bands (Table 2), CLM bands (Table S1), Fu-Liou bands (Table S2), and RRTM bands (Table S3). We further compare this parameterization of dust-snow external mixing with that developed in Dang et al. (2015) for spherical snow grains, which reveals reasonably good agreement, except for slightly stronger albedo reductions for small snow sizes (e.g., R e = 100 μm) and weaker reductions for large snow sizes (e.g., R e = 1,000 μm) from our parameterization ( Figure S15). Such differences could be due to the uncertainties and biases in parameterization development and dust-snow albedo modeling in both studies (e.g., different radiative transfer models and dust imaginary refractive indices at wavelengths of >1 μm). Figure 9 shows the broadband snow albedo reduction caused by dust-snow internal and external mixing as a function of dust concentration in snow or snow grain effective radius (R e defined in equation (4)). We find a clear indication that snow albedo reductions are significantly enhanced by dust-snow internal mixing compared to external mixing, with stronger absolute increases in albedo reduction for larger dust concentrations and/or snow grain sizes. The relative increases in albedo reduction caused by internal mixing relative to external mixing tend to be larger for nonspherical snow grains than spherical grains, which are the largest for Koch snowflake, while larger snow grain sizes also lead to stronger relative increases (Figures 9 and S16-S17). Moreover, due to the weaker absolute albedo reduction caused by nonspherical snow grains than spherical grains (Figure 9 and sections 3.2.1 and 3.2.2), such snow nonsphericity effect could effectively offset the enhancement in albedo reduction induced by dust-snow  internal mixing. This highlights the importance of accounting for dust-snow internal mixing and snow grain shape concurrently in snow and climate modeling.

Dust-snow internal mixing
To quantitatively understand the albedo reduction enhancement caused by dust-snow internal mixing relative to external mixing, we further define an enhancement ratio (E Δα ) as the ratio of snow albedo reduction caused by dust-snow internal mixing (Δα int ) to that caused by external mixing (Δα ext ): Figures 10 and S18 show the enhancement ratios for different dust concentrations, snow sizes, and snow shapes. We find that the enhancement ratios are similar for various dust concentrations with fixed snow grain size and shape, which varies from 1.1 to 1.3 at the visible band and from 1.1 to 3.3 at the NIR band depending on snow size and shape (Table 3). In general, the enhancement ratio increases as snow grains become larger or nonspherical, for example, with the all wavelength (0.3-5.0 m) values of 1.29 (1.12) and 1.40 (1.19) for Koch snowflake (sphere) with R e of 100 and 1,000 μm, respectively ( Figure 10). As a result, the effect of snow grain shape on the enhancement ratio is comparable to that of snow grain size. We note that the albedo reduction enhancement caused by internal mixing relative to external mixing in dust-snow interactions is also comparable to that in BC-snow interactions (Flanner et al., 2012;He, Liou, Takano, Yang, et al., 2018). Thus, estimates of dust-induced snow albedo reductions without accounting for the reduction enhancement by internal mixing may result in important biases. The enhancement ratio (Table 3) in the present study provides a useful scaling factor for a first-order adjustment of dust-induced albedo reductions in snow/climate models to account for dustsnow internal mixing, since current models almost rely exclusively on the assumption of dust-snow external mixing.  2-4)) for dust-snow internal (solid lines) and external (dashed lines) mixing with a snow grain effective radius (R e , defined in equation (4)) of 100 μm. Results shown here are for snow sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange) at the visible (0.3-0.7 μm; first column), near-infrared (NIR, 0.7-5.0 μm; second column), and all wavelength (all, 0.3-5.0 μm; third column) bands. (Bottom row) same as Figure 9 (top row), but for snow albedo reduction as a function of snow grain effective radius (R e ) with a dust concentration in snow (C dust ) of 100 ppm. See Figure S16 for results under the overcast cloud condition and Figure S17 for results with a R e of 1,000 μm and a C dust of 1,000 ppm.

Radiative Implications
The dust-snow internal mixing strongly enhances snow albedo reduction compared with external mixing, which has important implications for the associated dust-induced snow albedo radiative effects.  . This estimate is generally consistent with their results from a sophisticated snow energy balance model . However, we find that dust-snow internal mixing can enhance the albedo reductions and hence surface radiative effects by about 20% (40%) for snow spheres (Koch snowflakes) compared with external mixing, which are equivalent to increases of 6.0-16.0 (11.7-24.4) W m -2 in dust-induced snow albedo radiative effects. Thus, the resulting absolute surface radiative effects are 47.0-93.3 (41.5-85.5) W m -2 for dust internally mixed with snow spheres (Koch snowflakes).
We note that the preceding calculation is not rigorous, which do not account for the variation and dynamic evolution of snowpack properties and environment conditions. Nevertheless, these estimates highlight that the dust-snow mixing state (internal versus external) together with snow grain shape (spherical versus nonspherical) plays a key role in quantifying dust-induced snow albedo radiative effects and needs to be accounted for in snow and climate modeling, particularly over snow-covered regions with substantial dust contamination. This could have further implications for water resource management and hydrological forecast (Painter et al., 2010). Figure 10. Comparisons of clear sky snow albedo reduction caused by dust-snow internal mixing (y-axis) and external mixing (x-axis) for snow sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange) with grain effective radii (R e , defined in equation (4)) of 100 (top row) and 1,000 μm (bottom row) at the visible (0.3-0.7 μm; first column), near-infrared (NIR, 0.7-5.0 μm; second column), and all wavelength (all, 0.3-5.0 μm; third column) bands. Circles are calculated results assuming dust concentrations of 0-1,000 ppm (with an interval of 10 ppm) in snow. Also shown are linear regression lines (solid lines) and regression slopes (enhancement ratio, E Δα defined in equation (5); colored numbers) for each snow grain shape (see also  Table 3). See Figure S18 for results under the overcast cloud condition.

10.1029/2019MS001737
Journal of Advances in Modeling Earth Systems 5. Sensitivities to Dust Refractive Index and Size Distribution

Dust Refractive Index
Previous studies (e.g., Balkanski et al., 2007;Dang et al., 2015;Wagner et al., 2012) have shown large variations in dust imaginary refractive indices, which strongly depend on dust composition (e.g., hematite/iron content), while the uncertainty in dust real refractive indices is much smaller. Thus, to assess the sensitivity of dust effects to its refractive indices, we conduct two sensitivity simulations by varying dust imaginary refractive indices. We use spectral dust imaginary refractive indices from the Wagner et al. (2012) database for dust with 2.7% hematite content as an upper limit (hereinafter the Wagner12 case) and the Balkanski et al. (2007) database for dust with 0.9% hematite content as a lower limit (hereinafter the Balkanski07 case). The dust imaginary refractive indices used in our standard simulations are approximately in the middle of this range. We find that the spectral patterns of snow single-scattering coalbedo enhancement (E 1-ω defined in section 3.1.1) in the two sensitivity simulations are similar to that in the standard simulation. However, the magnitude of E 1-ω is significantly affected by dust imaginary refractive indices, with differences of up to one order of magnitude between the Wagner12 and Balkanski07 cases for dust-snow external and internal mixing, particularly at visible wavelengths ( Figure 11). The Wagner12 case shows the highest E 1-ω , while the Balkanski07 case has the lowest E 1-ω , due to much smaller dust imaginary refractive indices in the Balkanski07 case. As a result, the Wagner12 case leads to much stronger dust-induced snow albedo reductions than the Balkanski07 case (e.g., by more than a factor of 2 at the visible band) for both dust-snow external and internal mixing with different snow grain shapes (Figures 12 and S19). The albedo reductions from the standard simulation lie in the middle of those from the two sensitivity simulations. The differences (absolute values <0.05) among these cases are much smaller at the NIR band (Figures 12 and S19).
Moreover, the snow albedo reduction enhancement (E Δα defined in equation (5)) caused by dust-snow internal mixing relative to external mixing in the two sensitivity simulations shows similar patterns and variations with snow grain size and shape as those in the standard simulation, where E Δα tends to be larger as snow grains become larger or nonspherical. The magnitude of E Δα only varies slightly (relative differences <10%) with dust imaginary refractive indices at the visible (0.3-0.7 μm) and all wavelength (0.3-5.0 μm) bands for different snow grain sizes and shapes (Figures 13 and S20-S22). Although the relative Figure 11. Spectral snow single-scattering coalbedo enhancement (E 1-ω ; see text for definition) caused by dust-snow internal (left column) and external (right column) mixing relative to pure snow, based on standard simulations (solid lines) and sensitivity simulations by varying dust refractive indices (top row) and size distributions (bottom row). Dust concentrations of 10 (red), 100 (blue), and 1,000 ppm (orange) are shown. Each line (solid, dashed, or dotted) is the mean value based on calculations using different snow grain shapes (i.e., sphere, spheroid, hexagonal plate, and Koch snowflake) and volume-equivalent sphere radii (i.e., 100 and 1,000 μm). Note that the variation caused by snow shape and size is very small for each line and thus not shown here (see also Figure 3 and text for detail).
In the top row, the dashed and dotted lines represent results using spectral dust imaginary refractive indices from the Wagner et al. (2012) and Balkanski et al. (2007) databases, respectively, compared with spectral dust imaginary refractive indices from the Dang et al. (2015) database used in standard simulations (solid lines). In the bottom row, the dashed and dotted lines represent results using dust effective radii (R dust ) of 2.5 and 5.0 μm, respectively, compared with a R dust of 1.1 μm used in standard simulations (solid lines).
differences in E Δα at the NIR band between the Wagner12 and Balkanski07 cases can be as large as 50%, the absolute NIR albedo reductions are much smaller compared with visible albedo reductions. Interestingly, we find that the albedo reduction enhancement (E Δα ) is generally larger for smaller dust imaginary refractive indices, with the strongest (weakest) enhancement for the Balkanski07 (Wagner12) case.

Dust Size Distribution
In addition to refractive indices, dust size distribution is also associated with large uncertainty and spatiotemporal variations, which is strongly affected by dust source and transport (Mahowald et al., 2014;Shao et al., 2011). In the standard simulation, we assume a lognormal dust size distribution with a geometric mean diameter of 0.65 μm and a geometric standard deviation of 2.0 (equivalent to an effective radius of 1.1 μm) Figure 12. Clear sky visible (0.3-0.7 μm; red) and near-infrared (NIR, 0.7-5.0 μm; blue) broadband snow albedo reduction as a function of dust concentration in snow for dust-snow internal (first and third rows) and external (second and fourth rows) mixing for sphere (first column), spheroid (second column), hexagonal plate (third column), and Koch snowflake (fourth column) with a snow volume-equivalent sphere radius of 100 μm, based on sensitivity simulations by varying dust refractive indices (top two rows) and size distributions (bottom two rows). In top two rows, the dashed and dotted lines represent results using spectral dust imaginary refractive indices from the Wagner et al. (2012) and Balkanski et al. (2007) databases, respectively, compared with spectral dust imaginary refractive indices from the Dang et al. (2015) database used in standard simulations (solid lines). In bottom two rows, the dashed and dotted lines represent results using dust effective radii (R dust ) of 2.5 and 5.0 μm, respectively, compared with a R dust of 1.1 μm used in standard simulations (solid lines). See Figure S19 for results with a snow grain radius of 1,000 μm.
typically for dust from long-range transport (Formenti et al., 2011;Maring et al., 2003), while dust near sources tends to have larger sizes (Kok, 2011). Thus, to assess the sensitivity of dust effects to its size distribution, we conduct two additional sensitivity simulations by varying dust size distributions. We assume lognormal dust size distributions with effective radii of 2.5 and 5.0 μm (i.e., geometric mean diameters of 1.5 and 3.0 μm and geometric standard deviations of 2.0 and 2.0), respectively, in the two sensitivity simulations (hereinafter the R dust 2.5 and R dust 5.0 cases).
We find that the snow single-scattering coalbedo enhancement (E 1-ω defined in section 3.1.1) caused by dustsnow internal/external mixing is affected by dust size distributions mainly at wavelengths smaller than 0.4 μm in the R dust 2.5 and R dust 5.0 simulations (Figure 11), with up to a factor of approximately three variations among the standard and two sensitivity simulations. The E 1-ω at the visible wavelengths decreases with increasing dust effective radii ( Figure 11). As a result, the dust-induced snow albedo reduction also decreases with increasing dust effective radii at the visible band for both dust-snow internal and external mixing (Figures 12 and S19), with negligible variations at the NIR band. This is consistent with the BC-snow interaction case , where BC-induced snow albedo reductions decrease with increasing BC effective radii. The variation (up to~40%) in visible albedo reductions tends to be smaller for larger snow sizes or nonspherical snow grains. However, the impact of dust size distribution on snow absorption and albedo reduction is much smaller than that of dust refractive indices.
Furthermore, the snow albedo reduction enhancement (E Δα defined in equation (5)) caused by dust-snow internal mixing relative to external mixing is only slightly affected by dust size distributions for different snow grain sizes and shapes (Figures 13 and S20-S22), with relative differences of <~10% at the visible and all wavelength bands among the standard and two sensitivity simulations (i.e., R dust 2.5 and R dust 5.0). We find that the albedo reduction enhancement (E Δα ) tends to be smaller for larger dust effective radii with relatively large snow radii (e.g., 1,000 μm), whereas it is the opposite with relatively small snow radii (e.g., 100 μm). This interesting feature requires further investigations.

Parameterization Application and Uncertainty
In this study, we have developed a set of parameterizations (equations (1-4)) for dust-induced changes in snow optical properties and albedo to account for internal and external mixing between dust and snow grains with different shapes for application to land/climate models. Specifically, for models with explicit snowpack radiative transfer schemes such as CLM/CESM, the snow single-scattering properties (i.e., single-scattering albedo, asymmetry factor, and extinction efficiency) are usually required. Thus, equation (1) can be used to account for the dust-induced changes in snow single-scattering albedo and hence the effect of dust-snow internal mixing (i.e., enhanced snow absorption), while dust-snow internal mixing has negligible impacts on snow extinction efficiency and asymmetry factor. The effect of dust-snow external mixing can be incorporated using either this study's single-scattering albedo parameterization for external mixing (equation (1)) or the weighted sums/averages of single-scattering properties from snow and dust, which have been used in previous models (e.g., Flanner et al., 2007;Warren & Wiscombe, 1980). The snow grain shape effect can be further included by using the parameterizations for snow single-scattering properties developed in He, Takano, Liou, Yang, et al. (2017). Recently, He, Flanner, et al. (2018) have successfully implemented similar parameterizations for internal mixing of BC and nonspherical snow grains into the SNICAR/CLM model, providing an example of feasible implementation. For models with bulk snow albedo schemes (i.e., no explicit snowpack radiative transfer computations) such as Noah-MP/WRF, equations (2-4) can be adopted for direct implementation to account for albedo changes in semi-infinite snowpack caused by dust internally/externally mixed with snow grains of different shapes.
We also realize that our parameterization development is associated with uncertainties and limitations. For example, we consider up to 1,000 ppm dust concentrations in snow that cover the majority of observed conditions (e.g., Huang et al., 2011;Qian et al., 2015;Wang et al., 2017), but extremely severe dust pollution in snow (>1,000 ppm) has also been measured occasionally , where the present parameterizations should be cautiously applied. Besides, dust properties (e.g., refractive index, size distribution, and particle morphology), strongly affected by its source, transport, and deposition, have high spatiotemporal heterogeneity and large uncertainty (Ginoux et al., 2012;Mahowald et al., 2014;Shao et al., 2011). We adopt the observed typical properties of dust particles from large-scale transport, which is reasonable for application in global climate models that rarely resolve subgrid processes of local soil dust. Furthermore, we assume monodisperse snow grain sizes, an assumption often used in previous modeling studies (Aoki et al., 2011;Liou et al., 2014;Warren & Wiscombe, 1980), whereas observations showed that polydisperse snow size distributions are common (Nakamura et al., 2001). Saito et al. (2019) recently developed snow albedo parameterizations that account for observed distributions of snow grain shapes and sizes, which tend to produce higher albedo than that using monodisperse snow spheres or Koch snowflakes (as used in this study) in the case of BC-snow internal mixing. Thus, it is important to incorporate observed snow size and shape distributions in modeling dust-snow interactions in future work.

Conclusions
In this study, we extended the SASAM capability to explicitly simulate dust internally/externally mixed with snow grains of different shapes and for the first time quantified the combined effects of dust-snow internal mixing and snow nonsphericity on snow optical properties and albedo. For application to land/climate models, we further developed a set of new parameterizations to account for dust-induced changes in snow optical properties and albedo. We also discussed important radiative implications for dust-snow internal mixing and associated uncertainties as well as sensitivities of dust effects to dust refractive indices and size distributions.
We found that dust-snow internal and external mixing significantly enhances snow single-scattering coalbedo and hence snow absorption relative to pure snow, primarily at wavelengths of <1.0 μm, with larger enhancements for internal mixing (relative to external mixing), higher dust concentrations, and shorter wavelengths. This enhancement varies by several orders of magnitude with different dust contents in snow but is very weakly affected by variabilities in snow grain size and shape. Due to the enhanced snow absorption, dust-snow internal mixing reduces snow albedo substantially at wavelengths of <1.0 μm relative to pure snow, with stronger reductions for higher dust concentrations, larger snow grain sizes, and spherical (relative to nonspherical) snow grain shapes. Compared with internal mixing, dust-snow external mixing shows similar spectral patterns of albedo reductions and effects of snow grain size and shape at wavelengths of <1.0 μm. Further analysis indicated that relative to external mixing, dust-snow internal mixing enhances snow albedo reductions by 10%-30% at the visible band and 10%-230% at the NIR band. This relative enhancement is generally stronger as snow grains become larger or nonspherical, with comparable influences from snow grain size and shape.
Moreover, for both dust-snow internal and external mixing, spherical snow grains lead to the strongest snow albedo reduction, while the nonspherical counterparts have up to~45% weaker reduction depending on snow size, dust concentration, and wavelength, with the smallest reduction for Koch snowflakes. The interactive effect of dust-snow mixing state (enhanced albedo reductions from internal mixing relative to external mixing) and snow nonsphericity (weakened albedo reductions from nonspherical grains relative to spherical grains) reveals the importance of accounting for these two factors concurrently in snow and climate modeling.
Further sensitivity studies showed that the dust impacts on snow absorption and albedo are significantly affected by dust imaginary refractive indices, with more than a factor of 2 variations at the visible wavelengths, whereas the sensitivity to dust size distribution is relatively weaker and mainly at wavelengths smaller than~0.4 μm, with less than~40% variations at the visible wavelengths. The dust-induced snow absorption enhancement and albedo reduction tend to be larger for larger dust imaginary refractive indices or smaller dust effective radii. However, the snow albedo reduction enhancement caused by dust-snow internal mixing relative to external mixing is only slightly affected by dust refractive indices and size distributions, with relative differences of <~10% at the visible and all wavelength bands.
For convenient application to land/climate models, we developed parameterizations for snow singlescattering coalbedo and snow albedo as a function of dust concentration in snow and/or snow grain size and shape, with high accuracies. Further, applying the present parameterizations to measured dust concentrations in Colorado snowpack as a case study, we suggested that dust-snow internal mixing can substantially (up to 40%) increase snow albedo reduction and associated surface radiative effects compared with dust-snow external mixing, depending on snow grain size, grain shape, and dust concentration. This highlights the essential role of dust-snow mixing state together with snow grain shape in estimating dust-induced snow albedo radiative effects, particularly over snow-covered regions with severe dust contamination.