The snow broadband albedo is defined as (Aoki et al., 2011) α = ∫ λ 1 λ 2 r ( λ ) F ( λ ) d λ ∫ λ 1 λ 2 F ( λ ) d λ (1) where r ( λ ) is the spectral snow albedo, F ( λ ) is the incident spectral solar flux at the snow surface. This work is aimed at the parametrization of snow albedo in relatively clean regions such as Arctic and Antarctica. The function F ( λ ) is determined by the solar spectral irradiance at the top of atmosphere, atmospheric transmittance and solar elevation. In this work we assume that the solar zenith angle is 60° in the calculation of the spectral dependence F ( λ ) . The function F ( λ ) depends on the solar zenith angle. However, this dependence only weakly influences the BBA calculations (Grenfell and Perovich, 2008) because this function appears both in the dominator and nominator of Eq. 1. Also we use the parameterization of the smoothed function F ( λ ) proposed by Kokhanovsky et al. (2020): F ( λ ) = f 0 + f 1 ⁡ exp ( − ψ λ ) + f 2 exp ( − γ λ ) (2) where the corersponding parameters are given in Table 1. The multiplication of the function F ( λ ) by the spectrally independent parameter does not influence the calculations of BBA (Eq. 1).

The spectral albedo of clean plane-parallel snow surfaces is given by (Kokhanovsky et al., 2020): r ( λ ) = exp ( − k ( λ ) s ) (3) where k ( λ ) = 4π χ / λ is the bulk ice absorption coefficient, χ   is the imaginary part of ice refractive index ,   s = u 2 ( μ 0 ) l (4) is the effective attenuation scale (EAS), u ( μ 0 ) is the photon escape function, μ 0 is the cosine of the incidence angle. The parameter l is related to the effective grain diameter (EGD) d (Kokhanovsky et al., 2019): l = ζ d (5) where the shape factor ζ = 16 B / 9 ( 1 − g ) depends on the shape of ice grains. Here, g is the asymmetry parameter and B is the absorption enhancement factor (Kokhanovsky, 2006; Libois et al., 2014). Kokhanovsky (2006) has found that the shape factor is in the range 13–20 with the largest value corresponding to the case of spherical ice grains. We assume that ζ = 16 in this study. The value of   ζ can be considered as an additional parameter of the parameterization in terms of EGD, which accounts for the shape of particles.

The following expression for the escape function proposed by Kokhanovsky et al. (2021) is used: u ( μ 0 ) = 3 5 μ 0 + 1 + μ 0 3 (6) Let us substitute Eqs. 2, 3 into Eq. 1 and account for the fact that k ( λ ) s ≪ 1 in the UV and the visible. Then it follows for the UV and visible albedo: α = 1 − p s , (7) where p = [ ∫ λ 1 λ 2 k ( λ ) F ( λ ) d λ ∫ λ 1 λ 2 F ( λ ) d λ ] 2 (8) and [ λ 1 , λ 2 ] depend on the spectral region studied. The f unction   y ( λ ) = k ( λ ) can be approximated by the polynomial of the second order in the UV and the visible: y ( λ ) = a + b λ + c λ 2 (9) where the coefficients for various spectral intervals are given in Table 2. We have used the data for the imaginary part of ice refractive index obtained by Picard et al. (2016) (http://pp.ige-grenoble.fr/pageperso/picardgh/ice_absorption/). The substitution of Eqs. 2 and (9) in Eq. 8 makes it possible to derive the analytical expression for the parameter p. Namely, it follows: p = ( a + b 〈 λ 〉 + c 〈 λ 2 〉 ) 2 (10) where 〈 λ n 〉 = ∫ λ 1 λ 2 λ n F ( λ ) d λ ∫ λ 1 λ 2 F ( λ ) d λ (11) This integral can be evaluated analytically. In particular, one derives: 〈 λ 〉 = ∑ n = 0 2 f n X n ∑ n = 0 2 f n A n , 〈 λ 2 〉 = ∑ n = 0 2 f n Y n ∑ n = 0 2 f n A n (12) where A 0 = λ 2 − λ 1 , A 1 = Q ( λ 1 , ψ ) − Q ( λ 2 , ψ ) ,   A 2 = Q ( λ 1 , γ ) − Q ( λ 2 , γ ) (13) X 0 = ( λ 2 2 − λ 1 2 ) / 2 , X 1 = M ( λ 1 , ψ ) − M ( λ 2 , ψ ) ,   X 2 = M ( λ 1 , γ ) − M ( λ 2 , γ ) (14) Y 0 = ( λ 2 3 − λ 1 3 ) / 3 , Y 1 = N ( λ 1 , ψ ) − N ( λ 2 , ψ ) ,   Y 2 = N ( λ 1 , γ ) − N ( λ 2 , γ ) (15) Q ( λ , v ) = exp ( − v λ ) / v ,   M ( λ , v ) = ( 1 + v λ ) exp ( − v λ ) / v 2 , N ( λ , v ) = ( 1 + ( 1 + v λ ) 2 ) exp ( − v λ ) / v 3 . (16) The values 〈 λ 〉 , 〈 λ 2 〉 , and p calculated using equations given above for several spectral intervals [ λ 1 , λ 2 ] are given in Table 2.

Eq. 7 is valid at small values of the product η = ps. To extend the applicability of Eq. 7 with respect to the value of the parameter η, we propose to use the following parameterization of the UV and visible albedo: α U V , V I S = exp ( − p s ) (17) where the value of p is given by Eq. 10. Eq. 7 follows from Eq. 17 at η → 0. One can state that the UV and visible broadband albedos depend on just one parameter–the effective attenuation scale s. It appears that the same is true for the NIR and SW albedo of pure snow. In particular, we have found that the NIR and shortwave albedo can be parameterized as follows: α = α 0 + α 1 ⁡ exp ( − p s ) . (18)

Eq. 17 follows from Eq. 18 at α 0 = 0   and   α 1 = 1.   The values of α 0 , α 1 and p in Eq. 18 for various bands [ λ 1 , λ 2 ] were derived from the numerical evaluation of integrals present in Eq. 1 with account for Eqs. 2, 3. The numerical fitting of the derived dependence of the broadband albedo on the parameter s to the function shown in Eq. 18 using the software package ORIGIN has been used. We have used the data for the imaginary part of ice refractive index given by Picard et al. (2016) (in the visible and UV) and data of Warren and Brandt (2008) at longer wavelengths. The values of the derived parameters α 0 , α 1 and p are given in Table 3. In Table 3 and also in the text below we consider three spectral ranges: 0.3–0.7, 0.7–2.5 and 0.3–2.5 μm. The first interval incorporates UV and visible wavelengths and the second interval incorporates NIR wavelengths.

The dependence of the visible (0.3–0.7  μ m , which also includes UV part), NIR (0.7–2.5 μm) and shortwave (0.3–2.5 μm) albedo on the effective grain diameter calculated using analytical Eq. 18 with account for data in Table 3 and the numerical calculation using Eq. 1 at μ 0 = 0.65   are given in Figure 1. The results derived using the parametrizations proposed by Dang et al. (2015) are shown in Figure 1 as well. It follows approximately that BBA(SW)=(BBA(VIS)+BBA(NIR))/2.

We have found that the difference between numerical calculations using Eq. 1 and analytical Eq. 18 with account for data given in Table 3 is smaller than 1% for the shortwave and visible albedo and it is smaller than 2% for the NIR albedo at the diameters d > 0.1 mm characteristic for terrestrial snow covers. This difference is smaller that the respective error of the BBA measurement. The difference between the parameterization given by Eq. 18 and former parametrization given by Dang et al. (2015) is below 2% in the visible, 3% in the shortwave, and 6% in the NIR regions with our parametrization providing smaller albedos for a given size of ice grains. This is mostly due to the fact that the new compilation of the ice refractive index (Picard et al., 2016) gives larger values of the imaginary part of ice refractive index in the VIS/NIR part of the electromagnetic spectrum as compared to the corresponding values given by Warren and Brandt (2008). Also our model is based on the assumption that ice grains are irregularly shaped as compared to the parameterization for snow BBA based on the model of ice spheres proposed by Dang et al. (2015). The effective ice grain diameter used by us is defined as d = 3 v 2 σ , where v is the average volume of grains and σ is there average area of the ice grain projection on the surface perpendicular to the incoming light (Zege and Kokhanovsky, 2004). It follows for the spherical particles that σ = π d 2 / 4 and the effective diameter used by ice coincides with that used by Dang et al. (2015) in the case of spherical ice grains.

Eq. 18 can be used in the simplifaction of the corresponding blocks in Global Circulation Models (GCMs). It follows from Eq. 18 that the values of BBA of snowpacks with various microstructure coinside, if the effective attenuation scale s is the same. Also our new snow BBA albedo parametrization can be used for the determination of the snow BBA for given sizes of ice grains and illumination conditions and also for the solution of inverse problems of snow optics. In particular, it follows for the effective pure snow grain diameter from Eq. 18: d = l n 2 z ζ p u 2 ( μ 0 ) (19) where z = ( α − α 0 ) / α 1 . Also the snow specific surface area (SSA) can be derived from the BBA measurements. It is defined as S S A = K / ρ i d , where ρ i is the density of ice, K = 3 Σ 2 σ , Σ   is the average  surface area of ice grains .   The parameter K is equal to 6 for the spherical ice grains (and also for randomly oriented ice convex particles of the same shape) because it follows in this case: Σ = 4 σ . As a matter of fact, the value of K can be derived from independent measurements of the effective diameter and SSA for a given snowpack.

The dependence of the effective ice grain diameter on the SW BBA derived using Eq. 19 at u = 1 corresponding to the case of overcast sky (spherical or white sky BBA (Kokhanovsky et al., 2019)) is given in Figure 2, where the correspondence between SW BBA values and international snow classification (Fierz et al., 2009) is also presented.

We show the temporal behaviour of the SW BBA as measured by the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) network of pyranometers (Fausto et al., 2021) at the East GRIP (EGP) location (75.6N, 36 W) in Greenland and also temporal variation of the ice grain diameter derived from Eq. 19 in Figure 3. Also in situ ground measurements of the grain diameter as reported by Kokhanovsky et al. (2019) for July 8, 9, and 13 (2018) are given. One can see that the SW broadband albedo at EGP does not change considerably for the time interval studied. It is close to 0.8. The grain diameter is in the range 0.1–0.4 mm for most of cases. The diameters of grains derived from in situ measurements are close to those derived from Eq. 19.

The snow specific surface area derived from shortwave BBA measurements is given in Figure 4. It is mostly in the range 15–45 kg/ m 2 , which is consistent with the values of the SSA for the wind packed snow (Domine et al., 2008) to be expected at the EGP site located at 2.66 km above the sea level far from the ocean.

The intercomparison of the SSA determined from the SW BBA observations and those performed in the vicinity of the EGP station using NIR hemispherical snow reflectance observations under artificial light illumination conditions (Gallet et al., 2009) is given in Figure 5. Further details on the measurements of the SSA at the EGP station using NIR observations are given by Kokhanovsky et al. (2019). It follows that SW BBA and NIR measurements provide similar values of the SSA. This is due to the fact that both BBA and spectral snow reflectance measurements provide the same quantity–the effective snow grain diameter, which is used to derive the snow specific surface area. The difference in the measurements is due to the local variation of the SSA as shown in Figure 5 and also due to different average light penetration depths for the NIR (at 1.31 μm) hemispherical reflectance and SW BBA measurements.

The aveage values of shortwave broadband albedo measured at EGP site for 2016–2018 and average values of the ice grain diameter and snow specific surface area (derived using daily SW BBA measurements) are given in Table 4. It appears that the interannual variations are quite small at the site.

We compare the average values of the SSA for the 2 months (June, July) derived at the EGP site from BBA measurements for several years and also similar results from NIR measurements in Table 5. One can see that both datasets produce similar results for the average values of BBA. Clearly, the determination of the SSA from routine BBA measurements requires a fraction of time as compared to the hemispherical snow reflectance measurements under the artificial light (a laser diode) illumination conditions.

Let us consider the polluted snow now. Then it follows for the spectral albedo (Kokhanovsky et al., 2021): r ( λ ) = exp ( − ( k ( λ ) + G ( λ λ 0 ) − x ) s     ) (20) where G = k p ( λ 0 ) B c (21) c = c p c i is the relative concentration of pollutants in a snow layer, c p   ( c i ) is the volumetric concentration of pollutants  ( ice grains ) ,   k p ( λ 0 )     is the volumetic pollutant absorption coefficient at the wavelength   λ 0 , x is the absorption Angstrom exponent, B is the snow absorption enhancement coefficient (Kokhanovsky and Zege, 2004; Libois et al., 2014). It has been assumed that the impurity volumetric absorption coefficient k p ( λ ) can be presented as (Kokhanovsky et al., 2021): k p ( λ ) = k p ( λ 0 ) ( λ / λ 0 ) − x (22) and light scattering effects by pollutants are ignored as compared to light scattering by ice grains. Eq. 3 follows from Eq. 20 at c = 0 as it should be.

One can derive using Eqs. 1, 20 for the shortwave BBA in the spectral range 0.3–2.5 μm: α = α v i s p o l + Q α n i r 1 + Q (23) where α v i s p o l is the visible BBA (0.3–0.7 μm) for the polluted snow, α n i r is the NIR BBA (0.7–2.5 μm) of the pure snow and Q = ∫ 0.7 μ m 2.5 μ m F ( λ ) d λ ∫ 0.3 μ m 0.7 μ m F ( λ ) d λ     (24) The parameter Q gives the ratio of incident light flux in the band 0.3–0.7 μm to that at the NIR band 0.7–2.5 μm. It appears that Q = 1.08 and, therefore, one can assume that the SW albedo is approximately equal to the average of VIS and NIR albedos (similar to the case of pure snow SW BBA discussed above). The NIR BBA α n i r in Eq. 23 is given by Eq. 18 under assumption that NIR BBA is not influenced by the pollutants. This is often the case (Kokhanovsky et al., 2020). Therefore, the problem is reduced to the parameterization of α v i s p o l (Eq. 23).

It follows from Eq. 20 for the UV and visible albedo of polluted snow: r = 1 − ( k ( λ ) + G ( λ / λ 0 ) − x ) s (25) assuming that b = ( k ( λ ) + G λ − x ) s → 0 , which is a reasonable assumption in the visible and UV. Therefore, one derives for the visible BBA of polluted snow using Eqs. 1, 25: α v i s p o l = 1 − W s , (26) where W = [ ∫ λ 1 λ 2 k ( λ ) + G ( λ / λ 0 ) − x F ( λ ) d λ ∫ λ 1 λ 2 F ( λ ) d λ ] 2 (27) The numerical simulations show that the parameter W can be parametried as follows in the spectral range 0.3–0.7 μm: W = p v i s + m G exp ( γ x ) (28) where the corresponding constants m and are γ given in Table 6. The constant   p v i s coinsides with constant p for the visible range given in Table 3. This ensures that Eq. 26 coincides with Eq. 7 at c = 0 (pure snow case).

To extend the area of applicability of Eq. 26 we shall use an exponential approximation: α p v i s p o l = exp ( − ( p + q ) s ) (29) where q = m G e x p ( γ x ) . (30)

Eq. 29 is the final parametrization of the polluted snow visible broadband albedo. One can see that the visible BBA depends on the effective absorption scale s (d,   μ 0 ), absorption Angström parameter x and also on the pollution impact parameter G (dependent on the relative concentration of pollutants c, the volumetric absorption coefficient of pollutants at the selected wavelength k 0 ≡ k p ( 1 μm ) ) and the solar light absorption enhancement factor B ) . Eq. 17 follows from Eq. 29 at c = 0 (pure snow).

The shortwave snow albedo of polluted snow can be derived using analytical Eqs. 18, 23, 29. Namely, it follows: α = exp ( − ( p v i s + m G e γ x ) s ) + Q ( α 0 n i r + α 1 n i r ⁡ exp ( − p n i r s ) ) 1 + Q , (31) where α 0 n i r = 0.2335 , α 1 n i r = 0.56 , p n i r = 32.7 m − 1 (Table 3).

The accuracy of this approximation as compared to the numerical integration using Eq. 1 is shown in Figure 6 for snow contaminated by soot aerosol and the same solar zenith angle as in Figure 1. It has been assumed in calculations that (Kokhanovsky and Zege, 2004): B = 1.8, x = 1, and k 0 = 10 μ m − 1 . The results similar to those reported in Figure 6 except for the dust–loaded snow are given in Figure 7. It has been assumed that B = 1.8, the solar zenith angle equal to 27°, and values of the pair (x, G) coincide with those derived by Kokhanovsky et al. (2021) for the alpine snow polluted by the Saharan dust. They are given in Table 7 together with the derived value of the effective ice grain diameter for each case. It follows from Figure 7 that the proposed parameterizations can be used to derive shortwave snow broadband albedo with a high accuracy.

The spectral albedo measurements in the range 400–900 nm as reported by Kokhanovsky et al. (2021) can be used to estimate the shortwave snow albedo. This is possible because the parameters q and s (Eq. 29) can be assessed from spectral measurements. The SW BBA for the polluted snow cases investigated by Kokhanovsky et al. (2021) is presented in Table 7. The corresponding values were derived using Eq. 31 and data presented in Table 7.