Black carbon-induced snow albedo reduction over the Tibetan Plateau: uncertainties from snow grain shape and aerosol–snow mixing state based on an updated SNICAR model

We implement a set of new parameterizations into the widely used Snow, Ice, and Aerosol Radiative (SNICAR) model to account for effects of snow grain shape (spherical vs. nonspherical) and black carbon (BC)–snow mixing state (external vs. internal). We find that nonspherical snow grains lead to higher pure albedo but weaker BC-induced albedo reductions relative to spherical snow grains, while BC–snow internal mixing significantly enhances albedo reductions relative to external mixing. The combination of snow nonsphericity and internal mixing suggests an important interactive effect on BC-induced albedo reduction. Comparisons with observations of clean and BC-contaminated snow albedo show that model simulations accounting for both snow nonsphericity and BC–snow internal mixing perform better than those using the common assumption of spherical snow grains and external mixing. We further apply the updated SNICAR model with comprehensive in situ measurements of BC concentrations in the Tibetan Plateau snowpack to quantify the present-day (2000–2015) BC-induced snow albedo effects from a regional and seasonal perspective. The BC concentrations show distinct and substantial subregional and seasonal variations, with higher values in the non-monsoon season and low altitudes. As a result, the BCinduced regional mean snow albedo reductions and surface radiative effects vary by up to an order of magnitude across different sub-regions and seasons, with values of 0.7–30.7 and 1.4–58.4 W m−2 for BC externally mixed with fresh and aged snow spheres, respectively. The BC radiative effects are further complicated by uncertainty in snow grain shape and BC–snow mixing state. BC–snow internal mixing enhances the mean albedo effects over the plateau by 30–60 % relative to external mixing, while nonspherical snow grains decrease the mean albedo effects by up to 31 % relative to spherical grains. Based on this study, extensive measurements and improved model characterization of snow grain shape and aerosol–snow mixing state are urgently needed in order to precisely evaluate BC–snow albedo effects. Published by Copernicus Publications on behalf of the European Geosciences Union. 11508 C. He et al.: Black carbon-induced snow albedo reduction over the Tibetan Plateau

The Tibetan Plateau (TP), also known as the Third Pole, is covered by the largest mass of snow and ice outside the Arctic and Antarctic (Kang et al., 2010;Yao et al., 2012).It is the source region of major Asian rivers, providing fresh water for billions of people (Qin et al., 2006;Immerzeel et al., 2010).Meanwhile, because of its thermal heating, the TP has profound dynamical influences on the atmospheric circulation in the Northern Hemisphere and long been identified to be critical in regulating the Indian and East Asian monsoons (Manabe and Terpstra, 1974;Yeh et al., 1979;Yao et al., 2012).The TP is very sensitive to the changes in snow albedo and cover, which alter surface heat and water balances and further disturb the Asian hydrological cycle and monsoon climate (Kang et al., 2010).Observations have shown substantial BC concentrations in snow over the TP and suggested that BC deposition is an important driver of strong albedo reductions and accelerated glacier retreat in the region (Ming et al., 2008(Ming et al., , 2013;;Xu et al., 2009;Qu et al., 2014;Ji et al., 2015;Niu et al., 2017;Li et al., 2017;Zhang et al., 2018).Recent studies found that BC particles over the TP are primarily from South and East Asia, while long-range transport from northern mid-latitudinal source regions outside Asia also has non-trivial contributions (Kopacz et al. 2011;Lu et al., 2012;He et al., 2014a, b;Zhang et al., 2015;Li et al., 2016;Yang et al., 2018).
To estimate BC-induced snow albedo effects over the TP, previous studies often used observed BC concentrations in snow/ice as inputs to snow albedo models by assuming spherical snow grains and BC-snow external mixing (e.g., Ming et al., 2013;Jacobi et al., 2015;Schmale et al., 2017;Zhang et al., 2018).This simplified treatment of BC-snow interactions has been widely used in snow albedo modeling over various snow-covered regions (e.g., Warren and Wiscombe, 1980;Flanner et al., 2007;Aoki et al., 2011).However, snow grains are usually nonspherical (Dominé et al., 2003) and internally mixed with BC particles (Flanner et al., 2012) in real snowpack, which could significantly affect BC-snow albedo effects.For example, Kokhanovsky and Zege (2004) pointed out that substantial errors could occur if assuming spherical snow grains in albedo modeling.Dang et al. (2016) found that, compared with spherical snow grains, the nonspherical counterparts lead to higher pure snow albedo but smaller BC-induced albedo reduction for BC-snow external mixing.In addition, Flanner et al. (2012) showed that there could be up to 73 % of BC in global snowpack internally mixed with snow grains, which increases BCinduced albedo effects by up to 86 % relative to purely external mixing for spherical snow grains.Moreover, recent studies (He et al., 2014b(He et al., , 2018a;;Liou et al., 2014), combining both effects of snow nonsphericity and BC-snow internal mixing, revealed that the enhancement in snow albedo reduction caused by internal mixing can be weakened by snow nonsphericity effect.Therefore, ignoring these two critical factors in previous studies could lead to biased estimates of BC-induced snow albedo effects over the TP and elsewhere, which highlights the necessity of accounting for the two features together in snow albedo modeling and assessing the associated uncertainty.
In this study, we implement a set of new BC-snow parameterizations (He et al., 2017b) into the widely used Snow, Ice, and Aerosol Radiative (SNICAR) model (Flanner et al., 2007) to consider the effects of snow nonsphericity and BCsnow internal mixing.We further apply the updated SNICAR model with a comprehensive set of in situ measurements of BC concentrations in the TP snowpack to estimate the present-day (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) BC-induced snow albedo effects and associated uncertainties from snow grain shape (spherical vs. nonspherical) and BC-snow mixing state (external vs. internal) from a regional and seasonal perspective.To the best of our knowledge, this is the first attempt to quantify BC-snow albedo effects over the TP by taking into account the aforementioned two factors concurrently with observational constraints.We describe BC observations in the TP snowpack in Sect. 2. We implement the BC-snow parameterizations and evaluate model results in Sect.3. We further quantify and discuss the BC-snow albedo effects and associated uncertainties in Sect. 4. Finally, we present conclusions, implications, and future work in Sect. 5.

BC observations in the Tibetan snowpack
We collect available in situ observations of BC concentrations in snow/ice over the TP and surrounding areas during 2000-2015 from historical measurements (see Table S1 for summary).Although the features of BC concentrations at each site have been described in detail by previous observational studies, the present analysis seeks to summarize all these measurements in order to understand the regional and seasonal characteristics of BC pollution in the TP snowpack and more importantly to estimate the corresponding BC-snow albedo effects and associated uncertainties due to snow grain shape and BC-snow mixing state using an updated snow model (see Section 3).
For detailed analyses, we divide the entire TP and surrounding areas into six sub-regions (Fig. 1 27-34 • N, 95-105 • E), the central TP (CTP; 30-36 • N, 78-95 • E), and the Himalayas (HIMA).We note that NOTP represents the Tianshan region.Due to its proximity to the TP, we analyze it together with the TP snowpack in this study.Moreover, BC concentrations in the TP snowpack show distinct altitudinal and seasonal variations within each sub-region (Fig. 1a-f), with much larger values at relatively lower altitudes (< 5200 m a.s.l.) and in the non-monsoon season (October-May; Xu et al., 2009), compared with higher altitudes (> 5200 m a.s.l.) and the monsoon season (June-September; Xu et al., 2009), respectively.Thus, we conduct analyses according to different altitudes (above or below 5200 m a.s.l.) and seasons (monsoon or non-monsoon).In addition, for any observational site with multiple measurements during the same season, we average the measurements to represent the mean BC pollution condition for this site during the season.As a rather limited number of sites provide vertically resolved BC measurements throughout snowpack, we 11510 C.He et al.: Black carbon-induced snow albedo reduction over the Tibetan Plateau average BC concentrations throughout snow layers at these sites, which may introduce some uncertainties.
Figure 1a-f show that BC concentrations in snow are generally much higher during the non-monsoon period than during the monsoon period by up to one order of magnitude, except for NWTP and NOTP.This is because the four subregions (NETP, SETP, CTP, and HIMA) are dominated by the strong BC emissions in the non-monsoon season (particularly winter and spring) over South and East Asia (Lu et al., 2012;Zhang et al., 2015;Yang et al., 2018) and the efficient wet removal of BC in Asia in the monsoon season (Xu et al., 2009;He et al., 2014a).In contrast, the high concentrations during the monsoon period over NWTP and NOTP are primarily caused by the enrichment of BC via sublimation and/or melting of snow (Ming et al., 2009;Yang et al., 2015) and emissions from Central Asia and Middle East (Kopacz et al., 2011;Schmale et al., 2017).
Furthermore, BC concentrations are consistently larger at low altitudes (< 5200 m) than at high altitudes (> 5200 m) by a factor of 2-10 in each sub-region (Fig. 1a-f), which is consistent with previous studies (Ming et al., 2009(Ming et al., , 2013) ) that suggested that BC concentrations decrease with increasing elevations.Such altitudinal contrast in BC concentrations are maximal (with differences larger than one order of magnitude) over HIMA and SETP.This elevational dependence can be attributed to the stronger local emissions at lower elevations, the reduced efficiency of BC transport to higher elevations, and the higher temperature at lower elevations leading to stronger snow melting and hence BC enrichment in snow (e.g., Ming et al., 2013;Niu et al., 2017;Zhang et al., 2018).
Among the six sub-regions, the high-altitude areas in HIMA and SETP show the lowest BC concentrations (5-30 ppb) throughout the year (Fig. 1d-f), while NETP (with only low-altitude sites) during the non-monsoon season is most severely polluted by BC (∼ 4300 ppb).The results further indicate that BC concentrations in low-altitude areas across different sub-regions are comparable (190-450 ppb) during the monsoon season but are much more variable during the non-monsoon season (Fig. 1d-f).The variation of BC concentrations across the sub-regions is a result of combined effects of the aforementioned factors (e.g., regionally and seasonally dependent impacts from BC source, transport, removal, and snow aging).We note that the current observations over the TP are still rather limited spatially and temporally, leading to questions of representativeness and introducing uncertainty in the analysis.Thus, the large sub-regional, altitudinal, and seasonal heterogeneity of BC concentrations in the TP snowpack highlights an urgent need for extensive measurements.
3 Model description, implementation, and evaluation 3.1 SNICAR model Flanner et al. (2007) developed a multi-layer Snow, Ice, and Aerosol Radiative (SNICAR) model, which has been widely used for snowpack simulations globally.It is also coupled to global climate models (e.g., Community Earth System Model, CESM) to investigate effects of impurity contamination, snow grain properties, and snow aging on snowpack albedo.A detailed model description has been presented by Flanner et al. (2007) and implementation in CESM is described by Oleson et al. (2013).Here we briefly summarize the key model elements related to the present study.SNICAR simulates snowpack radiative transfer based on the theory from Wiscombe and Warren (1980) and the multi-layer twostream radiative transfer solution from Toon et al. (1989).It resolves vertical distributions of snow properties, impurity distributions, and heating throughout the snowpack column, as well as impact of underlying ground surfaces.The number of snow layers can be specified by users according to research objectives.The default SNICAR model assumes spherical snow grains and external mixing between impurities and snow grains.As inputs to radiative transfer calculations, the optical properties (extinction cross section (Q ext ), single-scattering albedo (ω), and asymmetry factor (g)) of snow grains and impurities, archived as lookup tables, are offline computed by the Mie theory based on particle size distributions and refractive indices.SNICAR utilizes clear-and cloudy-sky surface incident solar flux typical of mid-latitude winter.The input parameters for SNICAR include incident radiation type (direct or diffuse), solar zenith angle, number of snow layers with thickness, density, and grain effective radius in each layer, underlying ground albedo, and aerosol concentrations in snow.In this study, we use the stand-alone version of SNICAR (available at http://snow.engin.umich.edu/snicarcode/) and implement new parameterizations of snow nonsphericity and BC-snow internal mixing into it (see Sects. 3.2 and 3.3).The updated SNICAR model is available at https://github.com/EarthSciCode/SNICARv2.

Implementation of nonspherical snow grains
Previous studies commonly used an effective spherical snow grain with an equal volume to area ratio (i.e., equal surface area-weighted mean radius; hereinafter effective radius, R e ) to represent its nonspherical counterpart (e.g., Fu et al., 1999;Grenfell et al., 2005).The equal-effective-radius representation works well in computing extinction efficiency and single-scattering albedo but is inaccurate for asymmetry factor (Dang et al., 2016).To explicitly resolve snow grain shapes, Liou et al. (2014) have developed a stochastic snow albedo model based on a geometric-optics surfacewave (GOS) approach (Liou et al., 2011;He et al., 2015He et al., , 2016;;Liou and Yang, 2016).Further, He et al. (2017b) devel-oped a parameterization to account for snow nonsphericity effects on asymmetry factors for three typical grain shapes, including spheroid, hexagonal plate, and Koch snowflake (see Fig. 1 of He et al. 2017b).They parameterized the asymmetry factor (g ns ) of nonspherical snow grains as follows: where a i (i = 0-2) is the wavelength-dependent coefficient available in He et al. (2017b).R s (unit: µm) is equal to the snow effective radius (R e ) for spheroid or hexagonal plate, and R e /0.544 for Koch snowflake due to its complex concave shape.f s,x and f s,hex are the shape factors (i.e., ratio of R s of a nonspherical grain to that of an equal-volume sphere) of a nonspherical grain (x: spheroid, hexagonal plate, or Koch snowflake) and hexagonal plate, respectively.C g is the correction factor, and g hex is the asymmetry factor for hexagonal shapes computed as follows (Fu, 2007): where ω is the snow single-scattering albedo, and g is the asymmetry factor related to geometric reflection and refraction.b i and c i (i = 0-2) are the wavelength-dependent coefficients available in Fu (2007).AR is the snow aspect ratio (i.e., ratio of grain width to length).
Here we implement the He et al. (2017b) parameterization (Eqs.1-4) for snow asymmetry factor into SNICAR to account for nonspherical shapes.Due to the coarse spectral resolution (6 bands) of the parameterization, we further use a piece-wise shape-preserved polynomial interpolation method (Fritsch and Carlson, 1980) to interpolate the parameterized results into 470 bands (0.3-5 µm with a 10 nm resolution) used in SNICAR.The same interpolation method is also applied to implementing the single-scattering co-albedo parameterization for BC-contaminated snow (see Sect. 3.3).We use the extinction efficiency and single-scattering albedo of equal-effective-radii spheres for those of the nonspherical grains.
Figure 2a-c show the spectral snow asymmetry factors for different grain shapes based on the updated SNICAR model.Compared with spherical snow grains, nonspherical grains (particularly Koch snowflakes) result in up to ∼ 17 % smaller asymmetry factors at wavelengths < ∼ 3.0 µm, consistent with previous studies (Liou et al., 2014;Dang et al., 2016).We note that the results slightly (< 3 %) overestimate the asymmetry factors at two spectral peaks within 1.5-2.5 µm for spheroids with large sizes (R e ≥ 500 µm), due to parameterization uncertainties (He et al., 2017b).
As a result of the smaller asymmetry factors, nonspherical snow grains lead to weaker forward scattering and hence higher albedo relative to their spherical counterparts (Figs. 3 and S1 in the Supplement).We find up to about 2 % and 27 % higher albedo for Koch snowflakes in the visible (0.3-0.7 µm) and near-infrared (NIR, 0.7-5 µm) bands, respectively, compared to equal-R e spheres (Fig. 3d and e).These results show good agreement with the conclusions from previous studies (Wang et al., 2017;He et al., 2018a).The results also have important implications for snow grain size retrievals via the use of albedo models to match observed spectral reflectance.For example, Dang et al. (2016) and He et al. (2018a) suggested that if a nonspherical grain is simulated as a sphere, its effective size has to be scaled down by a factor of 1.2-2.5 to obtain the correct snow albedo.

Implementation of BC-snow internal mixing
Flanner et al. (2012) showed that the effect of BC-snow internal mixing can be equivalent to applying an enhancement ratio to BC absorption cross sections with the external mixing assumption and developed a lookup table for the enhancement ratio.Recently, He et al. (2017b) explicitly resolved the structures of BC-snow internal mixtures for different snow shapes and found that inclusions of BC increase snow single-scattering co-albedo (1-ω) and hence absorption but have negligible effects on snow asymmetry factor and extinction efficiency.They further parameterized the effect of internal mixing on 1-ω as follows: where E 1−ω is the co-albedo enhancement defined as the ratio of single-scattering co-albedo for BC-contaminated snow to that for pure snow, which is a function of BC mass concentration in snow (C BC , unit: ppb).d i (i = 0-2) is the wavelength-dependent parameterization coefficient available in He et al. (2017b).
Here we implement the He et al. (2017b) parameterization (Eq.5) for snow single-scattering co-albedo to account for BC-snow internal mixing.We note that the BC mass absorption cross section (MAC) at 550 nm used in He et al. ( 2017b) is 6.8 m 2 g −1 with a BC density of 1.7 g cm −3 and an effective radius of 0.1 µm.Thus, to obtain a BC MAC of 7.5 m 2 g −1 at 550 nm recommended by Bond and Bergstrom (2006), we adjust the BC size and density in this study.We assume a lognormal BC size distribution with a geometric mean diameter of 0.06 µm following Dentener et al. (2006) and Yu and Luo (2009) and a geometric standard deviation of 1.5 following Flanner et al. (2007) and Aoki et al. (2011).Then, we tune the BC density to 1.49 g cm −3 to match the MAC (7.5 m 2 g −1 ).The resulting BC size effect on E 1−ω is quantified using a parameterization developed by He et al. (2018b) as follows:   with where E 1−ω,R BC and E 1−ω,R BC =0.05 are the E 1−ω for a certain BC effective radius (R BC , unit: µm) and a R BC of 0.05 µm (reference case), respectively.k λ,R BC and f λ,R BC are empirical parameters relying on wavelength and BC size.m λ and n λ are wavelength-dependent coefficients available in He et al. (2018b).We should note that BC MAC could vary significantly in reality (e.g., from 2 to 15 m 2 g −1 at 550 nm) due to uncertainties from particle density, size, structure, and refractive index (Bond and Bergstrom, 2006).Thus, we use the recommended value (7.5 m 2 g −1 ) derived from a comprehensive review of measurements to reduce the potential uncertainty from BC MAC in this study.Compared with the current estimates, using a smaller BC MAC (e.g., 6.8 m 2 g −1 at 550 nm as used in He et al., 2017b) would lead to weaker BC-induced snow albedo reductions, the quantification of which is beyond the scope of this study and will be investigated in future work.In addition, the adjusted BC density and size used in the present study are still within the observed ranges, with 1.2-1.9g cm −3 for densities (Bond and Bergstrom, 2006;Long et al., 2013) as well as 0.01-0.15µm and 1.2-2.0 for geometric mean diameters and standard deviations (Bond et al., 2006), respectively.Figure 2d-f show the spectral single-scattering co-albedo of snow internally mixed with BC based on the updated SNICAR model.The strongest co-albedo enhancement (up to about 4 orders of magnitude for 1000 ppb BC) is in the visible band, with negligible effects at wavelengths > 1 µm.As a result of the enhanced snow absorption, snow albedo reduces about two-fold more due to internal mixing than external mixing (Figs. 4 and S2-S4 in the Supplement).In contrast, BC decreases snow albedo much less for nonspherical snow grains than spherical grains (Figs. 4 and S3-S4 in the Supplement), suggesting an important interactive effects of snow grain shape and BC-snow mixing state on snow albedo reductions.For example, BC-sphere external mixing leads to similar visible albedo reductions with BC-hexagonal plate internal mixing.This is consistent with our previous findings (He et al., 2018a).Although the internal mixing effect dominates at the NIR wavelengths (Fig. 4e), the NIR albedo reduction is a factor of 3-5 lower than the visible reduction.Thus, both snow nonsphericity and BC-snow internal mixing play comparably important roles in determining allwavelength albedo reductions (Fig. 4f).This highlights the significance of simultaneously accounting for these two factors in accurate estimates of BC-snow albedo effects.
Table 1.Parameter values used in SNICAR simulations when comparing with observed snow albedo (see Figs. 6 and 7).The observed snowpack properties are used in each case when they are available.Four types of snow shapes (sphere, spheroid, hexagonal plate, and Koch snowflake) and/or two types of BC-snow mixing (internal and external) are assumed in the simulations.Moreover, to cross-validate model results, we compare the simulated snow albedo and its reduction for BC-snow internal mixing using the He et al. (2017b) parameterization with those using the Flanner et al. (2012) lookup table.We find very good agreement (mean differences < 3 %) between the two schemes for different snow sizes and shapes (Figs. 5 and S5-S6), although the He et al. (2017b) parameterization leads to slightly stronger and weaker albedo reductions for higher (> 1000 ppb) and lower (< 1000 ppb) BC concentrations, respectively.Compared with the lookup table method, the newly implemented parameterization in this study can be applied to a wider range of snow grain size, shape, and BC concentration scenarios without sacrificing computational efficiency.

Pure snow albedo
We evaluated spectral pure snow albedo from SNICAR simulations by comparing with observations (Fig. 6) from laboratory measurements (Hadley and Kirchstetter, 2012), open-field experiments (Brandt et al., 2011), and field measurements in the Rocky Mountains (Painter et al., 2007) and at the South Pole (Grenfell et al., 1994).To conduct reasonable comparisons, we used the observed snow density, grain size, thickness, snowpack layer, direct or diffuse radiation, solar zenith angle, and underlying ground albedo in model simulations for each case (see Table 1 and Fig. 6 for details), except for underlying ground albedos in the Brandt et al. (2011) and Painter et al. (2007) cases and the grain size of the second snow layer in the Brandt et al. (2011) case because of unavailable measurements.Thus, we assumed black underlying grounds (albedo = 0) in the two cases, which has negligible effects on albedo estimates due to thick snow optical depths.In the Brandt et al. (2011) case, we further assumed an effective radius of 500 µm (typical for aged snow) in the second snow layer to make it optically semi-infinite, which is consistent with the observed condition.We also assumed four types of snow shapes (sphere, spheroid, hexagonal plate, and Koch snowflake) in the simulations to investigate shape effects, due to the lack of measurements.Wavelength (µm) Figure 6.Comparisons of spectral pure snow albedo from observations (black) and SNICAR simulations using observed snowpack properties (see Table 1 and text for details) and assuming sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange).
(a) Observations are obtained from laboratory measurements (Hadley and Kirchstetter, 2012).(b) Observations are obtained from open-field experiments in New York (Brandt et al., 2011).The effective radii (R e ) for each snow shape are obtained to best match observations at wavelengths of 1-1.3 µm.(c) Observations are obtained from field measurements over Rocky Mountains (Painter et al., 2007).(d) Observations are obtained from field measurements at the South Pole (Grenfell et al., 1994).
We find that model simulations generally capture the observed patterns of spectral snow albedo in all cases (Fig. 6).However, assuming spherical grains tends to underestimate snow albedo in the NIR band, while using nonspherical grains improves model results.For example, compared with the observations (Painter et al., 2007), simulations assuming snow spheres show a systematic underestimation of up to ∼ 0.1 at wavelengths > 0.7 µm, particularly at 1.0-1.2µm (Fig. 6c), while simulations assuming hexagonal plates well match the observations.Similarly, in the observational case of Grenfell et al. (1994), assuming hexagonal plates and Koch snowflakes substantially reduces model underestimates at 1.5-2.5 µm relative to assuming spheres, though leading to a slight overestimate at 0.9-1.3µm (Fig. 6d).In contrast, in comparison with the laboratory measurements from Hadley and Kirchstetter (2012), the spherical assumption works reasonably well, particularly for large sizes, with only slight (< 0.05) underestimates.This is because the snow grains created in those experiments tend to be spherical.Nevertheless, using spheroids and hexagonal plates in this case still leads to slightly better model results for large (R e = 65 and 110 µm) and small (R e = 55 µm) grain sizes, respectively (Fig. 6a).
In the observational case of Brandt et al. (2011), they determined snow effective sizes by matching model results with the measured NIR (1.0-1.3 µm) albedo.We find that assuming different snow shapes results in drastically different grain sizes retrieved by matching their measured NIR albedo (Figs.6b and 7d), with effective radii of 80 and 160 µm for spheres and Koch snowflakes, respectively.This implies the necessity of accounting for realistic grain shapes in snow grain size retrievals.
We note that model results in all cases show slight but consistent albedo overestimates around 400 nm compared with observations (Fig. 6), probably due to the uncertainty of ice refractive indices.In this work, we used ice refractive indices from the most widely used database (Warren and Brandt, 2008) obtained from measurements in the Antarctic, which shows a very low ice absorption coefficient around 400 nm.However, based on more recent measurements in Antarctic snow, Picard et al. (2016) found a much higher ice absorption coefficient around 400 nm than that from Warren and Brandt (2008), which suggested that the uncertainty in ice visible absorption is probably larger than generally appreciated.Therefore, the weak snow absorption caused by refractive indices used in this study could lead to the overestimates in modeled albedo around 400 nm.

BC-contaminated snow albedo
We further compared BC-contaminated snow albedo from SNICAR simulations with observations (Fig. 7) from laboratory measurements (Hadley and Kirchstetter, 2012), openfield experiments (Brandt et al., 2011;Svensson et al., 2016), and field measurements in the Arctic (Meinander et al., 2013;Pedersen et al., 2015).Similar to Sect.3.4.1,we used the observed BC concentration in snow, snow density, grain size, thickness, snowpack layer, direct or diffuse radiation, solar zenith angle, and underlying ground albedo in model simulations for each case (see Table 1 and Fig. 7 for details), except for the snow density in the Pedersen et al. (2015) case and the underlying ground albedo in the Meinander et al. (2013) case because of unavailable measurements.Thus, we assumed a typical fresh snow density of 150 kg m −3 in the former case and a black underlying ground (albedo = 0) in the latter case.Compared with assuming a black underlying ground, we find that using a non-black underlying ground albedo typically observed over the Tibetan Plateau (Qu et al., 2014) only leads to very small (< 5 %) relative differences in albedo calculations in the Meinander et al. (2013) case.Due to the lack of measurements, we further assumed BC internally or externally mixed with different snow shapes in the simulations to quantify the combined effects of BC-snow mixing state and snow grain shape.
Compared with the widely used assumption of BC externally mixed with spherical snow grains, we find that accounting for both internal mixing and snow nonsphericity improves model simulations (Fig. 7).For example, assuming BC-sphere external mixing leads to a systematical underestimate of polluted snow albedo for < 2000 ppb BC compared with the observations from Svensson et al. (2016), while assuming BC-Koch snowflake internal mixing reduces the model underestimate (Fig. 7b), with the normalized mean bias (NMB) and root-mean-square error (RMSE) decreasing from −0.04 to 0.01 and from 0.033 to 0.019, respectively.Similarly, in the observational case of Pedersen et al. (2015), simulations assuming BC-spheroid external mixing perform better than those assuming BC-sphere external mixing (Fig. 7a), reducing the NMB from −0.012 to −0.003 and RMSE from 0.028 to 0.025.Compared with the observations of Meinander et al. (2013), model results using spherical snow grains underestimate the spectral snow albedo contaminated by BC (Fig. 7c), regardless of model assumptions of BC-snow mixing state.Using nonspherical grains instead increases the simulated albedo and reduces model biases in this case, although it is still unable to fully capture the observed pattern (Fig. 7c).Considering that snow grains tend to be spherical in the observations from Hadley and Kirchstetter (2012), we assumed BC-sphere external/internal mixing in the comparisons.The model results with external mixing are systematically biased high, particularly for large BC concentrations (>110 ppb), while using internal mixing effectively reduces the albedo overestimates (Fig. 7e).As such, the observations fall between the results of external and internal mixing, suggesting a combination of partial external and internal mixing would best match the observations.Compared with the way of increasing BC MAC for BC-snow external mixing to reduce model overestimates in polluted snow albedo, which was used in Hadley and Kirchstetter (2012), the present study provides a physically based alternative (i.e., internal mixing) for model improvements.In fact, it is very likely that a large portion of BC is internally mixed with snow grains in the experiments of Hadley and Kirchstetter (2012), as they produced the BC-contaminated snow via freezing of aqueous hydrophilic BC suspensions.
We note that the snow grain sizes reported by the aforementioned field studies are retrieved by different methods, including matching snow model results with measured albedo spectra (Painter et al., 2007;Hadley and Kirchstetter, 2012;Pedersen et al., 2015) and visual estimates with tools (Grenfell et al., 1994;Meinander et al., 2013;Svensson et al., 2016) that are not equivalent to the snow effective size (i.e., surface area-weighted mean radius) defined in SNICAR.This could introduce uncertainties to snow albedo calculations and model-observation comparisons.

BC-snow albedo effects and uncertainties over the Tibetan Plateau
Based on the observed BC concentrations in snow (see Sect. 2), we applied the updated SNICAR model (see Sect. Figure 7. Comparisons of BC-polluted snow albedo from observations and SNICAR simulations using observed snowpack properties (see Table 1 and text for details).(a) Observations (x axis) are obtained from field measurements in the Arctic (Pedersen et al., 2015).Model results (y axis) for spheres (circles) and Koch snowflake (triangles) are shown as lower and upper limits for shape effects, along with BC-snow external (blue) and internal (orange) mixing.Also shown is the best case (red asterisks; BC-spheroid external mixing) that matches observations.(b) Observations (red asterisks; broadband albedo for 0.285-2.8µm) are obtained from open-field experiments in Finland (Svensson et al., 2016).Model results for spheres (circles) and Koch snowflake (triangles) are shown as lower and upper limits for shape effects, along with BC-snow external (blue) and internal (orange) mixing.(c) Observations (black lines) are obtained from field measurements in the European Arctic (Meinander et al., 2013).Model results assuming sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange) along with BC-snow external (dashed lines) and internal (solid lines) are shown.(d) Observations (black) are obtained from open-field experiments in New York (Brandt et al., 2011).BC is assumed to be externally mixed with snow spheres (blue), spheroids (red), hexagonal plates (green), and Koch snowflakes (orange).The effective radii (R e ) for each snow shape are obtained to best match observations at wavelengths of 1-1.3 µm.(e) Observations (circles) are obtained from laboratory measurements (Hadley and Kirchstetter, 2012).BC is assumed to be externally (solid lines) and internally (dashed lines) mixed with snow spheres.
reduction and associated surface radiative effects over the TP.We conducted albedo simulations at each observational site using the measured snowpack thickness and density (see Table S1 in the Supplement) concurrently with BC measurements.If the snow property measurements are missing at certain site, the data from nearby sites are used instead.We then computed the regional mean values by averaging across all sites within each sub-region and season.We used typical effective radii of 100 and 1000 µm for fresh and aged snow, respectively, to demonstrate snow aging and size effects.Due to the lack of measurements for snow grain shape and BC-snow mixing state, we considered eight simulation scenarios with the combination of four snow shapes (sphere, spheroid, hexagonal plate, and Koch snowflake) and two mixing states (internal and external).In the simulations, the underlying ground albedo over the TP is 0.1 at the visi-ble band (0.3-0.7 µm) and 0.2 at the NIR band (0.7-5 µm), following observations (Qu et al., 2014).We adopted a solar zenith cosine of 0.65 (i.e., an angle of 49.5 • ), which is equivalent to the insolation-weighted solar zenith cosine in the sunlit hemisphere.The effect of solar zenith angle on snow albedo can be approximated via changing snow effective size (Marshall, 1989).Previous studies (e.g., Aoki et al., 2003;Dang et al., 2016) indicated that the impact of snow shape and BC contamination decreases with an increasing solar zenith angle.Following Dang et al. (2017), we compute all-sky snow albedo via averages of clear-and cloudysky albedo weighted by cloud cover fraction.The mean cloud cover fraction and all-sky surface downward solar radiation in different sub-regions and seasons (see Table S2 1), including the northwestern TP (NWTP), north of the TP (NOTP), the northeastern TP (NETP), the southeastern TP (SETP), the central TP (CTP), and the Himalayas (HIMA).(a-c) Box plots of mean snow albedo reductions within each sub-region based on SNICAR simulations using the observed BC concentrations in snow (Fig. 1), snow thicknesses, and snow densities (see text for details).Results for altitudes > 5200 and < 5200 m are shown as left and right box plots within each sub-region, respectively, with circles and triangles indicating mean values.Model results assume BC externally and internally mixed with spheres, spheroids, hexagonal plates, and Koch snowflakes for fresh (blue, R e = 100 µm) and aged (red, R e = 1000 µm) snow.Each data point used for the box plot is the sub-regional average assuming a type of snow shape and BC-snow mixing, and hence the box plot indicates the variation caused by effects of snow shape and BC-snow mixing state.Note that some sub-regions only have BC observations at altitudes > 5200 or < 5200 m (see Fig. 1).(d-f) Same as (a-c), but for BC-induced all-sky surface radiative effects caused by the snow albedo reductions shown in (a-c).Calculations use the surface downward solar radiation and cloud cover fraction from the MERRA-2 reanalysis fields (see text and Table S2 in the Supplement for details).
search and Applications version 2 (MERRA-2) reanalysis meteorological fields (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/) with a spatial resolution of 0.5 • × 0.625 • .Figure 8a-c show the regional mean BC-induced snow albedo reductions in different sub-regions and seasons.The spatiotemporal distribution of albedo reductions generally follows that of BC concentrations in snow (Fig. 1d-f), with stronger albedo reductions in low-altitude areas and the non-monsoon period.We find that snow albedo decreases by a factor of 2-3 more for aged snow (Table S3 in the Supplement) than for fresh snow (Table 2), due to larger grain sizes for aged snow.This aging and size effect dominates the albedo reductions in most of TP sub-regions, particularly during the monsoon season (Fig. 8a-c).However, in severely polluted sub-regions including the lowaltitude areas of NETP, SETP, CTP, and HIMA during the non-monsoon season, the effects of snow grain shape and BC-snow mixing state are comparable to those of snow size/aging (Tables 2 and S3 in the Supplement).For example, BC-sphere internal mixing leads to an albedo reduction of 0.114 for fresh snow in low-altitude CTP during the nonmonsoon season, while BC-Koch snowflake external mixing leads to a reduction of 0.119 for aged snow.
Moreover, BC-snow internal mixing enhances the mean albedo reductions by 30-60 % (relative difference) across all the sub-regions and seasons, with similar enhancements for different snow shapes and sizes (Tables 2 and S3 in the Supplement).For example, assuming BC-sphere external mixing leads to an annual albedo reduction of 0.066 (0.164) for fresh (aged) snow in NETP, while the internal mixing counterpart results in a reduction of 0.095 (0.225).Our results are partially different from those in He et al. (2018a), which showed a stronger enhancement (relative difference) in albedo reduction caused by internal mixing for nonspherical grains than Table 2. Regional and seasonal mean BC-induced all-sky snow albedo reductions for fresh snow over the Tibetan Plateau during 2000-2015.See Table S3 for results of aged snow.), the northwestern Tibetan Plateau (NWTP), the northeastern Tibetan Plateau (NETP), the southeastern Tibetan Plateau (SETP), and north of the Tibetan Plateau (NOTP).Each sub-region is further divided into high (> 5200 m) and low (< 5200 m) altitudes.spherical grains, due to different environmental conditions and snow albedo models used in the two studies.We further find that nonspherical snow grains weaken the mean albedo reductions by up to 31 % relative to spherical grains in different sub-regions and seasons, with the strongest weakening for Koch snowflakes (Fig. 8a-c).The nonsphericity effect is smaller for aged snow compared with fresh snow (Tables 2   and S3 in the Supplement), consistent with our previous findings (He et al., 2018a).
Although the BC concentrations in the TP snowpack tend to dominate the regional and seasonal pattern of snow albedo reductions for fresh/aged snow (Figs.1d-f and 8a-c), the combined effects of snow grain shape and BC-snow mixing state can complicate the picture.For example, with the widely used assumption of BC externally mixed with snow spheres, the non-monsoon albedo reductions are 0.034 and 0.067 for high-altitude CTP and low-altitude SETP with BC concentrations of 332 and 1111 ppb in fresh snow, respectively.However, if BC particles were internally mixed with snow spheres in CTP and externally mixed with Koch snowflakes in SETP, the albedo reductions in the two areas would become the same (0.047), regardless of the substantially different BC concentrations.This points toward an imperative need for both extensive measurements and improved model characterization of snow grain shape and aerosol-snow mixing state for accurate quantification of BCinduced snow albedo reductions over the TP and elsewhere with strong heterogeneity of snowpack properties and contamination.
Figure 8d-f show the regional mean surface radiative effects caused by BC-induced snow albedo reductions, which vary from 0.7 to 11.2 W m −2 across different sub-regions during the monsoon season and from 1.2 to 30.7 W m −2 during the non-monsoon season for BC externally mixed with fresh snow spheres.The sub-regional variation increases to 1.4-37.7 and 3.5-58.4W m −2 for aged snow during the monsoon and non-monsoon periods, respectively (Tables 3 and  S4 in the Supplement).In general, the spatiotemporal distribution of surface radiative effects follows that of snow albedo reductions (Fig. 8a-f).The impacts of snow nonsphericity and BC-snow internal mixing on the surface radiative effects are similar to those on the albedo reductions discussed above.The maximum surface radiative effect over the TP can reach up to 45.4 (79.9)W m −2 in NETP during the non-monsoon season for BC internally mixed with fresh (aged) snow spheres (Tables 3 and S4 in the Supplement).The mean BC-induced snow albedo effects in the relatively clean TP areas (e.g., high-altitude HIMA and SETP) are comparable to those over the Arctic and North American snowpack (Dang et al., 2017;He et al., 2018a), while the effects in the contaminated TP areas (e.g., low-altitude HIMA, CTP, SETP, and NETP) are generally similar to those in the low-elevation Alps (Painter et al., 2013) and northern China snowpack (Wang et al., 2017).
Previous studies have shown accelerated snowmelt caused by BC-snow albedo effects in the TP.For example, Yasunari et al. ( 2010) estimated that BC-induced albedo reductions over Himalayan glaciers could result in an extra snowmelt of 1-7 mm day −1 during the melting/summer season.Qian et al. (2011) found a BC-induced snowmelt of up to 1.3 mm day −1 in late spring and early summer averaged over the entire TP.Our results further suggest that the uncertainty associated with snow shape and BC-snow mixing state could lead to a substantial variation in BC-induced albedo reduction and hence snowmelt, which has significant implications for runoff and water management in Asia.Accurate quantifications of the impact of snow grain shape and BCsnow mixing state on snowmelt and subsequent hydrologi-cal processes require interactive land surface and/or climate modeling, which will be investigated in future work.
We note that the present estimates of BC-induced snow albedo effects have uncertainties and limitations.For example, different techniques have been used to measure BC concentration in snow and ice, which may lead to discrepancies and inconsistency among observations and in modelobservation comparisons (Qian et al., 2015 and references therein).Additionally, BC measurements across the TP are from various sample types, such as surfaces of snowpack (with fresh and aged snow) and glacier (with both snow and firn and granular ice), which may introduce uncertainty to the understanding of BC contamination patterns (Zhang et al., 2017a;Li et al., 2018).In addition, in the model, we do not account for the vertical variability of BC and snow grain properties in the TP snowpack as well as some complex snowpack processes, including dynamic snow aging and melting, post-depositional enrichment, and melting water scavenging, which may exert non-trivial effects on BC-snow albedo effects (e.g., Flanner et al., 2007;Qian et al., 2014;Dang et al., 2017).These uncertainties associated with modeling and measurements may decrease the signal-to-noise ratio for the detection of BC effects on snow albedo, particularly in relatively clean regions with small BC-induced albedo reductions (e.g., < 0.01).Thus, improved and robust estimates require both accurate snow albedo modeling and snowpack measurements.

Conclusions, implications, and future work
We implemented a set of new BC-snow parameterizations into SNICAR, a widely used snow albedo model, to account for the effects of snow nonsphericity and BC-snow internal mixing.We evaluated model simulations by comparing with observations.We further applied the updated SNICAR model with a comprehensive set of in situ measurements of BC concentrations in the Tibetan Plateau (TP) snowpack (glacier) to quantify the present-day BC-induced snow albedo effects and associated uncertainties from snow grain shape and BCsnow mixing state.
Based on the SNICAR model updated with new BCsnow parameterizations, we found that nonspherical snow grains tend to have higher pure albedos but lower BCinduced albedo reductions compared with spherical snow grains, while BC-snow internal mixing substantially enhances albedo reductions relative to external mixing.Compared with observations, model simulations assuming nonspherical snow grains and BC-snow internal mixing perform better than those with the common assumption of snow spheres and external mixing.The results suggest an important interactive effect from snow nonsphericity and internal mixing, and highlight the necessity of concurrently accounting for the two factors in snow albedo and climate modeling.
We further applied the updated SNICAR model with comprehensive in situ observations of BC concentrations in snow and snowpack properties over the TP to quantify the presentday (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) BC-induced snow albedo effects.We found that BC concentrations show distinct sub-regional and seasonal variations.The concentrations are generally higher in the non-monsoon season and low-altitudes (< 5200 m) than in the monsoon season and high-altitudes (> 5200 m), respectively.The spatiotemporal distributions of snow albedo reductions and surface radiative effects generally follow that of BC concentrations.As a result, the BC-induced mean albedo effects vary by up to an order of magnitude across different sub-regions and seasons, with values of 0.7-30.7 (1.4-58.4)W m −2 for BC externally mixed with fresh (aged) snow spheres.
Moreover, the BC-snow albedo effects over the TP are significantly affected by the uncertainty in snow grain shape and BC-snow mixing state.We found that BC-snow internal mixing enhances the mean albedo effects by 30-60 % relative to external mixing across different sub-regions and seasons, while nonspherical snow grains reduce the albedo effects by up to 31 % relative to spherical grains.These effects become comparably important with the snow aging and size effect over polluted areas.Therefore, the combined effects of snow grain shape and BC-snow mixing state can complicate the spatiotemporal features of BC-snow albedo effects over the TP, with significant implications for regional hydrological processes and water management.
In summary, this study points toward an imperative need for improved measurements and model characterization of snow grain shape and aerosol-snow mixing state in order to accurately estimate BC-snow albedo effects.In future work, we will incorporate the new features of the updated SNICAR model into land surface and climate models, including CESM-Community Land Model (CLM) for global modeling and WRF-Noah-MP for regional modeling, to account for the effects of snow grain shape and aerosol-snow mixing state and to assess the associated uncertainties and hydrological feedbacks in global/regional climate system.

Figure 1 .
Figure 1.Observed BC concentrations in snow over the Tibetan Plateau (TP) during (a, d) monsoon, (b, e) non-monsoon, and (c, f) annual periods in 2000-2015 (see TableS1for details).(a-c) Spatial distributions of seasonal mean BC concentrations at altitudes > 5200 (circles) and < 5200 m (triangles) in six sub-regions, including the northwestern TP (NWTP), north of the TP (NOTP), the northeastern TP (NETP), the southeastern TP (SETP), the central TP (CTP), and the Himalayas (HIMA).(d-f) Box plots of observed BC concentrations in snow (shown in a-c) within each sub-region, with medians (middle bars), interquartile ranges (between 25th and 75th percentiles; boxes), and maxima/minima (whiskers) within ±1.5 × interquartile ranges.Some box plots are shrunk due to limited samples.Results for altitudes > 5200 and < 5200 m are shown as left and right box plots within each sub-region, respectively, with circles and triangles indicating mean values.Note that some sub-regions only have observations at altitudes > 5200 or < 5200 m.

Figure 4 .
Figure 4. (a-c) Spectral (0.3-1.5 µm) direct-beam albedo of semi-infinite snow layers with effective radii (R e ) of (a) 100, (b) 500, and (c) 1000 µm for pure snow (dotted lines), snow externally mixed with 100 ppb BC (dashed lines), and snow internally mixed with 100 ppb BC (solid lines) with shapes of sphere (blue), spheroid (red), hexagonal plate (green), and Koch snowflake (orange) based on the updated SNICAR model.The results for 1000 ppb BC and diffuse snow albedo are shown in Fig.S2in the Supplement.(d-f) Same as (a-c), but for broadband albedo reduction as a function of BC concentration in snow with R e of 100 µm at (d) visible (VIS, 0.3-0.7 µm), (e) near-infrared (NIR, 0.7-5 µm), and (f) all (0.3-5 µm) wavelengths.The results for snow with R e of 500 and 1000 µm and diffuse albedo reductions are shown in Figs.S3 and S4 in the Supplement, respectively.

Figure 8 .
Figure 8. Regional and seasonal mean BC-induced all-sky snow albedo reductions and surface radiative effects during (a, d) monsoon, (b, e) non-monsoon, and (c, f) annual periods in 2000-2015 over six Tibetan Plateau (TP) sub-regions (see Fig.1), including the northwestern TP (NWTP), north of the TP (NOTP), the northeastern TP (NETP), the southeastern TP (SETP), the central TP (CTP), and the Himalayas (HIMA).(a-c) Box plots of mean snow albedo reductions within each sub-region based on SNICAR simulations using the observed BC concentrations in snow (Fig.1), snow thicknesses, and snow densities (see text for details).Results for altitudes > 5200 and < 5200 m are shown as left and right box plots within each sub-region, respectively, with circles and triangles indicating mean values.Model results assume BC externally and internally mixed with spheres, spheroids, hexagonal plates, and Koch snowflakes for fresh (blue, R e = 100 µm) and aged (red, R e = 1000 µm) snow.Each data point used for the box plot is the sub-regional average assuming a type of snow shape and BC-snow mixing, and hence the box plot indicates the variation caused by effects of snow shape and BC-snow mixing state.Note that some sub-regions only have BC observations at altitudes > 5200 or < 5200 m (see Fig.1).(d-f) Same as (a-c), but for BC-induced all-sky surface radiative effects caused by the snow albedo reductions shown in (a-c).Calculations use the surface downward solar radiation and cloud cover fraction from the MERRA-2 reanalysis fields (see text and TableS2in the Supplement for details).

Table 3 .
Regional and seasonal mean BC-induced all-sky surface radiative effects (W m −2 ) for fresh snow over the Tibetan Plateau during 2000-2015.See TableS4in the Supplement for results of aged snow.