Effect of small-scale snow surface roughness on snow albedo and reflectance

The primary goal of this paper is to present a model of snow surface albedo accounting for small-scale surface roughness effects. The model is based on photon recollision probability and it can be combined with existing bulk volume 15 albedo models, such as TARTES. The model is fed with in situ measurements of surface roughness from plate profile and laser scanner data, and it is evaluated by comparing the computed albedos with observations. It provides closer results to empirical values than volume scattering based albedo simulations alone. The impact of surface roughness on albedo increases with the progress of the melting season and is larger for larger solar zenith angles. In absolute terms, small-scale surface roughness can decrease the total albedo by up to about 0.1. As regards the bidirectional reflectance factor (BRF), it is found that surface 20 roughness increases backward scattering especially for large solar zenith angle values.


Introduction
The global energy budget is affected by surface albedo, which describes the level of brightness of the surface. Due to its central role for climate, it has been defined as an essential climate variable (ECV) by the Implementation Plan for Global Observing System for Climate in support of the United Nations Framework Convention on Climate Change (GCOS Secretariat 2006). 25 The large areal coverage of seasonal snow, together with the high reflectivity of snow, contributes to the relevance of snow albedo on the global energy budget (Flanner et al., 2011;Mialon et al., 2005). Snow component is also important for the liveability of dry and cold areas for both humans and ecosystems by providing a source of melt water in spring and shelter and insulation in winter. Changes in the duration of snow cover and snow type are vital for people and ecology of these areas.
Accurate large-scale monitoring of snow properties over large areas is only feasible in practice using satellite data-based 30 methods. Prior to that, it is required to obtain a detailed understanding of the reflectivity and scattering properties of snow.

Grain size and density profiles of snowpack
Measurements of snow depth, total density, water equivalent (SWE), humidity profile, temperature profile, grain size profile, surface roughness and surface impurity content were carried out at snow pits located in Sodankylä in an area with corner coordinates (67. 36°N, 26.63°E; 67.45°N, 26.86°E) in and in March 2010(Manninen and Roujean, 2014. 100 In addition crystal size photos of the snow layers, surface roughness photos and photos of the top surface impurities were taken. In this study we concentrate on the values measured in 2009. The air temperature in March was mostly below 0°C, whereas in April it was almost all the time above 0°C (Table 1). Hence, April represents the melting season and March is still midwinter. This is also clear from the increase in median density of the snowpack and the decrease of median snow water equivalent value from March to April. About 40 snow pit measurement points were located in a larger area (67. 42°N, 26.04°E; 105 67. 85°N, 26.91°E), where the maximum measured snow depth was 92 cm in March and 76 cm in April. The total density varied in the range 0.18 … 0.32 g/cm 3 in March and in the range 0.27… 0.57 g/cm 3 in April. The corresponding variation ranges for the snow water equivalent were 20 … 250 mm and 34 … 239 mm, but in April there was plain water in several places in the snowpack. Hence, the area covered by the snow pit measurements represents the local variation of snow properties to a large extent. 110 The traditional snow grain size (the largest dimension of the snow grains, Fierz et al, 2009) was visually estimated using graded plates, collecting the snow crystals from 10 cm-thick snow layers from the bottom to the top of the snowpack. For each analysed sample of snow crystals, in addition to the typical value of the largest grain dimension also its minimum and maximum value were provided. The measured snow grain sizes differ from the optically equivalent snow grain size (Mätzler, 1997;Neshyba 115 et al., 2003). To partly compensate for this, the minima of the largest grain diameter were applied in the radiative transfer calculations as the effective diameter. Although this causes some uncertainty in the interpretation of the computed absolute albedo values (particularly for the cases of fresh snow) it has much less impact on the derived effect of small-scale surface roughness on snow albedo. The density profile of the snowpack was measured for the same layer structure using the snow fork (Toikka, 1992). The variation range of the grain size and density is shown for the surface layer and for the whole snowpack in 120 Table 2.

Surface roughness from plate measurements
The surface roughness up to 1 m scale was measured in March and April, 2009 and in March 2010 by taking photos of a graded plate placed perpendicularly in the snowpack (Figure 2). The snow surface profiles were automatically calculated from the 125 photos using an image processing technique and the scale at the edge of the plate (Manninen et al., 2012). Control points at the scales were used both for the removal of the barrel distortion of the camera optics and transformation of the pixel coordinates to millimetres with photogrammetric methods. The plate surface roughness measurements were carried out at the https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. same sites as the snow pit measurements (Section 2.2). At each site, profiles were measured in two perpendicular directions with a 1 m interval along 50 m to 100 m distance. 130 The surface profiles were used to derive the rms height and correlation and their distance dependence (Manninen, 2003;Anttila et al., 2014). The snow surface roughness is close to a Brownian fractal surface  so that the logarithm of the rms height  depends linearly on the logarithm of the distance x and the correlation length L is linearly related to the distance 135 log = + log (1) where a, b, k0 and k are constants and k0 = 0 for an ideal Brownian surface. For each profile the values of the constants were 140 calculated by linear regression using varying sliding window sizes .
In addition, the rms slopes  were calculated for the measured spatial resolution, which was on the average 0.26 mm. The vertical resolution was about 0.1 mm and the horizontal resolution 0.04 mm (Manninen et al., 2012).

Surface roughness from laser scanning
In addition to the plate measurements, laser scanning data for snow roughness was utilized. The laser scanning data used in this study has been acquired using the FGI ROAMER system (Kukko et al. 2007). The system, including a FARO Photon 120 laser scanner, a NovAtel SPAN GPS-IMU system, and data synchronizing and recording devices, was mounted on a sledge, which was towed by a snow mobile. The data acquisition covers a 2.5 kilometres zone at each side of an official snow mobile 150 track (see Kukko et al. 2013 for more details). The landscape covered sparse pine forests and open bogs. The absolute precision of these measurements was analysed by Kaasalainen et al. (2011) to be better than 5 cm, while the relative accuracy required to observe the changes in the snow roughness were found to be 0.7 mm-2 mm for a static system, and less than 10 mm when the snowmobile was moving. The best repeatability was achieved at ranges closer to the scanning system, i.e., below 5 m. The data quality and precision were controlled using control points measured with a VRS-GPS (Leica SR530 receiver + AT502 155 antenna) precision GPS.
The laser profiles (about 16 profiles per 1m at 3 m/s snow mobile velocity) measured on March 18, 2010 were used to analyse the variation of the slope angles in a larger area than was possible using the plate profiles. The dependence of slope angles on distance between successive points was studied by calculating the root mean square (rms) slopes per size of horizontal distance 160 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
increment within the about 2.4 km long and 2 x 3.2 m wide scanned area. The scanned profiles on both sides of the snowmobile route were used.

Snow impurity content
The snow impurity was measured by filtering a melted sample of snow. The Quartz filters were analysed using NIOSH 5040 165 protocol. The increase of the median amount of impurities from March to April is obvious from Table 1. (Meinander et al., 2013;Meinander et al., 2014;Meinander et al. 2020).

Surface albedo
The surface albedo was operationally measured at Sodankylä (67. 36664°N, 26.628253°E downward, 67.36695°N, 26.62973°E 170 reflected) with a one minute interval using Kipp & Zonen CM11 Pyranometers. The site is surrounded by trees and houses, so that shadowing takes place in certain azimuth directions, when the solar elevation is very low. Hence, the measured white-sky albedo values are considered more reliable than the blue-sky values. The least shadowed azimuth direction in early March corresponded to the solar zenith angle value of 73° in the afternoon. Thus, the blue-sky albedo values used in the analysis were all taken from the afternoon, when the solar zenith angle equalled 73°. This means that the azimuth direction used increased a 175 bit during the spring, but it did not cause any additional shadowing problem. Yet, the clear-sky albedo of March 12 was replaced with the diffuse albedo dominating that day, because the clear-sky albedo value of a narrow time window seemed unrealistically small. Albedo values measured at the NorSEN mast on April 21, 2009 using a portable Kipp & Zonen CM 14 albedometer were used to calibrate the operationally measured albedo data in order to correct for the slight difference in location of the upward and downward looking pyranometer used for operational measurements. 180 The portable albedometer was used in April 2009 to measure the snow surface albedo in the same areas where the snow pits were made (Section 2.2). The instrument was installed on a short boom affixed on a lightweight camera tripod for easy transport. The tripod legs affect somewhat the reflected radiation measurements, and therefore a first-order correction was applied. It was based on estimation of the solid angle blocked by the tripod legs from fisheye lens photograph with a camera 185 mounted onto the albedometer position on the tripod, assuming a constant albedo of 0.1 for the dark carbon fibre legs. The albedometer was calibrated against a reference pyranometer at FMI-Helsinki prior to each campaign. The albedometer was carefully levelled on the tripod before measurements at each location and the stability of levelling monitored regularly, as melting snow may become unstable during the day.

Spectral reflectance
The spectral reflectance of snow was measured using the ASD FieldSpec Pro JR spectrometer on several days, specifically in the perfectly overcast conditions on March 13 and on the perfectly clear-sky day April 22, during the campaign in 2009. The irradiance spectra were measured as well. Every spectrum is an average of 30 individual spectra. The spectrometer was calibrated by the manufacturer prior to the campaigns. The instrument was powered on at least 15-20 minutes before each 195 measurement to ensure an even operating temperature. A Spectralon panel was used as white reference for the reflectance measurements. Most measurements took place from a height of 50-60 cm with an 8-degree foreoptic. The spectrometer was optimized before each measurement.

BRF 200
The bidirectional reflectance factor BRF of snow was measured using the Finnish Geodetic Institute's Field Goniospectrometer FIGIFIGO (Peltoniemi et al., 2005(Peltoniemi et al., , 2015(Peltoniemi et al., , 2014. FIGIFIGO consists of a motorized arm of length of 2 m, moving the optics head +-90 degree around nadir, and an ASD Field Spec PRO FR spectrometer recording the spectrum in the range of 350 -2400 nm. The azimuth is turned manually, and all angles and coordinates are recorded automatically, based on inclination, direction, and position sensors. The footprint is around 10 cm in diameter. FIGIFIGO gives spectrally resolved BRF data, 205 relative to Spectralon reference standard, from which also spectral albedo can be evaluated, but as the system is not absolutely calibrated in the field setup, real broadband albedos and BRF need external solar spectrum. In the results shown, a mean solar spectrum is used that may differ several per cent from the real time one.
In Mantovaaranaapa, 3 sets of rough snow were measured, and one set of smoother snow formed by a thin layer of windblown 210 grains. Another set of thin and rough snow was measured in Korppiaapa, but were not used in the present study. The sunlight measurements were complemented by set or artificial light measurements of smoother snow near the NorSEN mast.

Smooth snowpack albedo modelling 215
The TARTES model (available at: https://snowtartes.pythonanywhere.com/) was used to estimate the snowpack white-sky and black-sky albedo values (Warren, 1984;Warren and Brandt, 2008;Kokhanovsky and Zege, 2004;Baldridge et al., 2009;Libois et al., 2013;Libois et al., 2016;Picard et al., 2016). It is a fast and easy-to-use optical radiative transfer model and represents the snowpack as a stack of horizontal homogeneous layers. Each layer is characterized by the snow grain size, snow density, impurities amount and type, and two parameters for the geometric grain shape: the asymmetry factor and the 220 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. absorption enhancement parameter. The albedo of the bottom interface can be prescribed (here 0.13), although the bottom interface only markedly impacts thin snowpacks (<5 cm depth). The model is based on the Kokhanovsky and Zege (2004) formalism. The required input values for the model (density and grain size profile) were provided by the snow pit measurements of the SNORTEX campaign (Manninen and Roujean, 2014; Section 2.2). The amount of impurities was temporally interpolated from the values of measured days. The black-sky albedo values of the bulk snowpack were derived by 225 weighting the spectral albedo with the standard top of atmosphere (TOA) spectrum ASTMG173. The white-sky albedo values of the bulk snowpack were derived by weighting the spectral albedo with measured diffuse irradiance spectra of the cloudy day March 13, 2009.
The black-sky albedo values were calculated for three local incidence angle values of each plate surface roughness profile: the 230 mean, mean minus one standard deviation and the mean plus one standard deviation of the individual local incidence angle values determined for each slope of the surface roughness profiles. The nominal incidence angle was set to the solar zenith angle value at the time of the measurements of the plate surface roughness profiles and the density and grain size values of the snowpack layers. The blue-sky albedo values were obtained from the black-sky and white-sky albedo values using the fraction of diffuse irradiance operationally measured at Tähtelä. 235

Rough snowpack albedo modelling
From the theoretical point of view there is a difference in scattering from a snowpack having an ideally planar surface and a rough surface, because the rough surface may have an incidence angle distribution that markedly differs from the Gaussian distribution of incidence angles produced by a random volume of spherical scatterers partly shading each other. In addition, 240 the roughness may cause a markedly higher amount of multiple scattering thus reducing the amount of radiation escaping the target.
Scattering from randomly rough continuous surfaces is related to the characteristic size of the surface roughness with respect to the wavelength used (Beckmann and Spizzicchino, 1963;Ulaby et al, 1982;Tsang et al., 1985, Fung, 1994. When the 245 surface roughness of a randomly rough continuous surface is large compared to the wavelength of the electromagnetic wave, the scattering of the wave from the surface can be approximated by scattering from random facets (i.e. using the Kirchhoff approximation), whose slopes determine the scattering directions. As the shortwave illumination covers the wavelength range of about 300 nm -2500 nm, all structures in the mm scale (or above) are large compared to the wavelength, so that a facetbased surface scattering calculation is reasonable. Each facet is taken to represent a volume of random scatterers and the local 250 incidence angle of the incoming radiation is the angle between the normal of the facet and the solar zenith angle (Figure 3).
The surface of a snowpack is not a continuous solid surface, but when the snowpack surface is rough, the incidence angle distributions of the scattering elements may deviate from that of a planar surface with randomly oriented scatterers. In addition, https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
it is possible that a photon escaping one facet hits another facet. The snowpack scattering can then be thought to have elements both of bulk volume scattering and surface scattering. The following 2D analysis demonstrates this idea. 255 Multiple scattering between facets can be taken into account using the photon recollision probability theory (Knyazikhin et al., 1998;Panferov et al., 2001;Smolander and Stenberg, 2005;Rautiainen and Stenberg, 2005;Stenberg et al., 2008;Stenberg and Manninen, 2015;Stenberg et al., 2016). The formulation is shown in Appendix A separately for diffuse and direct irradiance. The essential equations are repeated here. Firstly, for the diffuse case, the white-sky albedo w (Lucht et al., 2000, 260 Schaepman-Strub et al., 2006 is related to the average number of facet-to-facet scattering events n where w0 is the white-sky albedo of the bulk volume.

265
Second, for direct illumination, the black-sky albedo b(i) (Lucht et al., 2000, Schaepman-Strub et al., 2006 is approximately related to the w0, n and the average number m(i) of facet-to-facet scattering events in direct illumination conditions by where b0(i) is the black-sky albedo of the bulk part of the snowpack (i.e. the albedo without the surface roughness contribution). In this study the bulk albedo values w0 are produced using the TARTES model (Section 3.1).
The albedo  in mixed illumination conditions is typically estimated using the weighted mean approximation of the two 275 extreme values w and b (Lucht et al., 2000, Schaepman-Strub et al., 2006Román et al., 2010) where f is the fraction of diffuse irradiance. According to Eqs. 3 -5 the blue-sky albedo is then estimated from 280 Obviously, surface roughness decreases the white-sky albedo and typically also the black-sky and blue-sky albedo. Only when m > n, surface roughness can increase the black-sky albedo of bright targets. The effect of surface roughness is non-negligible 285 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. even when the roughness is not large. On the other hand, the larger the roughness is (i.e. the larger n and m are), the larger the effect is for darker targets. Hence, for snow the effect is larger in the near-infrared than in the visible wavelengths and in midwinter roughness alters the broadband albedo only slightly, whereas during the melting season, when the snow is darker, the effect of the roughness may be much larger. Thus, the effect of surface roughness would be larger for old snow than for new snow, which may explain part of the quick darkening of snow during the melting season. 290

Ray tracing analysis of surface roughness
Scattering from the snow profiles measured using the plate was analysed by a ray tracing method using 1000 equally spaced rays per profile per direction. The number of hits ns on the surface (unity for single reflection and larger for multiple surface scattering) and the direction of the escaping reflected ray was calculated as a function of the zenith angle of the incoming ray 295 with an interval of two degrees from 0° to 80°. Scattering from the surfaces was calculated assuming completely smooth facets of continuous material, so that scattering from a facet was determined as mirror reflection. For each angle the case was calculated separately for rays coming from the left and from the right and the two results were unified to improve statistics. This choice was motivated by the known fact that, even when the snow is highly forward scattering, it is with respect to the direction of the incoming solar radiation, not the wind direction. Thus, the dominant scattering direction moves with the sun 300 during the day. In some cases, the ray was trapped to infinite reflection from facet to facet. But these quite rare (< 1 %) cases were excluded from the statistical analysis.
Late in spring the scattering from a snowpack often also contains a component that is something between volume and surface scattering, namely deep narrow pits generated by impurities that have sunk downwards in the snowpack due to melting caused 305 by absorption of solar radiation. This kind of an effect is not easy to take into account either in volume scattering or surface scattering, because they don't affect the roughness or density in a random way. Their contribution to albedo is beyond the scope of this study.

Surface roughness and inputs for albedo modelling
To start with, we consider the snow surface roughness from the plate measurements. The rms slope angle calculated from the plate roughness measurements had an increasing trend from March to April ( Figure 4). The surface height distributions developed towards a more Gaussian distribution from March (R 2 = 0.97) to April (R 2 = 0.99). This was evident also from the change of the mean skewness value that was 0.4 in March and -0.12 in April. 315 The mean number of individual reflections ns per ray on the surface has an increasing trend during the melting season ( Figure   5). It correlates well (R 2 = 0.83) with the rms slope angle () of the profile. A good general fit to all measured plate profile data was 320 = 1 + 0.355332 cos 4 − 1.08275(1 − exp(1.75 cos 1/4 4 )) .
The mean zenith angle to which the radiation escapes from the surface (when mirror reflection from the facet is assumed) is even more strongly correlated (R 2 = 0.93) with  ( Figure 6) All angles ( i, o) are given in radians in the above equations, and o > 0 (o < 0) for forward-scattering (backward scattering). Obviously, the probability for backward scattering increases with increasing incidence angle and increasing .
Since  increases with time, also the probability for stronger backscattering from the snow cover increases with time. Indeed,

330
it is well known that older snow cover is less strongly forward scattering than new midwinter snow (Peltoniemi et al., 2010).
Unfortunately, the  values depend strongly on the scale of the measurements. The laser scanning data is well suited to demonstrate this, because the horizontal distance between successive data points increases from the beginning to the end of the scan line. On March 18, 2010 plate profile measurements and laser scanning were carried out in the same relatively flat 335 wetland area Mantovaaranaapa (67.4°N, 26.7°E). The laser scanning data covered a 2.4 km long and about 3.2 m wide area on each side of the snowmobile route. Altogether 10 plate profile measurements were taken directly after the scanning, at about 100-200 m intervals starting from the western edge of the scan route (Kukko et al. 2013). The average rms slope angle of the 10 plate profiles was 30.7° ( = 0.54) with an 80 % variation range of 24.7°… 34.3° ( = 0.43 … 0.60). Consequently, ns would then vary in the range 1…1.5 ( Figure 5). The rms slopes were calculated also from the laser scanning profiles as a 340 function of horizontal increments, which were within the range 5 mm … 100 mm. The number of points per distance varied in the range 4 thousand … 2.2 million. Nonlinear regression to the 36 points in the range 5 mm … 100 mm produced an exponential curve that approaches the mean rms slope value obtained from the 10 plate profiles of the same area and ( Figure   7). Shorter increments of the laser data could not be reliably used in the analysis.

345
As the laser scanner covered a much larger area than the plate profiles, that data gives an estimate of the rms slope variation in a larger area and is based on a larger number of individual slope values. For the laser data set, the median difference between the 90 % quantile curve of the slope angle values and the rms slope angle value curve was 6.1°. The corresponding median https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. difference for the 10 % quantile curve from the rms slope angle curve was -5.9°. For the plate profiles, 90% and 10 % quantile values of the rms slope angle differed from the mean rms slope value by 3.6° and -6.0°. Thus, the larger area covered by the 350 laser scanner shows a larger variation of the rms slope angle, as one could expect. Obviously, the laser scanner data can be extrapolated to estimate  at higher horizontal resolution than the measurements directly enable, but then  has to be analysed as a function of the horizontal distance increment (Figure 7). However, the strong variation of  with spatial resolution suggests that using less scale-dependent surface roughness descriptors would be desirable, if they just can provide the information needed. 355 The relationship between other surface roughness parameters (such as rms height  and correlation length L) and  is in general not strong even for a Gaussian surface height distribution (Beckmann & Spizzicchino, 1963). For the whole period (March 3 -April 28, 2009) the ratio of  /L (determined for 60 cm distance) correlated however relatively well with  the R 2 values being 0.70, 0.62 and 0.67 for the whole data range, March and April, respectively, but the best descriptor of  was found to be 360  /b (Eq. 1) its R 2 values for the linear correlation being 0.78, 0.68 and 0.82 for the whole data range, March and April, respectively.  tends to increase with the progress of the melting season (0.002 radians per day). Likewise, its correlation with  /b increases during the melting season. It was therefore examined, whether the measured surface albedo correlates well with the measured surface roughness parameters. Using just the rms height (derived for a 60 cm horizontal scale) as an explanatory variable of the albedo the coefficient of determination was R 2 = 0.81. The relationship between the albedo and surface 365 roughness parameters that are scale-independent in a large range (Manninen, 2003)

Snow albedo spectra: measured vs. modelled
Two examples of snow nadir reflectance spectra measured with the ASD spectroradiometer are shown in Figure 9. On March 13, the sky was completely overcast, whereas on April 22 it was perfectly clear. For comparison, corresponding albedo spectra modelled using TARTES are shown. The ASD reflectance spectra were scaled so that the derived broadband reflectance value 375 matched the calibrated operationally measured broadband albedo value. No BRF was available for the clear-sky case for the location of ASD, but for old rough snow in the area it was relatively flat (see Section 4.5), so that the comparison of the spectral reflectance and albedo values seems reasonable enough to enable the choice of the grain shape to be used in the TARTES model calculations. The spectra modelled with TARTES accounted for the empirical grain size and density values, as well as for black carbon but not for organic carbon (Table 1), which included needles and various tree trash deposited on snow. In March, the impurity content in surface snow was very low, while in April, it was roughly 3 times higher (Table 1). In March a better fit is obtained using fractal grains, which result is also supported by photos taken of snow grains and the fact that the snow was fresh. In April the modelled albedo favoured the use of spherical grains rather than fractals in the calculations. The photos taken about the snow grains also supported the use of spherical grains in April. Thus, the TARTES results seem to represent the snowpack in March well, but in April when the melting has been going on for a longer time the modelled albedo 385 values are higher than the empirical reflectance values both in the visible and NIR wavelengths (less than 1 µm), which dominate the value of the broadband albedo of snow.

Broadband albedo: measured vs. modelled
The evolution of the operationally measured broadband albedo in Tähtelä is shown for the periods March 12-19 and April 21-390 28 in Figure 10 together with the corresponding albedo values simulated using the TARTES model and the grain size and density measurements of the day in the Sodankylä area. Following the justifications outlined above, fractal grain shapes were used in March and spheres in April. The simulated values tend to exceed the measured ones, especially in April. Besides, the variation range of simulated albedo values is rather small compared to that of the empirical broadband reflectance values.
However, it will be demonstrated next that taking into account the surface roughness decreases the difference between 395 simulated and empirical albedo estimates.
The albedo model taking into account both volume and surface scattering (Eqs. 3 -5) was applied so that w0 was the value provided by the TARTES model based on the measured density and grain size profile and impurity content. The values for n and m were derived from the empirical values of ns. Namely, m = ns -1 for the local solar zenith (il) angle range derived using 400 the ray tracing method. And n is the weighted mean of ns -1, where the weights are cos(il)sin(il). The results are shown in Figure 11 and Figure 12.
First, the ratio w/w0 of the total white-sky albedo and the bulk white-sky albedo provided by the TARTES model is considered in Figure 11. In March, this ratio varies mainly between 0.97 and 0.99, indicating that small-scale surface roughness decreases 405 the snow albedo typically by 1-3%. With the progress of snow melt, the effect of surface roughness increases markedly. On 26-27 April (Julian days 116-117), the median of w/w0 falls below 0.9, indicating an over 10% decrease in snow albedo. The relative difference between the total and bulk albedo values is about the same for the black-sky case as for the white-sky case, but the solar zenith angle naturally slightly complicates that relationship (see Eqs. 3 and 4).

410
The modelled albedo values are further compared to observations in Figure 12. The variation range of the simulations shown with the background shading is based on the variation of the grain size, density and local incidence angle of the measured https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
profiles. Overall, the inclusion of surface roughness improves significantly the agreement of modelled albedos with the observed ones. A notable overestimation remains, however, on Julian days 77-78, and will be discussed in Section 5. All in all, the model is robust enough to be reasonably applied to empirical data of grain size, density and surface roughness. 415 Obviously, taking into account the surface roughness contribution improves the match of empirical and modelled results, but still it is very clear that the grain size, grain shape and SWE (or density) are dominant parameters, because the amount of volume scattering also affects directly the amount of surface scattering. The variation range of the modelled albedo is much larger for clear-sky cases than diffuse cases, which is understandable as some of the variation in clear-sky cases comes from 420 the solar zenith angle variation during the day.

Albedo during snow metamorphosis on April 22
The Hence, the modelled albedo results were averaged to get one value per surface sample. The light grey band in Figure 13 shows 430 the variation range of the broadband-converted reflectance spectra measured with the ASD spectrometer in the same area at the same time. Obviously, the modelled profile albedos fit in that range. The mean of the empirical albedo values is 0.66 and the modelled mean values are 0.72 and 0.68 for volume scattering only and for both volume and surface scattering, respectively.
Taking into account the snow surface roughness thus improved the average modelled albedo estimate.

BRF
Since the contribution of surface roughness to the total albedo is markedly smaller than the contribution of the bulk volume scattering, it is clear that the volume scattering dominates also the BRF. However, the contribution of surface roughness is not negligible and the BRF of the surface scattering component may differ markedly from the bulk volume component, resulting in a complex total BRF. The ray tracing based surface BRFs (without any volume scattering contribution) were compared with 440 empirical BRFs measured using FIGIFIGO (Figure 14) in Mantovaaranaapa in April 22, 2009. The area is an aapa mire, which late in spring affects the snow properties markedly (Figure 15). Sporadically the snow had melted and refrozen. Since the surface roughness profiles produce only 2D information and scattering angles differ markedly in 2D and 3D, the calculated ray tracing based 2D BRFs were converted to 3D versions assuming that each facet has besides the measured vertical angle https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. also an azimuth angle obeying a random uniform distribution between 0 … 180°. In fact, calculations were made for the range 445 0 … 90°, assuming the case to be symmetrical with respect to azimuth angle, like in constructing the FIGIFIGO based BRFs.
The 3D conversions make the peaks of the 2D scattering angle distributions slightly less distinct. No atmospheric contribution was included in either data set.
The comparison of the ray tracing and surface roughness based BRFs with the empirical FIGIFIGO-based BRFs were made 450 in the principal plane, since the azimuth information of the former ones was just a statistical assumption to convert the 2D principal plane BRF to 3D principal plane BRF. The surface scattering BRFs are typically peaked to the direction  of forward scattering matching mirror reflection of the surface plane of the snowpack. In addition, a backward peak in the direction  -90° is also strong in most cases ( Figure 16). The balance between forward and backscattered intensity varies with incidence angle so that large incidence angles favour surface backscattering due to roughness ( Figure 17, Table 3). Although mere surface 455 scattering would lead to dominantly backward scattering in Mantovaaranaapa due to the large incidence angle values (55° … 64° for FIGIFIGO BRFs), the empirical BRFs measured using FIGIFIGO dominated by the volume scattering are still dominantly forward scattering. The balance between forward/backward volume scattering is related to the grain shape (Peltoniemi et al., 2010). But indeed, the smoothest snow sample produced least backscattering, the ratio of the backward to the forward scattered amount of radiation was 0.58, whereas the corresponding ratio was on the average 0.73 for the rough 460 BRFs. This result is quite in line with the general ray tracing analysis results that surface roughness increases the fraction scattered backwards (Table 3). Also, in a previous theoretical study of Gaussian surfaces it has been shown that roughness affects the maximum direction of backward scattering (Jämsä et al., 1993). For the profiles measured in Mantovaaranaapa on April 22 the ratio of the backward and forward scattered radiation amounts in the incidence angle range of the FIGIFIGO measurements varied from 1.0 to 1.46. One has to take into account that the plate profiles register roughness in 1 m scale, 465 whereas FIGIFIGO measures samples of 10 cm diameter. Hence the largest spatial roughness may not necessarily show up as strongly in the FIGIFIGO results.

Discussion
In this study, the equations combining the volume scattering and surface scattering were derived using the photon recollision 470 theory (Appendix A), because this theory could be extended to include the surface scattering effect. However, to describe properly the volume scattering of real snowpacks, it is essential to pay attention also to layer structure, grain shapes and various types of impurities etc., so that a more complex description is typically needed for realistic volume scattering estimation. In principle, the photon recollision theory could be extended to take the layer structure of the snowpack into account by just letting p to be a function of the depth, i.e. p = p(z), where z would be the distance from the bottom or top surface. However, 475 that is beyond the scope of this study. Hence, the surface scattering part is developed so that in principle it can be combined https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
with any volume scattering method. One just applies the estimates of w0 and b0 derived with the chosen volume scattering model in Eqs. 3 -5. Hence, to obtain more realistic volume scattering estimates for the snowpack, we used the TARTES model for volume scattering in the simulations.

480
The ray tracing analysis showed that the backward scattering increases with increasing surface roughness and increasing incidence angle of the illumination. However, that analysis concentrates only on surface scattering without any volume scattering contribution. Combining the surface and volume scattering contributions to BRF, perhaps using an adding procedure for radiative transfer, would be an interesting topic for future work.

485
Although the surface scattering model used here includes multiple scattering from the surface, it may be that a surface layer containing very deep pits would benefit from some special attention. Namely, surface roughness measurements methods are usually designed for typical roughness of about the same variation range horizontally and vertically, not for extremely deep pits. To some extent the pit structure will be taken into account by the volume scattering models, since they affect the density of the surface layer of the snowpack. However, their very anisotropic (vertical) orientation is not well described by random 490 scattering of a layer with reduced density. The pits act like illumination traps so that a larger part of illumination reaches lower layers of the snowpack before it is absorbed or scattered upwards. Therefore, the bulk albedo is smaller than for a completely random volume of the same density.
An example case is offered to illustrate this effect. Slight snow precipitation took place on March 16 so that a very fluffy 495 surface structure of large dendritic snow crystals was formed on the snowpack (Figure 18 and Figure 19). The rms slope based on the plate method showed an increase with time from about 0.38 to about 0.63 with R 2 = 0.68. The rms height and correlation length also manifested clear evolution during those days ( Figure 20). A related change is obvious also using the roughness parameters a and b ). Yet, the surface roughness measurements based on the plate method or laser scanning are not able to catch the deep pit structure of the surface, because of shadowing effects. Therefore, the simulated albedo is 500 higher than the empirical one ( Figure 12, Julian days 76-78). However, even if the 2D-surface roughness were characterized properly, the surface scattering model based on a statistical approach of random scatterers would not be ideal for a case with a distinct periodic surface structure. For example, one could analyse separately the percentage of illumination that will be completely trapped by the deep pits and reduce the simulated total albedo with that fraction.

505
The plate profile method has the advantage of high spatial resolution. Its main drawback is that it can be used only, when the snow is relatively soft. It has been successfully used in Finnish Lapland  and at Greenland Summit (Manninen et al., 2016), but Antarctic snow is typically so hard that it is not possible to immerse the plate in it. In addition, icy and crusty snow surface of Finnish Lapland in 2018 and 2019 turned out to be too hard for the plate. For laser scanning, however, the hardness of the snowpack does not cause any problems. Indeed, laser scanning shows great potential for measuring snow surface roughness as it can cover large areas with high point precision accuracy. It is a particularly good method for measuring larger-scale roughness from 5-10 cm upwards. The limiting factor in finer scale roughness measurements is the data resolution and footprint size of the laser beam. So far, the scanners with the highest point density, accuracy and smallest spot size are meant for indoors use, but as the technology improves, smaller and smaller features become measurable also outdoors. In addition, the fractal nature of snow surfaces enables extrapolation of surface roughness from cm 515 scale to mm scale (Kukko et al., 2013). Another benefit of using laser scanning for surface roughness measurements is that it leaves the surface intact. This enables repeatable measurements of the same surface giving a means to study the evolution of surfaces in time. The backscattering intensity of the laser beam is typically stored for each point measured by laser, and in the most modern scanners also the range deviation is stored. These features have so far not been widely used, but it could potentially be used in the future for surface scattering property measurements and snow surface classifications. 520

Conclusions
A method was developed to model the effect of surface roughness on albedo besides the volume scattering. It can be combined with any volume scattering model. Applying measured surface roughness values to the model produced results closer to measured values than only volume scattering simulations made with the TARTES model. The surface roughness is described 525 by the average number of surface scattering events per ray, which is currently estimated from the rms slope values of the measured surface roughness profiles. High empirical correlation (R 2 = 0.9) of albedo with just two surface roughness related parameters supports the importance of surface roughness to albedo.
The albedo modelling results taking into account also the surface roughness indicate that it may decrease the albedo by about 530 1-3 % in midwinter and even more than 10 % during late melting season. The effect is largest for low solar zenith angle values and lower bulk snow albedo values. Hence, the effect is larger early and late times of day everywhere and it increases during the melting season especially at high latitudes, where the sun elevation is lower. Increasing surface roughness also favours more backwards scattering. 535 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.

Scattering of diffuse radiation
The scattering of light in canopies has for several years successfully been described with spectral invariants and the so-called 540 photon recollision theory, p-theory, (Knyazikhin et al., 1998;Panferov et al., 2001;Smolander and Stenberg, 2005;Rautiainen and Stenberg, 2005;Stenberg et al., 2008;Stenberg and Manninen, 2015;Stenberg et al., 2016). The central parameter, the photon recollision probability p is spectrally invariant and depends on the amount of scattering surface in the volume. Canopies don't have distinct upper surfaces, hence the p-theory is developed so far only for a scattering volume, but it has already been successfully combined with forest floor scattering also for snow covered cases (Manninen and Stenberg, 2009). Here the p-545 theory is applied to snowpack scattering taking into account that the snowpack has a distinct surface, which may be rough.
A simple way to take into account the surface roughness effect on scattering is to consider every facet of the snowpack as a separate volume of scatterers. When the irradiance i0 = i0() first arrives at the facet, it enters a volume scattering sequence, which can be described with the spectrally invariant photon recollision probability p of the bulk part of the snowpack and the 550 single-scattering albedo of the snow grains  = () (Knyazikhin et al., 1998;Panferov et al., 2001). The radiation absorbed and scattered by the volume of the facet a0 and s0, respectively, are (Smolander and Stenberg, 2005;Rautiainen and Stenberg, 2005;Stenberg and Manninen, 2015) For simplicity the dependence of , i0, a0 and s0 on the wavelength  is not shown explicitly in the equations. The radiation escaping the volume of the facet either escapes altogether or hits another facet and experiences another volume scattering sequence. The probability of the latter case is defined to be the surface photon recollision probability ps. Because the snow grains in the bulk part are completely surrounded by other snow grains while in the surface almost only half of the surrounding volume may contain snow grains, it is essential to assume that p and ps are not identical. 560 The radiation escaping the snowpack altogether without hitting another facet, r0, is where q = q() is the fraction of the volume scattering escaping upwards. Essentially it corresponds to Q defined by Stenberg et al. (2016, Eq. 24), but as it does not contain the fraction scattered upwards by the surface, it is not the total upwards scattered fraction of light. Hence q is used here instead of Q. Theoretically the values for q are in the range 0 … 1 being the larger the 565 thicker and denser the scattering layer is. The radiation hitting another facet is https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
The absorbed (a1) and scattered (s1) amounts of radiation by the second volume scattering sequence are The amounts of radiation escaping (r1) and entering the following scattering sequence (i2)  The total amounts absorbed and scattered by the surface and volume are then (A18) And the total upwards escaping radiation is (A19) When the surface of snow is very smooth, the number of additional facet scattering sequencies (n) the photon has before it 590 escapes altogether may not be very large. When n = 0 and ps = 0, there is no facet-to-facet scattering, i.e. the case is the normal https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. volume scattering case. The total absorbed and upwards escaping radiation of the snowpack are derived as infinite geometrical sums and are The total white-sky (diffuse) albedo w of the snowpack is then (A22) The white-sky (diffuse) albedo of the bulk part of the snowpack (without facet-to-facet scattering) w0 is Hence, the total white-sky albedo is simply 600 (A24) Estimation of the parameter ps is not trivial, but when the average value for n is known from measurements, then a reasonable estimate can be obtained from the total scattered energy (Eq. A18), which can then also be calculated using a finite geometric sum of n terms with recollision probability of unity. This must then equal the probabilistic infinite sum formulation of Eq.
(A18). Taking into account Eq. (A23) the following relation is obtained 605 Now the value of ps can be solved and is The relationship between the total albedo and the bulk albedo is then When the surface does not cause additional scattering, n = 0 and the total albedo equals the bulk albedo. The larger the n is the smaller is w.
For n = 1, q equals the bulk white-sky albedo and for larger (smaller) values of n it is slightly smaller (larger) for medium albedo values.

Scattering of direct radiation 620
For the direct component of solar illumination one has to take into account that the irradiance depends besides the wavelength also on the solar zenith angle i, i.e. i0 = i0(, i). The photon recollision probability in the bulk snowpack will be denoted for the first scattering sequence p1 = p1(il) and for latter sequences p, assuming that the incidence angle dependence is lost after the first volume scattering sequence. Here the local incidence angle of the facet is denoted by il. Then, the absorbed radiation 625 of the volume will be (Stenberg and Manninen, 2015) and the radiation scattered by the volume will be When p1 = p the above formulas reduce to Eqs. A1 and A2 as they should. 630 The first surface scattering sequence shall be treated separately, because the volume scattering depends on the incident solar zenith angle. The surface photon recollision probability is denoted by ps1 = ps1(il) for the first facet-to-facet scattering sequence and by psr for the latter facet-to-facet scattering sequences. One should notice that psr is not necessarily equal to ps of the diffuse case, although they certainly approach each other asymptotically. The reason for not taking them to be identical immediately 635 is the typically small number of surface scattering events per ray of relatively smooth midwinter snow surfaces. The radiation escaping the facet altogether after just volume scattering, r0, is where q0 = q0(, i) denotes the fraction of the scattered radiation escaping upwards from the snowpack during the first volume scattering sequence. The radiation hitting another facet is 640 The further scattering sequencies are assumed to be independent of the solar zenith angle of the original incident radiation.
The absorbed (a1) and scattered (s1) amounts of radiation by the second volume scattering sequence are then The amounts of radiation escaping (r1) and entering the following scattering sequence (i2) are Formulas for corresponding radiation components in the following round are The components corresponding to the jth round are The total amounts absorbed and scattered by the surface and volume are then 660 And the total upwards escaping radiation is where m is the number of facet-to-facet scattering events. The total radiation escaping the snowpack upwards is derived again 665 as an infinite geometrical sum The total black-sky (direct) albedo b of the snowpack is then (A50) 670 The black-sky (directional) albedo of the bulk snowpack b0 is Taking into account also Eq. (A23) the relationship between the total black-sky albedo and the bulk snowpack black-sky and white-sky albedo is https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.
Estimating ps1 or psr (whether equal to ps or not) is not trivial, but like in the case of diffuse irradiance one can benefit from measured value of m by requiring that the total scattered energy (Eq. A47) is the same, whether it is calculated from the probabilistic infinite sum or a deterministic sum with recollision probability of unity, i.e.
Unfortunately, there are two variables, ps1 and psr, to solve, but only one equation. Hence, the theory does not provide an exact 680 solution for both ps1 and psr, but only psr remains explicitly in the equation of the black-sky albedo When the assumption that psr = ps is valid, the following formula is obtained for the black-sky albedo using Eq. (A26) When m = n, the black-sky albedo reduces to b0 w0 n . Further on, when b0 = w0, the total black-sky and white-sky albedo 685 values are equal as well, if m = n. The above approximation of the black-sky albedo is reasonable only, when the assumption psr = ps is good. Hence, further studies are needed for proper estimation of psr and ps1.
The number of facet-to-facet scattering rounds m is estimated to be .
When ps1 → ps, q0 → q and b0 → w0, m → n, as it should. It should be noticed that the number of facet-to-facet scattering 690 rounds in direct illumination (m) and in diffuse illumination (n) are not necessarily equal, the difference being related to the difference of b0 and w0 and q and q0. The estimate for q0 can now be derived from Eq. (A56).

Data availability 695
Some data of the SNORTEX campaign is directly listed in the report (Manninen and Roujean, 2014). The plate snow profile data and the ASD data are available upon request from Finnish Meteorological Institute. The FIGIFIGO measurements and the laser scanning data are available upon request from Finnish Geospatial Research Institute/National Land Survey.  Appl. Phys., 81, 1509Phys., 81, -1517Phys., 81, , 1997Phys., 81, . https://doi.org/10.5194/tc-2020 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.                 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. Figure 7. Relationship between rms slope and the horizontal distance between the points used for its determination according to the laser scanning data of March 18, 2010. For comparison the mean value and variation range of the rms slope values derived from the 10 plate profiles for horizontal resolution 0.25 mm in the same area in the same day are shown in red. The dashed black curve is the regression to the 36 laser scanning based points (black polyline) and the grey shaded area covers the 80% variation range of rms slope at the distance in question. 1125 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License. The corresponding variation range of simulated albedo values based on simultaneous grain size and density measurements and temporally interpolated impurity content in the same day are shown in grey. The minimum, mean and maximum values are indicated with darker grey curves. The vertical 'error bars' marked for some of the measured albedo values are based on the variation range of the measured ASD spectra based broadband reflectance values during the same day in the larger test area in Sodankylä. 1150 https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.  https://doi.org/10.5194/tc-2020-154 Preprint. Discussion started: 30 July 2020 c Author(s) 2020. CC BY 4.0 License.