Snow albedo sensitivity to macroscopic surface roughness using a new ray tracing model

Most models simulating snow albedo assume a flat and smooth surface, neglecting surface roughness. However, the presence of macroscopic roughness leads to a systematic decrease in albedo due to two effects: 1) photons are trapped in 10 concavities (multiple reflection effect) and, 2) when the sun is low, the roughness sides facing the sun experience an overall decrease in the local incident angle relative to a smooth surface, promoting higher absorption, whilst the other sides has weak contributions because of the increased incident angle or because they are shadowed (called the effective angle effect here). This paper aims to quantify the impact of surface roughness on albedo and to assess the respective role of these two effects, with 1) observations over varying amounts of surface roughness, and 2) simulations using the new Rough Surface Ray Tracer 15 (RSRT) model, based on a Monte Carlo method for photon transport calculation. The observations include spectral albedo (400-1050 nm) over manually-created roughness surfaces with multiple geometrical characteristics. Measurements highlight that even a low fraction of surface roughness features (7 % of the surface) causes an albedo decrease of 0.02 at 1000 nm when the solar zenith angle (Өs) is larger than 50°. For higher fractions (13 %, 27 % and 63 %), and when the roughness orientation is perpendicular to the sun, the decrease is of 0.03 – 0.04 at 700 nm and of 0.06 – 20 0.10 at 1000 nm. The impact is 20% lower when roughness orientation is parallel to the sun. The observations are subsequently compared to RSRT simulations. Accounting for surface roughness improves the model observation agreement by a factor two at 700 nm and 1000 nm (errors of 0.03 and 0.04, respectively), compared to simulations considering a flat smooth surface. The model is used to explore the albedo sensitivity to surface roughness with varying snow properties and illumination conditions. Both multiple reflections and the effective angle effect have more impact with low SSA (< 10 m kg). The effective angle 25 effect also increases rapidly with Өs at large Өs. This latter effect is larger when the overall slope of the surface is facing away the sun and with a roughness orientation perpendicular to the sun. For a typical alpine snowpack in clear sky conditions, a broadband albedo decrease of 0.05 causes an increase of the net short wave radiation of 80 % (from 15 W m to 27 W m). This paper highlights the necessity to consider surface roughness in the estimation of the surface energy budget. 30


Introduction
Spectral albedo quantifies the proportion of solar energy reflected by a surface for each wavelength, and governs the quantity of solar radiation absorbed in the snowpack. Because snow has an overall high albedo in the solar spectrum, a small decrease in albedo (e.g. from 0.85 to 0.75) drastically increases the proportion of absorbed energy (from 25% to 15%; Genthon, 1994).
Thus, a reduction in albedo has important consequences on the surface energy budget, impacting surface temperature (Mondet 35 and Fily, 1999, Picard et al. 2012, Fréville et al., 2014, and the hydrology of watersheds (e.g. Flanner et al., 2009;Painter et al., 2010;Oaida et al., 2015). Several studies have investigated the spatial and temporal variability of snow albedo using insitu data (Brock et al., 2000;Wuttke et al., 2006;Dumont et al., 2017) or satellite observations (Atlaskina et al., 2015;Naegeli & Huss, 2017). Snow spectral albedo generally depends in a complex way on several factors, including 1) the snow physical https://doi.org/10.5194/tc-2019-179 Preprint. Discussion started: 30 August 2019 c Author(s) 2019. CC BY 4.0 License. and chemical properties, mainly the Specific Surface Area of snow grains (SSA, Gallet et al., 2009) and the concentration of 40 light absorbing particle (LAP, Skiles et al., 2018), 2) the spectral and angular characteristics of the incident radiation (Warren, 1982), 3) the presence of macroscopic surface roughness (Kuhn, 1985;Warren et al., 1998;Mondet and Fily, 1999). The first two points have been thoroughly studied, showing that for a smooth surface, snow albedo decreases as SSA lowers (coarsening snow granularity) and with a higher sun elevation (i.e. a decrease in solar zenith angle), both of which lead to an increased absorption (Warren et al., 1998, Kokhanovsky andZege, 2004). Nevertheless, the effects of roughness are often neglected due 45 to the difficulty to characterise the actual surface roughness within the footprint of the sensor.
Snow-covered surfaces often exhibit macroscopic roughness, resulting from snow transport or erosion by the wind or snow melting (Filhol and Sturm, 2015). In Antarctica, roughness height ranges from a few centimetres to a few meters (Warren et al., 1998;Wuttke et al., 2006), and the features' axis is usually aligned along the prevailing wind direction (Furukawa et al., 1996), whereas in alpine areas the spatial distribution of macroscopic roughness mainly depends on topography, which drives 50 wind direction and its intensity (Naaim-Bouvet et al., 2011). Kuhn (1974) was the first to report a reduction of the forward peak of the Bidirectional Reflectance Distribution Function (BRDF) over a sastrugi field, and attributed this fact to shadows when the solar azimuth angle is perpendicular to the sastrugi. This motivated further studies that showed a systematic albedo decrease in presence of roughness (Caroll and Fitch, 1981;Leroux and Fily, 1998;Corbett and Su, 2015). The amplitude of the reduction in albedo depends on illumination conditions, snow properties, the size and the orientation of roughness features 55 (Hudson and Warren, 2007;L'Hermitte et al., 2014). For instance, in high altitude mountain glaciers, the presence of penitentes, which can reach several meters in height (Lliboutry, 1953), causes a measured albedo decrease of 8-10% (Corripio and Purves, 2006). These studies underlined the difficulty of precisely quantifying the impact of roughness since the illumination conditions and snow properties also vary during albedo measurements, making it difficult to evaluate the reduction in albedo due to roughness only. A protocol was proposed by Kuchiki et al. (2011) using a controlled environment where the precise roughness 60 shapes, orientation and dimensions, snow properties and illumination conditions were known. Over a manually-created artificial roughness field, they showed that the hemispherical-directional reflectance (HRDF) factor varies by more than ± 50% relative to a smooth surface. Nevertheless, they did not acquire albedo measurements, i.e. bi-hemispherical reflectance. Warren et al. (1998) showed that the albedo decrease over a roughness field is controlled by two effects: 1) a decrease in the insolation-weighted average incidence angle relative to a flat surface (further referred to as the effective angle effect), and 2) 65 multiple reflections in the concavities. The first effect is explained by the fact that the sides of the roughness shapes facing the sun experience stronger radiation with a smaller angle than the solar zenith angle which enhances absorption in the case of snow surface (Warren, 1982), and the sides facing away from the sun receive less radiation due to shadows or grazing angles.
The insolation-weighted average albedo is therefore reduced relatively to a flat and smooth surface (Warren, 1982;Kokhanovsky and Zege, 2004). The effective angle effect varies with the shape, size, and orientation of the roughness features 70 (Carroll and Fitch, 1981;L'Hermitte et al., 2014), and is significant under direct illumination and for low sun elevations only (Warren et al., 1998). The second effect of roughness involves multiple reflections cause by the trapping of photons between roughness shapes (Pfeffer and Bretherton, 1987). Over a smooth surface, a photon only hits the surface once and is either absorbed or reflected to the sky. Over a rough surface, photons can not only be absorbed or reflected to the sky, but they can also be reflected back to the surface. In this latter case, they have another probability to be absorbed, at every hit. This results 75 in a systematic increase in absorption, and thus a decrease in albedo. The impact is maximal when the probabilities of reflection and absorption are balanced, i.e. for intermediate values of albedo (close to 0.5 in the near infrared at 700-1100 nm). Instead in the visible where albedo is close to 1, the probability of absorption is too low to trap the photons, and oppositely in the midinfrared where the albedo is close to 0, the impact of multiple reflections is negligible. This trapping effect operates under direct and diffuse illumination. Although these two effects have never been quantified separately, Warren et al. (1998) 80 suggested to acquire measurements in diffuse illumination to estimate the impact of multiple reflections only. Photometric models based on analytical equations were developed to simulate the effects of roughness on albedo using idealized geometric shapes (Carroll, 1982;Pfeffer and Bretherton, 1987;Wendler and Kelley, 1988;Leroux and Fily, 1998;Cathles et al., 2011;Zhuravleva and Kokhanovsky;. Leroux and Fily (1998) predicted a decrease in albedo over a sastrugi field of 5 -9 % at 900 nm, depending on the sastrugi orientation with respect to the sun position. Despite their 85 interest to draw general conclusions on the albedo sensitivity to roughness characteristics, these models are of limited interest for real roughness features due to the idealization of the shapes (Warren et al., 1998). In addition, they use the Lambertian approximation to represent the surface reflectivity, and do not consider the intrinsic BRDF of the snow, meaning that they cannot simulate the effective angle effect. To explore the real impact of surface roughness, a 3D radiative transfer model is needed. Monte Carlo photon transport algorithms are convenient approaches (Lafortune, 1995;O'Rawe, 1991;Iwabuchi, 2006;90 Kuchiki et al., 2011). However, most studies using these numerical methods aim to evaluate the BRDF or HRDF instead of albedo, as their application domain was remote sensing (Kuchiki et al., 2011;L'Hermitte et al., 2014;Corbet et al., 2015).
The aims of this paper are two-fold: 1) to quantify the impact of surface roughness on snow albedo, as a function of roughness features, illumination conditions and snow properties, and 2) to assess the respective roles of the effective angle effect and multiple reflections with a new model able to represent surface roughness. Firstly, we collected albedo measurements in 95 controlled experiments following the idea of Kuchiki et al. (2011). We produced various artificial rough surfaces during four field campaigns in the French Alps in 2018 and 2019 (Sect. 2). In each experiment, albedo measurements were acquired for several illumination conditions and with numerous geometrical characteristics at the surface. Observations were also acquired over nearby smooth surfaces to serve as references. Secondly, we developed a new model based on the Monte Carlo photon transport method, the Rough Surface Ray Tracing model (RSRT), to simulate albedo by considering surface roughness (Sect. 100 3). RSRT was evaluated using the albedo observations (Sect. 4.1). In Section 4.2, the model was used to explore the albedo sensitivity to surface roughness according to SSA, terrain slope, roughness orientation and solar zenith angle. The model was applied to assess the respective roles played by the effective angle effect and multiple reflections (Sect. 4.3). At last, the sensitivity of the net short wave radiation to the presence of surface roughness is discussed to estimate the potential impact on the surface energy balance (Sect. 4.4). 105

Field experiments
In situ measurements of albedo were acquired in the French Alps over smooth and rough snow surfaces. This section details how the rough surfaces were created, and measurements acquired in the field.

Artificial rough snow surfaces
Artificial rough snow surfaces were created by delineating squares of 2.5 x 2.5 m 2 . Roughness features were manually created 110 on natural smooth surfaces, by varying their number and orientation. The features were produced parallel to each other, regularly spaced with a period Λ, and with an azimuth angle φr, taken clockwise from the North. The roughness orientation with respect to the solar azimuth angle (φs) was defined by Δφr, the difference φs -φr. Figure 1 shows the experimental setup and the variables involved. Each surface was characterised by its aspect, its slope, and its roughness properties (number, shape, size and orientation). Two types of experiments were performed: 115 a) Sensitivity to the fraction of roughness features: The fraction of roughness features in the 2.5 x 2.5 m 2 area is described with the width-to-period ratio η (i.e. η = W/Λ, expressed in percentage, where W is the width of roughness shapes). The albedo sensitivity to η was studied during two experiments at the Col du Lautaret site (45°2' N, 6°2' E, 2100 m a.s.l.) over two different dates (06 April 2018 and 17 April 2018), respectively called experiments A and B. Figure 2 illustrates the field experiment A, and Table 1  parameters for each studied surface. Snow albedo was first measured over the smooth surface (called A-smooth and B-smoothdry in Table 1), and then the roughness shapes were created in the smooth surface by uniformly pressing a rectangular metal bar into the snow (H = 2 cm depth and W = 4 cm width), in the North-South direction (φr = 0°-180°). The rectangular shapes were created with a period Λ = 55 cm (5 shapes over 2.5 m, η = 7 %). After albedo measurements were acquired, identical rectangular shapes were added to reach a period Λ = 30 cm (10 shapes over 2.5 m, η = 13 %), then Λ = 15 cm (20 shapes over 125 2.5 m, η = 27 %) ( Figure 2 and Table 1). Because it takes approximately one hour to make a series of measurements, the increasing fraction of roughness features is correlated to solar zenith angle (Өs) variations that also change albedo. To attempt to decouple the two effects, experiment A was conducted when the sun was going down (Өs went from 56.6° to 63.7°), whereas in experiment B, the sun was going up (Өs went from 56.0° to 40.0°).
Other changes may also occur during that time. In experiment B for instance, melting was observed on the B-η27% surface 130 (the sun was close to the nadir), which leads to an increase in snow wetness and a decrease in surface SSA compared to the Bsmooth-dry surface analysed at the beginning of the day. To allow more reliable comparisons, we simultaneously measured albedo over the B-η27% surface, and a nearby smooth surface (called B-smooth-wet in Table 1). b) Sensitivity to the roughness orientation: The albedo sensitivity to roughness orientation was studied with two experiments at the Arcelle site (45°6' N, 5°52' E, 1729 m 135 a.s.l.) over two dates (11 January 2019 and 22 February 2019), respectively called experiments C and D. The roughness shapes were triangular, H = 6 cm depth and W = 7 cm width, and created with a period Λ = 11 cm (η = 63 %). Fig. 1b and Table 1 detail the experimental setup.
In experiment C, measurements were simultaneously acquired every 20 minutes over a surface with roughness features oriented at φr = 90° (called C rough 90°), another one with roughness features at φr = 0° (called C rough 0°), and a smooth 140 surface for reference (called C smooth). In experiment D, only two surfaces were compared every 20 minutes: a rough surface with roughness features at φr = 90° (called D rough 90°), and a smooth surface (called D-smooth). For both experiments, studied surfaces were close enough to consider that snow properties evolved with the same dynamics. Note that it took about up to 5 minutes to acquire one set of albedo measurements, and to move to the next surface. Measurements were acquired all day in the experiment C (sun going up and down), and during the morning in experiment D (sun going up only). 145 The albedo sensitivity to roughness features is quantified by comparing rough and smooth surfaces for each experiment.   (Table 1), from left to right: A-smooth, A-η7%, A-η13%, and A-η27% sites.

Spectral albedo measurements
Spectral albedo, or more precisely the bi-hemispherical reflectance (Schaepman-Strub et al., 2006), is the ratio of the upwelling and the downwelling spectral irradiance. Snow spectral albedo measurements were acquired with the Solalb instrument, a manual version of the albedometer AutoSolexs described by Picard et al. (2016). Solalb is a hand-held instrument using a single light collector with a near-cosine response and equipped with an inclinometer located at the end of a 3 m boom. 160 The boom was rotated by the operator to successively acquire the downward and upward solar radiation with a horizontal sensor (± 0.1° accuracy). This operation usually takes up to a maximum time of 30 seconds. Variations of incident illumination caused by clouds between two acquisitions were also measured with a photodiode receiving ambient radiation. Only spectra with stable incident illumination within 1 % were selected. Spectra were acquired over the 400-1050 nm wavelength range with an effective resolution of 3 nm. The height of the sensor impacts the measured roughness effects, by changing the footprint 165 of the sensor (L'Hermitte et al., 2014). To study this sensitivity, albedo was measured with sensor heights of 45 cm, 55 cm and 65 cm, in experiments A and B (not shown). We found a weak influence on measured albedo (0.4 ± 0.5 % of differences between spectra), showing that this sensitivity was negligible given the type of roughness considered here, and the sensor's height. Therefore, the sensor was set to 65 cm high for all experiments. At this height, the footprint is about 2.3 x 2.3 m 2 (99 % of the signal is coming from a viewing angle of 60°, Picard et al., 2016). The ratio of diffuse-to-total irradiance (rdiff-tot) was 170 also measured shortly after the albedo measurement by screening the sun to record the diffuse irradiance, the total irradiance being measured with the sensor looking upward. Post-processing was applied to each acquired spectrum following Picard et al. (2016). This includes dark current correction, considering the integration time, and the correction of the collector angular responses. The observed apparent albedo, called αobs, is the processed spectrum measured with Solalb, considering the sensor in a horizontal position (Sicart et al., 2001). The 175 accuracy of αobs mainly depends on that of the levelling of the arm. To estimate αobs uncertainties, measurements were duplicated three times for 6 different sites. A maximal variation of 1.6 % was estimated between the αobs spectra acquired in same field conditions.

Snow surface properties
Snow SSA was measured at the surface using the Alpine Snowpack Specific Surface Area Profiler (ASSSAP) instrument 180 that has an accuracy of 10 % . For the two experiments A and B, only one surface SSA was measured, in the middle of the experiment (corresponding to η = 13%), and the SSA was assumed to be constant throughout the experiments (3 hours). The albedo sensitivity to SSA variations and associated uncertainties is discussed in Sect. 4.2 in order to untwine these contributions from those of roughness. For the two experiments C and D, SSA were measured at the surface at each albedo acquisition. 185 Moreover, albedo are highly impacted by the surface slope because of the change in the local incident solar angle (Grenfell et al., 1994). It has been shown that a slope of 2° facing the sun induces a variation in albedo by up to 5 % over a smooth surface (Dumont et al., 2017). An inclinometer fixed at the end of a 2 meter ruler was used to measure the slope in the sensor's footprint. The aspect of the slope is defined as the azimuth of the direction of the steepest slope, clockwise in degrees from North. However, as the studied surfaces were chosen as flat as possible, the steepest inclination was not visually detectable. 190 Thus, a first inclination measurement was acquired with the ruler parallel to the roughness features, and a second one by rotating the ruler by 90°, in order to estimate the normal of the surface ( ⃗ = (nx, ny, nz)). The slope and the aspect were deduced as follow:

θn = arccos(nz)
[1], 195 where θn is the steepest slope angle, and φn is the aspect of the slope. In this study, all surface slopes were below 5°. The uncertainty of slope measurements was estimated of ± 1° due to natural ripples of the studied surfaces. The impact of this uncertainty in our roughness analysis is discussed in Sect. 4.2.
To limit the scope of this study, snow LAPs were not measured although they strongly lower the spectral signature in the visible range, especially at the end of the season when the concentration of impurities is high at the surface (Flanner et al., 200 2009). It was the case for experiments A and B (measurements acquired in April). Figure 3 shows the spectrum measured over the A-smooth surface. The albedo decrease in the 400-600 nm range is a clear signature of a high LAP concentration. This sensitivity is well described in Dumont et al. (2017) and is out of the scope of this paper. However, to minimize this contribution, we chose to quantify effects of roughness in the 600-1050 nm wavelength range. Statistical results are given at 700 nm and at 1000 nm, to represent the visible and the NIR domain, respectively.   (Table 1).

A 3-D Monte Carlo radiative transfer Model 210
The RSRT model was developed to simulate snow albedo considering macroscopic surface roughness. This combines both 1) the asymptotic radiative transfer theory (Sect. 3.1) to compute the spectral albedo each time a photon hits the modeled surface and 2) a Monte Carlo technique (Sect. 3.2) to estimate the geometric effects introduced by roughness and represented with a 3-D mesh of the studied area. Section 3.3 details the simulation framework and the sensitivity analysis. A simple approach is applied to estimate the impact of roughness on the quantity of energy absorbed in the snowpack (Sect. 3.4). 215

Asymptotic radiative transfer theory
In the RSRT algorithm, an ensemble of photons is launched over a modeled surface. This surface is represented with a triangular mesh composed of small facets. Both the spectral albedo and the BRDF distribution are computed for each facet hit by a photon. The Asymptotic Radiative Transfer theory (ART) provides analytical equations to estimate spectral albedo for highly reflective materials, which applies well to snow in the visible and the NIR domains, typically from 400 nm to 1100 nm 220 (Zege et al., 1991;Kokhanovsky and Zege, 2004). Several models use this theory (Negi et al., 2011;Libois et al., 2013;Wang et al., 2017), which is based on three assumptions: 1) the snowpack is represented with vertically and horizontally homogeneous plane-parallel layers, 2) the surface is perfectly smooth and horizontal (flat), 3) single scattering albedo and the snow phase function are described with the asymmetry factor, g, the absorption enhancement parameter, B, and the SSA of the snow. The albedo simulated with the ART theory have shown a good accuracy compared to observations over smooth surfaces 225 (Dumont et al., 2017;Wang et al., 2017). The facets of the mesh are small enough to be considered as smooth surfaces. The direct and the diffuse part of the albedo at the wavelength λ and θs, αdir(λ, θs) and αdiff(λ), are estimated with Eq. (3) and (4): where B and g are assumed to be constant (B = 1.6 and g = 0.86, Libois et al., 2014b), ρice = 917 kg m −3 is the bulk density of 230 ice at 0°C, and γ(λ) is the wavelength-dependent absorption coefficient of ice, taken from Picard et al. (2016) here.
The albedo of a flat smooth surface obtained with ART (αflat) at wavelength λ and at θs is deduced as follows: where rdiff-tot(λ, θs) is the ratio of diffuse-to-total illumination at wavelength λ and at θs.
These formulations apply to a strictly leveled terrain (better than 0.5°). In the case of a tilted surface, the apparent albedo 235 acquired with a sensor placed horizontally is different from αflat because the surface receives sun radiation with a different https://doi.org/10.5194/tc-2019-179 Preprint. Discussion started: 30 August 2019 c Author(s) 2019. CC BY 4.0 License. incidence angle and is viewed with a reduced solid angle by the sensor. To account for the slope and compute the apparent albedo of a titled smooth surface, called αsim,smooth, a K factor is applied (Dumont et al., 2017), such as: and: 240 where ̃ is the effective θs modified with the slope.
Following the ART theory, Kokhanovsky and Zege (2004) (further referred as the KZ04 approximations) estimated the snow BRDF distribution by calculating reflectance over a hemisphere with the reflection function of a semi-infinite medium: where the function R0( , , ) is the reflection function at 0 = 1 (Kokhanovsky, 2013), with 0 the single scattering albedo. is the relative azimuth angle, is the cosine of the viewing zenith angle, is the cosine of the solar zenith angle, and A is estimated as follows: ks and kv are called the escape functions, and are given by Kokhanovsky (2003)  (1 + 2 ) [11]..

Algorithms and model architecture
The Monte Carlo photon light transport algorithm propagates a large number of photons from their source to termination (i.e. 255 escape from the scene). This algorithm is described in detail in the following. A photon is a particle of light carrying a flux and described by its power (intensity), its origin , and its propagation direction . Each photon starts its trajectory with an intensity equal to 1 (unitless quantity of energy), and a direction described with the couple (θs, φs) given as input. Photons are either absorbed or reflected at each hit according to the facet albedo (Iwabuchi, 2006). A flow chart of a photon path as computed with RSRT is presented in Figure 4. This is computed in four main steps. 260 Step 1: Estimate the next intersection of the photon with the mesh of the surface (called "hit"). The Bounding Volume Hierarchies (BVH) technique (Ize, 2013) is used to efficiently search for the first facet in the photon propagation direction.
The precise intersection point within the facet is determined by applying the watertight ray/triangle intersection algorithm (Woop et al., 2013). If the photon hits a facet, its origin is updated on the intersected facet. The normal of the facet is estimated. If there is no hit, the photon escapes from the mesh, and depending on its direction (upward or downward), its 265 intensity is added to the down or up welling radiation bin (Fig. 4).
Step 2: Update the intensity. The photon intensity at hit n (called ip,n) is weighted by the spectral albedo accounting for the incoming direction angles αflat(λ, θi) as follows: ip,n+1 = ip,n αflat (λ, θi) [12], Two configurations are possible (Sect. 3.3). With the KZ04 configuration, the hit facet is considered as a snow surface and 270 αflat(λ, θi) is estimated by considering the local incident angle θi and snow properties (SSA, B, g) (i.e. with ART, Eq. (5)). With the Lambertian configuration, the hit facet is a Lambertian surface (i.e. isotropic diffusion), and the αflat(λ, θi) is a constant value equal to αflat(λ), given as input.
Step 3: Sample the outgoing direction. The most likely outgoing direction of the photon after a hit is estimated from the by the (cosθv, φv) couples. The outgoing direction is sampled from the BRDF distribution using a rejection algorithm as follows: in a first step, the azimuth φv is sampled from a uniform distribution between 0 and 2π, and cosθv with a uniform distribution between 0 and 1, so that the hemisphere is sampled with a cosine weighting distribution (Greenwood, 2002). In a second step, a probability of acceptance is given to each direction (θv, φv). This probability of acceptance is estimated by the BRDF value in this direction, normalized by the maximum value of the BRDF distribution. 280 Step 4: Update the direction . The new photon direction +1 after the hit n is updated as follows: With , and , the photon directions in the x and y axis before the hit n, respectively. 285 The algorithm returns to step 1 until the photon escapes from the scene (Fig. 4), or until its intensity is lower than a threshold (set to 0.01 in RSRT). To ensure an unbiased termination in the latter case, a 'Russian roulette' method is applied (Iwabuchi, 2006), which consists in: accepting or rejecting the termination with probabilities 1-p and p, respectively (p = 0.2 in RSRT).
In case of rejection, the weak intensity of the photon is rescaled by the factor 1/p, and the algorithm goes again to step 1. As explained by Iwabuchi (2006), the total energy is conserved for any p value, and this approach can be applied at any step of 290 the algorithm.
At the end of its path, the photon intensity is counted in: 1) the total upward intensity ( ↑ ) if the photon escapes with an upward direction, 2) the intensity lost downward if its final z axis direction is downward (this is possible with a tilted surface for instance). If the latter contribution is higher than 10 -3 for a horizontal rough surface, we consider that too many photons have been lost to output a realistic albedo, meaning that the simulation used a too wide radiation source or conversely a too 295 small mesh area.

Model simulations
RSRT is run here either by considering the snow surface 1) as Lambertian (Lambertian configuration), the albedo is not sensitive to the incident angle and each photon hitting the mesh is reflected with a constant facet albedo equal to αflat(λ), or 2) as a snow surface using KZ04 analytical equations (referred to as the KZ04 configuration). In this latter configuration, each photon hitting the mesh is reflected with αflat(λ, θi), which depends on the incident angle θi, SSA, B and g values, i.e. by 305 considering the intrinsic BRDF of the snow. The two configurations are compared in Sect. 4.3. The KZ04 configuration is used by default in for all other simulations to compare with observations. RSRT inputs are described in Table 2. Triangular meshes of rough surfaces are modelled by reproducing same linear shapes as those created in the field with an orientation φr, a height H, a width W, and spaced by a constant distance defined with the period Λ (as shown in Figure 1, with same values as in Table 1). Meshes have a spatial resolution of 1 cm and are produced 310 large enough to be considered as infinite (no edge effects). When a RSRT simulation is started, an ensemble of photons is first created on a horizontal plane above the surface mesh and distributed quasi-randomly to produce a parallel source. The size of the photon ensemble is set to 10 6 photons as a compromise between the computing time and a good representation of the emission source. The direction of propagation of the ensemble of photon is initialized with the solar zenith and azimuth angles given as inputs. 315 RSRT outputs the snow spectral albedo, either in direct or diffuse illumination conditions: αdir,rough(λ, θs) and αdiff,rough(λ), respectively, considering that the plane of the mesh is perfectly flat. To take into account the slope, αdir,rough(λ, θs) and αdiff,rough(λ) are combined with Eqs. (6) and (7) to simulate the apparent snow albedo of a titled rough surface, called αsim,rough(λ, θs).

Evaluation of simulations
The evaluation of simulations was treated over a set of N observations using the Root Mean Square Deviation (RMSD), defined as follow: with αsim,i(λ, θs) the i th simulation (either αsim,smooth or αsim,rough) at the wavelength λ and θs, and αobs,i(λ, θs) the i th measured albedo. 325 The performances of αsim,smooth(λ, θs) and αsim,rough(λ, θs) are compared to evaluate the accuracy gain acquired by taking into account surface roughness. We further called ̅̅̅̅̅̅̅̅ , the RMSD(λ, ) averaged over the 600-1050 nm range for one spectra.

Impact of uncertainties
Albedo observations may have been affected by uncertainties or unmeasured variations. To investigate the potential impact, we conducted the following simulations. Firstly, SSA may have varied over time in experiments A and B, whereas albedo was 330 simulated with a constant SSA. In order to estimate these variations, we retrieved SSA at the beginning of the experiments from albedo observations over the smooth surfaces, by fitting αsim,smooth with αobs using the same approach as described by Libois et al. (2015). RSRT was then run by considering retrieved SSA (SSAr) for simulations over A-smooth and B-smoothdry surfaces, and the measured SSA (SSAm) for simulations over the rough surfaces. Results are studied at 1000 nm where the sensitivity to SSA is larger (Dominé et al., 2006). Secondly, the difference between retrieved and measured SSA may be related 335 to the uncertainty in SSA measurements. We explored the impact of SSA uncertainties with RSRT simulations by varying SSAm by ± 10 % over the rough surfaces at 1000 nm. Thirdly, the impact of slope uncertainties was studied with RSRT simulations by varying the slope of the rough surfaces by ± 1° in experiments C rough 90° and D rough 90° at 1000 nm. We used C and D experiments only since observations over the rough and smooth surfaces were acquired simultaneously, with similar θs values (to not influence ̃ the effective θs modified with the slope). 340

Analysis of processes introduced by surface roughness
The variations of illumination conditions and SSA may attenuate or accentuate roughness effects by playing a role either in the effective angle effect, or in multiple reflections. We thus investigated separately these effects as a function of illumination conditions and SSA to better characterize roughness effects.
The effective angle effect is the alteration of the local incident angle over roughness shapes. It was simulated with RSRT using 345 the KZ04 configuration (albedo varying with θs) and by requiring that photons hit the surface only once, i.e. without multiple reflections. The total upward and downward intensities were then added to count all the photons that have not been absorbed after the first hit. We also conducted same simulations with the Lambertian configuration to check there was no angular dependence. These simulations were performed with various illumination conditions. The effect of multiple reflections caused by the photon trapping depends on the albedo value. While the effective angle effect 350 is significant under direct sunlight only, this second effect is significant both under direct and diffuse illumination (Warren et al., 1998). Therefore, it was simulated by running RSRT under diffuse sunlight. Simulations were conducted for various SSA.

Results and discussion
We explore the albedo sensitivity to macroscopic surface roughness through three questions: 1) is it possible to quantify the at each hit) relative to a smooth surface (only one hit), causing the observed albedo decrease. Note that the pattern of the measured spectra between 600-700 nm is probably led by a high concentration of snow LAPs (not visible to the naked eye on 365 the field).
By neglecting the presence of roughness, αsim,smooth spectra increase due to larger θs from A-smooth to A-η27% (the sun went down). Indeed, photons penetrate less deeply into the snowpack as they enter with a grazing angle (large θs). They encounter the first scattering event near the surface and has more probability to escape compared to a photon penetrating deeper with a low θs (Carroll and Fitch, 1981;Warren, 1982). The bias between αobs and αsim,smooth becomes higher while η increases ( ̅̅̅̅̅̅̅̅ = 370 0.06 when η = 27%). On the contrary, by considering the presence of roughness, αsim,rough spectra follow the observed trend, and decrease while η increases (even if θs becomes larger). Simulations are strongly improved compared to αsim,smooth that neglects surface roughness ( ̅̅̅̅̅̅̅̅ = 0.02 when η = 27%). Nevertheless, reductions in albedo are not clearly quantifiable from one rough surface to another (from A-η7% to A-η13%, for instance). Results are thus further plotted at two wavelengths only, 700 nm and 1000 nm. 375
In experiment A, the stronger Δαobs decrease is of 0.03 at 700 nm, and of 0.07 at 1000 nm, between A-η27% (η = 27 %, θs = 63.6°) and A-smooth (η = 0 %, θs,0 = 56.7°) (Fig. 6a and 6c). Experiment A-η7% shows that even a low fraction of roughness features causes an albedo decrease of 0.02 (0.03) compared to that of a smooth surface at 700 nm (1000 nm whereas θs increases, showing that albedo is more sensitive to roughness effects than to θs variations here. This result highlights 390 the need to consider the presence of roughness in albedo simulations. While θs becomes larger, Δαsim,smooth increases by 0.02 at 700 nm, and by 0.03 at 1000 nm, between A-η27% and A-smooth ( Fig. 6a and 6c). By considering roughness in simulations, Δαsim,rough decreases by 0.01 at 700 nm, and by 0.03 at 1000 nm between A-η27% and A-smooth. Simulations are thus improved by accounting for surface roughness but RSRT (i.e. αsim,rough) underestimates by almost a factor 2 the observed albedo reduction. 395 In experiment B, Δαobs shows a strong decrease of 0.11 at 700 nm, and 0.15 at 1000 nm, between B-η27% (η = 27 %, θs = 39.9°) and B-smooth-dry (η = 0 %, θs = 55.4°) (Fig. 6b and 6d). In this experiment, αobs decreases due both to the θs decrease (the sun went up, Sect. 2.1) and the η increase. To remove θs contribution, we use the Δαsim,smooth trend that depends on θs variations only: Δαsim,smooth lowers when θs decreases and the reduction is half of that of Δαobs (0.05 at 700 nm, and 0.07 at 1000 nm when η = 27%). Therefore, half of the αobs decrease is attributable to the decrease in θs, and the other half to the presence 400 of roughness.
Δαsim,rough decreases by 0.07 at 700 nm, and 0.11 at 1000 nm, between B-η27% and B-smooth-dry ( Fig. 6b and 6d). Simulations are more consistent with observations by considering the presence of roughness, but the simulated decrease is still underestimated compared to measurements, as for experiment A.
To accurately quantify roughness effects on albedo, it is important to compare rough and smooth surfaces for similar snow 405 and illumination conditions. This is why we simultaneously measured albedo over B-η27% (η = 27%, θs = 39.9°) and a nearby smooth surface (the B-smooth-wet surface: η = 0%, θs = 36.4°, Table 1 and Fig. 6b and 6d). The concurrent measurements show a decrease by 0.05 at 700 nm, and 0.07 at 1000 nm. This reduction is solely attributable to the presence of roughness. It is similar to the Δαobs decrease by subtracting the Δαsim,smooth that is caused by the θs decrease only (Fig. 6).
For both experiments, observations show that the albedo decrease is stronger when 1) the number of roughness features is 410 larger, and 2) at the longer wavelengths. As albedo is lower in the NIR domain, the impact of multiple reflections is stronger.
Indeed, the effect of multiple reflection is more important for intermediate values than for albedo close to 0 or 1 (i.e. systematic absorption or reflection, Warren et al., 1998). Table 3 shows that RMSDs estimated for αsim,smooth increase with the fraction of roughness features, and are higher at 1000 nm than at 700 nm. By considering roughness in the RSRT model, the accuracy of αsim,rough is improved by about a factor 2 at 415 700 nm and 1000 nm compared to αsim,smooth (average αsim,rough RMSD of 0.02 at 700 nm and 1000 nm).

Sensitivity to the roughness orientation
In experiments C and D, albedo measurements are simultaneously acquired over a rough surface and a nearby smooth surface for multiple illumination conditions, every 20 minutes. A spectral analysis of the impact of roughness is first discussed, using 430 a randomly chosen illumination condition for each experiment. Then, the albedo sensitivity to the roughness orientation with respect to the solar azimuthal angle (Δφr) is investigated at two wavelengths. αobs spectra show a significant decrease caused by the presence of surface roughness (~ -0.05 on average, Fig. 7a and 7d), more 435 pronounced in the NIR domain. αsim,smooth spectra are identical over the smooth and rough surfaces since illumination conditions are similar (dotted lines in Fig. 7b and 7e). For these illumination conditions, αsim,smooth spectra over the rough surfaces show a strong deviation compared to αobs, with a ̅̅̅̅̅̅̅̅ 0.04 and 0.05 for C rough 0° and C rough 90°, respectively, and of 0.06 for D rough 90°. By considering roughness, αsim,rough spectra reproduce the observed decrease (~-0.05, Fig. 7c and 7f). Simulations are improved compared to αsim,smooth, with a ̅̅̅̅̅̅̅̅ of 0.02 and 0.04 for C rough 0° and C rough 90°, and of 0.04° for D-440 rough90.  Effects of roughness with varying illumination conditions: Figure 8 shows the change in albedo as a function of Δφr at 700 nm and 1000 nm. Δαobs is the difference between αobs simultaneously acquired over the rough and the smooth surfaces. Similarly, Δαsim,rough is the difference between αsim,rough 450 simulated over the rough surface and αsim,smooth simulated over the smooth surface, at same illumination conditions. Thus, the change in albedo is not correlated to θs here, but to φs that leads the roughness orientation with respect to the sun position.
When roughness features are parallel to the sun (when Δφr = 0° for C rough 0°, Fig. 8a and 8d), αobs decreases of 0.01 at 700 nm, and of 0.08 at 1000 nm, relative to a smooth surface ( Fig. 8a and 8d). The impact becomes larger when the roughness orientation is perpendicular to the sun (when Δφr = 90° for C rough 90°, Fig. 8b and 8e), with an αobs decrease of 0.03 at 700 455 nm and of 0.10 at 1000 nm. For this case, the reduction in albedo is 20 % stronger when roughness features lie perpendicular to the sun than when they are parallel. When the sun elevation is low, if the roughness orientation is perpendicular to the sun, the effective incident angle over sides facing the sun is decreased compared to that of a smooth surface. In addition, the fraction of shadow is higher when Δφr = 90°. This effective angle effect leads to an average decrease in snow albedo relative to a smooth surface. However, for the C rough 90° experiment, Δφr varies from 50° to 122° and Δαobs does not show a strongest 460 albedo reduction around 90°. The SSA is particularly high for this experiment (~100 m 2 kg -1 ) and it may explain the albedo insensitivity to weak variations of roughness orientation. In experiment C, Δαsim,rough variations reproduce well the Δαobs decrease, with the same order of magnitude: the average decrease is of 0.01 at 700 nm and 0.08 at 1000 nm for C rough 0°, and of 0.02 at 700 nm and 0.10 at 1000 nm for C rough 90° ( Fig. 8a and 8d). 465 In experiment D rough 90°, measurements were acquired in the morning, so Δφr varies from 42° to 72° (Fig. 8c and 8f). We measured an average Δαobs decrease of 0.02 at 700 nm and 0.09 at 1000 nm, which is in agreement with results found for C rough 90°. The observed pattern of Δαobs shows oscillations, probably caused by differences in snow properties between the smooth and the rough surfaces. SSA was measured over the smooth surface to be representative, while SSA values over rough surfaces evolved with spatial variations according to the received illumination. In addition, for this experiment, melting was 470 observed at the surface, resulting in a smoothing of our roughness shapes during the day. This is certainly why in Fig. 8c and   8f, the Δαobs increases from 42° to 72°, while in theory it should decrease when Δφr approaches 90°. Fig. 8c and 8f shows that αsim,rough overestimates by almost a factor 2 the reduction in αobs: the average Δαsim,rough decrease is of 0.06 at 700 nm, and of 0.15 at 1000 nm. By considering roughness shapes constant along the day, Δαsim,rough decreases from 42° to 72° (i.e. Δφr gets closer to 90°). 475 Considering all observations of experiments C and D, αsim,rough is significantly improved compared to αsim,smooth, with an averaged RMSD of 0.03 at 700 nm and 0.04 at 1000 nm (against an averaged RMSD of 0.07 at 700 nm and 0.09 at 1000 nm for αsim,smooth, Table 2).
To sum up, observations show that an increase of the number of roughness features leads to a larger reduction in αobs, with a 480 higher sensitivity in the NIR domain. Roughness effects are also larger when the roughness orientation is perpendicular to the sun rather than parallel. The maximal αobs decrease is of 0.03 at 700 nm, and 0.10 at 1000 nm, for C rough 90°. Considering all observations, αsim,rough are improved by a factor 2 at 700 nm and 1000 nm compared to αsim,smooth, with an average RMSD of 0.03 at 700 nm and 0.04 at 1000 nm (Table 2). αsim,rough shows an underestimation of the observed albedo decrease, but observations may have been affected by uncertainties or unmeasured variations. 485

Analysis of uncertainties
In a first step, we explore the possible SSA variations in experiments A and B, and the impact on snow albedo. In a second step, we integrate SSA and slope uncertainties in our roughness analysis.

Sensitivity to SSA
By fitting αsim,smooth and αobs over the smooth surface (Sect. 3.3), we estimate an SSA (written SSAr) of 9.4 m 2 kg -1 over A-495 smooth and 5.3 m 2 kg -1 over B-smooth. Measured SSA (SSAm) are equal to 7.4 m 2 kg -1 over A-η13% and 4.5 m 2 kg -1 over B-η13%. Hence, for both experiments, there is a decrease in SSA from the beginning (smooth surface, η=0%) to η = 13%, which is compatible with the observation of melt at the surface during these two experiments performed in April. Grenfell and Maykut (1977) explained that snow albedo decreases when liquid water replaces air between ice grains. The refractive index of the water being very close to that of ice, this results in an increases the effective grain size. 500 To explore the impact of a decreasing SSA on albedo, RSRT is run by considering SSAr for simulations over A-smooth and Bsmooth-dry surfaces (SSAr equal to 9.4 m 2 kg -1 and 5.3 m 2 kg -1 , respectively), and SSAm for simulations over the rough surfaces (from η = 7 % to 27 %, SSAm equal to 7.4 m 2 kg -1 and 4.5 m 2 kg -1 , Table 1). Results are presented in Figure 9, where Δαsim,rough,ssa is the difference αsim,rough (λ, θs, SSAm) -αsim,smooth(λ, θs,o, SSAr). Compared to Δαsim,rough (i.e. constant SSA), the Δαsim,rough,ssa decrease is multiplied by a factor two by considering both the increase in the fraction of roughness features, and the SSA 505 decline from 9.4 to 7.4 m 2 kg -1 (-15 %) for experiment A, or from 5.3 to 4.5 m 2 kg -1 (-21 %) for experiment B (Fig. 9a and 9b).
Δαsim,rough,ssa reproduces well the Δαobs decrease, with the same order of magnitude. Therefore, both SSA variations and roughness effects overlap and lower snow albedo in these two experiments, making it difficult to accurately isolate roughness effects.
Differences between retrieved and measured SSA may be explained by the uncertainty in SSA measurements. The impact of 510 SSA uncertainties is investigated by varying SSA by ± 10 % in RSRT simulations. Obtained values range within the grey shade shown in Figure 9. Experiment C has a large SSA (~ 100 m 2 kg -1 ), typical of fresh fallen snow, and SSA uncertainties weakly affect Δαsim,rough (Fig. 10c). On the contrary, a variation of ± 10 % in SSA strongly impacts the experiments with a low SSA: Δαsim,rough,ssa varies between 0.05-0.10 in experiment A-η27% (Fig. 9a), between 0.11-0.16 in experiment B-η27% (Fig. 9b), and between 0.13-0.18 in experiment D when Δφr = 72° (Fig. 9d). The reduction in albedo is stronger when SSA is lower due 515 to higher absorption. More precisely, the grains at the surface control the first scattering event and large-coarse grains (i.e. low SSA) are both more absorptive and more forward scattering relative to fine grains since photons have to pass through longer paths in ice before being potentially scattered at the ice-air interfaces (Warren et al., 1998). Dominé et al. (2006) have shown that the SSA-albedo relationship is non-linear and that albedo varies weakly in the NIR domain when SSA > 30 m 2 kg -1 , while it is highly sensitive to SSA variations for SSA values below 10 m 2 kg -1 . Hence, while a large SSA leads to a weaker impact 520 of multiple reflections (high albedo), the impact of the photon trapping is more important at low SSA (<10 m 2 kg -1 ). Experiment D highlights that the impact of SSA uncertainties in albedo is linked to the roughness orientation (Fig. 9d). Albedo is twice as sensitive to SSA when Δφr = 72° as when Δφr = 42°. This is caused by the effective angle effect introduced by roughness: photons penetrate deeper over sides facing the sun when the roughness orientation is perpendicular to the sun (lower incident angle) than if it was oblique or parallel. When SSA is low, absorptions increase and a photon has larger probability to 525 be absorbed by penetrating deeply in the snowpack. Hence, the effective angle effect is more pronounced when roughness orientation is perpendicular to the sun and for low SSA. compensate for the albedo reduction caused by the presence of roughness. We calculated SSA retrievals for experiments C and D at each Δφr by fitting αsim,smooth and αobs acquired over the rough surfaces. Compared to measure SSA taken over the smooth surface, results demonstrate that roughness introduces a significant underestimation of the retrieved SSA, reaching 21 % for the roughness features considered here. Thus, it is important to use a model considering roughness to retrieve accurate SSA 535 from albedo observations.

Sensitivity to the surface slope
The impact of slope uncertainties is explored by varying the slope by ± 1° for simulations over the rough surfaces for experiments C and D at 1000 nm (Sect. 3.3). Obtained values range within the grey shade shown in Figure 10. The albedo sensitivity to the slope depends of the slope aspect φn, with respect to the solar azimuthal angle φs, since the aspect controls the change in the incident angle ( ) relative to Θs. Eq. (6) and (7) show that slopes have no impacts on albedo if the aspect is 550 perpendicular to the solar azimuthal angle ([φs -φn] = 90° or 270°), while impacts change rapidly when φn becomes parallel to φs ([φs -φn] = 0° or 180°). In addition, between [φs -φn] = 0° or 180°, the albedo sensitivity to the slope is opposite whether the surface is smooth or rough. Indeed, over a titled smooth surface, if the slope aspect is opposite to that of the sun (φs -φn = 180°), the surface experiences a larger effective incident angle relative to a flat smooth surface, i.e. a higher albedo. Conversely, over a titled rough surface with roughness orientation perpendicular to the sun and φs -φn = 180°, roughness sides facing the 555 sun experience a lower effective incident angle relative to a flat rough surface, leading to a lower albedo. Fig. 10b illustrates this point for experiment D rough 90°: the albedo sensitivity is twice as strong when the slope direction is closer to 180° ([φs sensitivity to a change of the incident angle (Eq. (4)). Therefore, for low SSA, the impact of roughness on albedo is accentuated when the slope direction is opposite to the sun, and attenuated when the slope is facing the sun. 560 In experiment C (Fig. 10a), the albedo is highly sensitive to slope uncertainties (variations of 0.05-0.15). However, due to high SSA there is a low albedo sensitivity to the φs -φn angle (the effective angle effect is negligible). Therefore, the observed albedo sensitivity may be explained by a larger effect of multiple reflections, accentuated by the fact that Θs is particularly large for this experiment (> 60°).

565
To sum up, we have shown that the albedo sensitivity to roughness is larges when the SSA is low (< 10 m 2 kg -1 ), when roughness features are perpendicular to the sun, and when the surface aspect is facing away the sun. The impact of roughness if strongly linked to SSA values which affect: 1) the effective angle effect, since the decrease of the incident angle on roughness sides facing the sun has more consequences on the albedo when SSA is low (high absorption), 2) the impact of multiple reflections, which is larger when the probability of a photon to be absorbed or reflected is well balanced. To accurately quantify 570 roughness effects, it is crucial to measure SSA regularly (a small variation may overlap the roughness effects) and to determine the slope. In our experiments C and D, where SSA was measured at each albedo acquisition, we have shown that even considering uncertainties of ± 10 % of SSA and of ± 1° of slopes, roughness effects are significant and cause at least an albedo decrease of 0.06 in experiment C rough 90°, and of 0.11 in experiment D rough 90°, at 1000 nm.
575 Figure 10. Same as Figure 9c and 9d, except that the grey shades represent the range of Δα obtained by varying the slope by ± 1° in RSRT simulations, for a) experiment C rough 90°, and b) D rough 90°. φn is the aspect of the slope and φs is the solar azimuthal angle, separately given in Table 1. Vertical black lines indicates [φsφn] angles at the beginning and at the end of experiments (and at Δφr =90° for experiment C).

Analysis of the two roughness effects 580
The two processes introduced by surface roughness are decoupled using RSRT to better characterise roughness effects as a function of snow properties and illumination conditions.

Effective angle effect
To simulate the effective angle effect, we count all photons that have not been absorbed after the first hit. RSRT is run at 1000 nm with a KZ04 configuration, and we sum the total upward and downward intensity considering one hit only. 585 Simulations are performed for various θs and Δφr. The initial conditions of the A-η27% experiment without slope and under direct sunlight are used. Roughness shapes are rectangular and the SSA is low (7.4 m 2 kg -1 ), which lead to a maximal effect of incident angle variations. The Lambertian configuration yields a constant albedo, as expected (no angle dependence). The albedo decreases of 0.04 due 590 to shadow areas that receive less radiations. With the KZ04 configuration, Fig. 11a and 11b shows that the effective angle effect is strongly linked to illumination conditions. Firstly, when θs > 50°, the model predicts a strong drop in albedo (Fig.   11a). With rectangular roughness shapes, the local incident angle of photons hitting the vertical sides facing the sun is lower than that of a smooth surface when θs > 45° if Δφr = 90°. Thus, photons penetrate deeply in the snowpack before being eventually redirected upward, which conduces to a stronger decrease in albedo relative to a smooth surface. Conversely, when 595 θs < 50°, the effective incident angle is higher over roughness sides facing the sun compared to that of a smooth surface. It leads to an increase in albedo, and this is why Δα is higher with the KZ04 configuration than with the Lambertian configuration when θs < 50° (Fig. 11a). Hence, the reduction in albedo depends on the slope of roughness sides (i.e. their shapes). Fig. 11a also illustrates that in presence of roughness, albedo decreases more rapidly with θs at large values of θs. Therefore, surface roughness play a more important role at grazing angle (large θs). Moreover, our results show that the effects of roughness 600 becomes negligible at 1000 nm when θs < 30°. The albedo decrease caused by the effective angle effect only is of 0.04 for experiment A-η27%, when θs = 63° (Fig. 11a, ). Secondly, by changing the incidence angle, the roughness orientation also plays an important role (Fig. 11b). The reduction in albedo caused by the effective angle effect goes from 0 when Δφr = 0° to 0.09 when Δφr = 90° for experiment A.
To sum up, the albedo decrease due to the effective angle effect becomes rapidly stronger with θs at large θs (θs > 50°) and 605 when Δφr = 90°. In experiment A, the model predicts a decrease in albedo of 0.7 when θs = 80°, caused by the effective angle effect only, i.e a drop 75% stronger compared to that of θs = 63°. Therefore, it is necessary to account for the intrinsic BRDF of the snow to simulate realistic albedo over rough surfaces, in particular in Polar Regions where θs is high.

Multiple reflections
RSRT is run by varying SSA values with the KZ04 configuration and under diffuse sunlight to simulate the trapping effect 615 of photons only (for the A-η27% experiment, Sect. 3.3). Simulations are performed over a smooth and a rough surface to compute Δαsim,rough . Results are shown in Figure 12 as a function of SSA. The impact of multiple reflections is higher for SSA between 8 m 2 kg -1 and 14 m 2 kg -1 , with a maximum effect at SSA = 9 m 2 kg -1 . For the experiment A-η27%, the measured SSA is of 7.4 m 2 kg -1 , and it induces a simulated albedo equal to 0.6 at 1000 nm. Fig. 12 shows that at SSA = 7.4 m 2 kg -1 , Δαsim,rough decreases of 0.035 with multiple reflections, which is significant. The impact of multiple reflections is larger for intermediate 620 values of albedo since photons have the same probability to be absorbed or reflected at each collision. Fig. 12  that multiple reflections are less sensitive at large SSA, as discussed in Sect. 4.2. Hence, it leads to albedo close to 1 and the absorption is too low to trap the photons. Similar results were found in the literature (O'Rawe, 1991;Warren et al., 1998).
Therefore, for the experiment A-η27%, we predict that albedo decreases by 0.04 with multiple reflections and by 0.04 due 625 to a change of the incidence angle, i.e. a total albedo decrease of 0.08 due to the presence of surface roughness. Effective angle effects increase with large θs and low SSA, while the impact of multiple reflections becomes larger when SSA correspond to intermediate value of albedo in the near-infrared wavelengths. Both effects are stronger when the roughness orientation is perpendicular to the sun. Figure 12. Δαsim,rough variations as a function of SSA (m 2 kg -1 ). RSRT simulations are computed with the KZ04 configuration at λ=1000 nm, with the initial conditions of the experiment A-η27% (rectangular shapes, θs = 63°, Δφr =71°, without slope, Table 1). The vertical dotted lines indicates the measured SSA (7.4 m 2 kg -1 ).

Impact on the radiative balance
In this study, the observed albedo change due to the presence of surface roughness may seem low, of the order of a few 635 percent. The decrease in albedo impacts the radiative balance by increasing the proportion of absorbed energy, estimated with the net short wave radiation (SWnet). To illustrate the importance of such an albedo decrease on the radiative balance, we compute SWnet using RSRT with the simple approach described in the following. The net short wave radiation in the 0.35 -4 µm range (in W m -2 ) is estimated with Eq. (17) where ( , ) and ( ) are the direct and diffuse albedo, and dir ( )and ( )the direct and diffuse solar spectral irradiance (W m -2 µm -1 ) computed with the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART, Ricchiazzi et al., 1998). SBDART is an atmospheric model computing radiative transfer within the Earth's atmosphere and at the surface, in clear-sky (direct illumination) and cloudy conditions (diffuse illumination).
The net short wave radiation is estimated with Eq. (17) over the C smooth and C rough 90° surfaces using αsim,smooth and 645 αsim,rough, respectively, at θs = 68°. SBDART is run with a mid-latitude winter atmospheric profile, at 1729 meters high (elevation of the site of experiment C), and at noon. Between these smooth and rough surfaces, the decrease in broadband albedo is equal to 0.05, resulting to an increase of the SWnet from 15 W m -2 to 27 W m -2 . Thus, the energy absorbed by the snowpack increases by almost a factor two (+80 %) with the presence of roughness. These results illustrate the necessity to consider surface https://doi.org/10.5194/tc-2019-179 Preprint. Discussion started: 30 August 2019 c Author(s) 2019. CC BY 4.0 License. roughness in the estimation of the surface energy budget. The RSRT was evaluated with artificial roughness here, and further 650 work will logically concern natural rough surfaces.

Summary and perspectives
Four controlled experiments using artificial roughness fields with various geometrical characteristics (fraction of roughness features, orientation, etc.) were studied. Our observations show that the presence of macroscopic surface roughness significantly decrease snow albedo. More specifically: 655 -Even a low fraction of roughness features (η = 7 %) causes a detectable albedo decrease up to 0.02 at 1000 nm relative to a smooth surface, -For higher fractions (η = 27 % and 63 %), and when the roughness orientation is perpendicular to the sun, the decrease ranges between 0.03 -0.05 at 700 nm and of 0.07 -0.10 at 1000 nm. The impact is 20% lower when the orientation is parallel to the sun. 660 -At low SSA (10 m 2 kg -1 ), the albedo sensitivity to surface roughness is twice as large at 1000 nm (NIR) than at 700 nm (visible) due to the higher intrinsic absorption of the snow.
We developed a new model to account for surface roughness in snow albedo simulations. RSRT considers both the 3-D geometric effects introduced by roughness and snow optical properties using a Monte Carlo photon transport algorithm. By considering roughness, albedo simulations are improved by a factor 2 compared to those assuming a smooth surface (RMSD 665 of 0.03 at 700 nm and 0.04 at 1000 nm).
Using RSRT, we analysed how the contributions usually affecting albedo interact with the effects of roughness.
Firstly, we investigated the impact of SSA and slope uncertainties in our roughness analysis. The amplitude of roughness effects is insensitive to SSA variations at high SSA. On the contrary, at low SSA, a SSA decrease of 50 % induces the same reduction in albedo that the one due to the presence of roughness. Hence, the albedo decrease due to the presence of roughness 670 is drastically accentuated when SSA is low (< 10 m 2 kg -1 ) and when the roughness orientation is perpendicular to the sun. This is explained by 1) when the sun elevation is low, the reduction of the local incident angle of roughness sides facing the sun has more consequences on the albedo when SSA is low (higher absorption of photons), and 2) the impact of multiple reflections is larger when the probability of a photon to be absorbed or reflected is well balanced, which is mainly controlled by a low SSA in the NIR (albedo ~ 0.6). In addition, the overall slope of the rough surface changes the local incident angle and 675 accentuates roughness effects when the surface aspect is facing away the sun. Therefore, to accurately quantify the effects of roughness, it is necessary to know SSA variations when albedo measurements are acquired and the slope of the surface.
Secondly, the two processes governing roughness effects were quantified separately with RSRT. We showed that the albedo decrease due to the effective angle effect becomes rapidly stronger with θs at large θs (θs > 50°) and when Δφr = 90°. For instance, the effective angle effect may cause a reduction in albedo 75% stronger when θs goes from 63° to 80° for roughness 680 shapes considered here. The impact of multiple reflections is larger for SSA between 8 m 2 kg -1 and 14 m 2 kg -1 . Thus, the impact of roughness is strongly linked to SSA, slope, the solar zenith angle and the roughness orientation. RSRT provides a useful tool to better characterize the albedo sensitivity to macroscopic surface roughness.
Roughness effects are significant and many biases are introduced by neglecting these contributions. For approaches considering a smooth surface and using simulated and observed albedo to retrieve SSA, the presence of roughness causes a 685 strong underestimation of SSA, which can be of the order of 20 % for roughness features perpendicular to the sun. Moreover, the albedo decrease leads to an increase of the absorbed energy in the snowpack. In one of our experience, we found that a decrease of the broadband albedo of 0.05 causes +80 % of additional net short wave radiations relative to a smooth surface.
This result highlights the necessity to take into account the roughness effects to compute the surface energy budget. RSRT was https://doi.org/10.5194/tc-2019-179 Preprint. Discussion started: 30 August 2019 c Author(s) 2019. CC BY 4.0 License.