Modelling surface temperature and radiation budget of snow-covered complex terrain

The surface temperature controls the temporal evolution of the snowpack playing a key role in many physical processes such as metamorphism, snowmelt, etc. It shows large spatial variations in mountainous areas because the surface energy budget is affected by specific radiative processes that occur due to the topography, such as the modulation of the irradiance by the local slope, the shadows and the re-illumination of the surface from surrounding slopes. These topographic effects are often neglected in large scale models considering the surface as flat and smooth. Here we aim at estimating the 5 surface temperature and the energy budget of snow-covered complex terrains, in order to evaluate the relative importance of the different processes that control the spatial variations. For this, a modelling chain is implemented to derive surface temperature in a kilometre-wide area from local radiometric and meteorological measurements at a single station. The main component is the Rough Surface Ray-Tracing model (RSRT), based on a photon transport Monte Carlo algorithm to quantify the incident and reflected radiation on every facet of a mesh, describing the snow-covered surface. RSRT is coupled to a 10 surface scheme in order to estimate the complete energy budget from which the surface temperature is solved. To assess the modelling chain performance, we use in situ measurements of surface temperature and satellite thermal observations (TIRS sensor aboard Landsat-8) in the Col du Lautaret area, in the French Alps. The satellite images are corrected from atmospheric effects with a single-channel algorithm. The results of the simulations show (i) an agreement between the simulated and observed surface temperature at the station for a diurnal cycle in winter within 0.3 °C; (ii) the spatial variations of surface 15 temperature are on the order of 5 to 10 °C between opposed slope orientations and are well represented by the model; (iii) the importance of the considered topographic effects is up to 1 °C, the most important being the modulation of solar irradiance by the topography, followed by altitudinal variations in air temperature, long-wave thermal emission from surrounding terrain, spectral dependence of snow albedo, and absorption enhancement due to multiple bounces of photons in steep terrain. These results show the necessity of considering the topography to correctly assess the energy budget and the surface temperature of 20 snow-covered complex terrain. 1 https://doi.org/10.5194/tc-2021-180 Preprint. Discussion started: 30 June 2021 c © Author(s) 2021. CC BY 4.0 License.

topography, the large number of possible approximations and the degrees of complexity used at both metre-scale roughness 90 and complex-terrain scales. These effects and approximations are summarized in Table 1 and illustrated in Fig. 1.  In addition, most models cited above consider only broadband SW fluxes, neglecting the spectral dependence of snow albedo and incident radiation. The probable consequence is an inaccurate estimation of the absorption enhancement due to neglecting the large difference of albedo between visible and near-infrared. In the former domain, the albedo is high (typically over 0.95), implying intense multiple scattering but extremely weak absorption. In contrast, in the near-infrared the albedo is lower and 95 closer to the optimal albedo for enhancement (0.5), where multiple scattering and absorption are balanced (Warren et al., 1998).
To investigate the relative importance of the topographic effects, this study aims at estimating the snow surface temperature in a mountainous area with a modelling chain that uses local in situ radiometric and meteorological measurements from a single station. The chain includes a 3D radiative transfer model, the Rough Surface Ray Tracer (RSRT) model , 100 which launches a set of photons to the snow surface described by a triangular mesh, i.e. a connected set of triangular facets, that can be derived easily from a Digital Elevation Model (DEM). Its spatial resolution limits the scope of this study to topography at the decametre to kilometre scale range. The simulation results are used in a surface scheme (RoughSEB model) to compute the short-wave and long-wave net radiation and the turbulent heat fluxes on each facet. Eventually the surface temperature is deduced. Satellite thermal infrared observations from Landsat-8 are used in order to evaluate the modelled spatial variations.
105 Table 1. Several effects relevant to short-wave (SW), long-wave (LW) and turbulent heat fluxes calculation in complex terrains and rough surfaces. Other effects such as those involving the atmosphere are beyond the scope of this study.

Effect
Spectral  (Arya, 1988): (i) the net radiation fluxes, which are split into the contributions of the short-wave radiation from 0.3 µm to 2 µm (SW net ) and the long-115 wave radiation from 2 µm to 100 µm (LW net ); (ii) the sensible heat flux (H), which measures the exchange of heat between the surface and the air just above; (iii) the latent heat flux (L) resulting from water state changes (sublimation or condensation) at the surface and exchange with the air above, and (iv) the ground heat flux (G), which is transferred to the snowpack and lead to a change of temperature or melt. Here, we consider only temperature estimation and assume below freezing temperature so that melt is not involved. The sum of these fluxes is null according to the first principle of the thermodynamics, taking into account that the surface has no internal enthalpy: In the RoughSEB model, all these terms are calculated for each facet (hereafter noted with the subscript f) of the modelled snow surface. The computation of the SW net radiation (Sect. 2.1.1) involves: (i) the Santa Barbara DISORT Atmospheric Radiative Transfer model (SBDART - Ricchiazzi et al. (1998)), an atmospheric model that simulates the incoming spectral 125 irradiance above the studied area (top of mountains -TOM), and (ii) the RSRT model , based on a Monte Carlo photon tracking algorithm that computes the path of the photons launched towards the surface and allows considering the modulation of SW radiation by the terrain slope and aspect. The simulations are run in both direct and diffuse illumination conditions (noted with subscripts dir and diff), and the atmospheric effects (i.e. atmospheric attenuation) are neglected within the studied area (between the surface and TOM). The scene is also considered to be covered by pristine, pure snow (i.e. no 130 impurities), which is applicable in winter. The calculation of the LW net radiation (Sect. 2.1.2) needs the downwelling LW flux, measured by a local station representative of the area, and subsequently updated by accounting for the reduced sky-view factor and the thermal emission of surrounding terrain as in Arnold et al. (2006). In order to consider the topographic effects affecting the SW radiation and the spectral dependence of snow albedo and incident 140 radiation, the calculation of the net SW radiation needs to be carefully performed. The aim is to compute it on every facet f due to the direct and diffuse irradiation: where the direct and diffuse absorbed SW radiation by every facet are computed using an absorption coefficient A dir, diff, f , derived from the RSRT model, and the spectral irradiance coming from the sky (I dir, diff (λ)), issued from the SBDART model.
The integration limits (λ 0 and λ 1 ) are respectively equal to 300 nm and 2000 nm.
A naive way to compute these terms would be to use the ray-tracing RSRT model to trace many photons for every wavelength in the solar spectrum to compute the absorbed SW radiation by each facet. Nevertheless, with millions of facets and hundreds 150 of wavelengths, this approach would imply an enormous computational cost, due to the inefficiency of the Monte Carlo ray tracing method. To overcome this, an alternative approach consists of taking advantage of the result of the RSRT model to account for multiple bouncing provided a few assumptions. The RSRT model can indeed compute the number of times a photon has hit a given facet regardless of the albedo (and so of the wavelength), according to the bounce order of the photon (first reflection, second reflection, ...). Noted n (i) hit, d, f , it corresponds to the proportion of photons that hit the facet on their i th hit 155 (d being dir or diff depending on the illumination conditions -direct or diffuse).
Assuming the area has an uniform albedo (same snow properties) and the reflection is Lambertian, the absorption coefficient is computed by taking into account (i) the spectral dependence of snow albedo, and (ii) the fact that the illumination received by each facet depends on the cosine of the local solar zenith angle and the absorption enhancement produced by multiple bouncing, with: where θ s is the local solar zenith angle and α dir, diff (λ, θ s ) is the snow spectral albedo in both direct and diffuse illumination.
Here, it is computed using the Asymptotic Radiative Transfer (ART) theory (Kokhanovsky and Zege, 2004). Its expression is presented in the Appendix A1, provided several assumptions about the snowpack (semi-infinite, vertical and horizontal homogeneous layers) and the surface (flat and smooth -facets are small enough to be considered as it). It considers the snow specific surface area (SSA), which is needed as input in the modelling chain.
The net broadband SW radiation per facet is therefore calculated by means of Eqs. (2, 5, 6), previously accounting for the spectral irradiance coming from the sky and integrating each component over the 300 -2000 nm wavelength range. that is then used in the second step to account for the emission of surrounding slopes by computing the average upwelling 175 long-wave from each facet as in Arnold et al. (2006), and using: where V f is the sky-view factor estimated with RSRT. V f is different for each facet so that facets in the valley receive more energy from the surrounding slopes than facets at the summits of the mountains. It is indeed equal to the proportion of the launched photons hitting a facet on their first bounce in diffuse illumination, namely V f = n For the upwelling long-wave radiation, LW u, f is determined by the Stefan-Boltzmann law: with = 0.98, σ the Stephan-Boltzmann constant and T s the snow surface temperature.

185
While the main focus of this work is on the role of the topographic effects controlling the radiative budget of the surface, the turbulent heat fluxes needs also to be assessed to compute the energy budget. We follow a very simple modelling approach at this stage, potential improvements are let to further work in the future. The sensible and latent heat fluxes are calculated following the one-level approach: where ρ air and c p,air are the density and heat capacity of the air, U is the wind speed, T air is the air temperature, L s is the sublimation heat, q sat (T s , P s ) is the specific humidity at snow surface temperature T s and pressure P s , and C H is a surface exchange coefficient. In principle this coefficient depends on atmospheric stability but for the sake of simplicity, a neutral situation is considered here. C H therefore depends only on the aerodynamic roughness length z 0 , assumed constant across the 195 scene and equal to 10 −3 m following previous works (Brock et al., 2006). The expression of this coefficient is provided in the Appendix A2, and the values of the symbols defined here are presented in the Table B1. The air temperature and wind speed are given from a meteorological station in the scene (Sect. 2.3). To account for the differences in altitude within the scene, the lapse rate Γ is taken into account: where T air, f is the air temperature over the facet f . We choose Γ = -6.5°C km −1 , the environmental lapse rate as defined in the International Standard Atmosphere (ISO 2533(ISO :1975. Wind speed and relative humidity are however considered constant.
Note that the specific humidity being deduced from the relative humidity and air temperature, it indirectly depends on altitude.
The ground heat flux (G) is here neglected as the snowpack is considered thermalized, meaning that no energy is exchanged with the snowpack. This assumption is checked in the discussion (Sect. 4.2) with regard to our results.

Snow surface temperature estimation
The equation of the surface energy budget (Eq. 1) is to be solved for T s for each facet f . However this equation is non-linear due to the T 4 s term (upwelling LW radiation) and the complex dependence of the surface specific humidity to T s . Solving such an equation for millions of facets can be computationally intensive. To avoid this issue, we linearize the specific humidity at saturation: 210 q sat (T s , P s ) = q sat (T air , P s ) + (T s − T air ) · q sat (T air , P s ) where: q sat (T air , P s ) = q sat (T air , P s ) − q sat (T air − ∆T, P s ) ∆T with a ∆T = 5 K. The simplified surface energy budget equation eventually become a quartic equation for T s of the form: where a, d, e are a constant (for each facet). The general solution is shown in the Appendix A3 for completeness. The evaluation of this equation for a large number of facets is computationally very efficient compared to using a non-linear optimization technique.

Description of the existing models
The modelling chain involves two existing models, namely the SBDART model (Ricchiazzi et al., 1998) and the RSRT model 220 . Here we briefly describe the parameters to run them in the present work.
The SBDART model provides the short-wave irradiance spectrum at the top of the mountains. It considers several atmospheric transmission models and the Mie theory to take scattering into account. The simulations are run in the spectral range 300 -2000 nm, by steps of 3 nm. In our case, a generic mid-latitude winter atmospheric profile is assumed and the aerosol optical depth at 550 nm is 0.08 ("rural" profile in SBDART). The elevation is 2052 m.a.s.l. (see Sect. 2.3) and solar zenith and 225 azimuth angles corresponding to the desired date and time.
RSRT model provides the number of photon hits on each facet according to its bounce order, noted n (i) hit, d, f . In this model, each launched photon has an initial propagation direction i, specified by the zenith and azimuth angles, and a random initial position over the mesh. Edge effects are avoided by excluding the outermost 15% of the mesh. The model calculates the propagation with the following main steps: (1) estimation of the next intersection between the photon path and the mesh; (2) determination 230 of the reflection direction (random cosine-weighted hemispherical distribution in our case to simulate Lambertian reflection), and (3) update of the direction i. The algorithm iterates until the photon escapes from the scene or up to a maximum number of bounces of 15. The number of facets in the mesh is typically 10 6 and the number of photons launched is 4·10 8 .  (EPSG: 2154), and the original spatial resolution was of 5 m. It was however resampled to 10 m due to computational limitations. To create the mesh, the centre of each pixel is a vertex that is connected to its closest neighbours, which eventually leads to two triangular facets shared by the same vertex.

Study area and in situ measurements
The study area also includes the measurement station FluxAlp (45.0413°N, 6.4106°E), on the Pré des Charmasses site.
This station collects meteorological and radiometric observations, and its description can be found in Dumont et al. (2017).

245
Radiation fluxes are measured with a Kipp & Zonen CNR4 radiometer. The data needed as input for the modelling chain is issued from here (meteorological -air temperature, wind speed, relative humidity and radiometric -downwelling long-wave radiation flux). The latter is assumed to be equal to that received above the area (TOM), as the LW re-illumination at the station is very low according to initial tests of the modelling chain. FluxAlp measurements are averaged over 30 min, i.e. a measurement at 10h30 UTC is the average of the period 10h00 UTC -10h30 UTC. The upwelling LW measurements at 250 FluxAlp are used to assess the modelled surface temperature temporal variations in a single point. They are complemented by a spatially distributed surface temperature dataset described in the following section.
In addition to the automatic measurements, manual measurements of SSA for two consecutive winter seasons (2016 / 2017 and 2017 / 2018) have been collected occasionally . The measurements were collected with the DUFISSS instrument (Gallet et al., 2009)

Surface temperature retrieval with Landsat-8 observations
Spatial variations of surface temperature are retrieved from satellite observations. The two thermal bands (TIRS -Bands 10 and 11) onboard Landsat-8 cover the spectrum between 10.6 µm to 12.51 µm, with a spatial resolution of 100 m (resampled 260 by Cubic Convolution methods to 30 m) and a 16 day repeat cycle. Different methods to correct atmospheric effects have been implemented to retrieve Land Surface Temperature (LST, hereafter), based on split-window methods (Jin et al., 2015), mono-window techniques (Tardy et al., 2016), or a single-channel approach (Jiménez-Muñoz and Sobrino, 2003). Soon after the launch of Landsat-8, stray light was observed on thermal data (Montanaro et al., 2014), coming from scattering of outer radiance in Band 11. Methods based on only one band are therefore suggested, so here we apply a single-channel approach, 265 which consists of approximating the atmospheric functions from atmospheric water vapour content (w, in g cm −2 ). Cristóbal et al. (2018) presented an improved single-channel method dependent not only on water vapour content, but also on nearsurface air temperature (T a ), which are available from reanalysis data. The recently available Landsat Collection 2 Surface Temperature product is also considered here and compared to the results of the aforementioned LST retrieval methods in our particular case of study. Figure 4 shows the flowchart to retrieve LST from satellite observations, following Cristóbal et al. Both the single-channel method (SC method -Jiménez-Muñoz and Sobrino (2003)) and the improved single-channel method (iSC method -Cristóbal et al. (2018)) have been implemented in this study. LST is calculated for each pixel by applying the radiative transfer equation to a sensor channel, which eventually leads to:

Results
The spatially-resolved LST observations from Landsat-8 are first assessed in the study area, before the evaluation of the model 285 simulations against the observations and the local measurements. 3.1 Surface temperature observations from Landsat-8 The surface temperature observations from Landsat-8 (LST) are compared to in situ FluxAlp measurements (Fig. 5a). The surface temperature from Landsat-8 is extracted from the pixel covering the location of FluxAlp from all 20 thermal images.
Both the single-channel (SC) method and the improved single-channel (iSC) method show better estimations at FluxAlp station 290 than the official Collection 2 Surface Temperature product (-3.5°C underestimation). As a result, we excluded this latter dataset from further analysis. The iSC method provides generally higher temperatures than the SC method, with a mean bias of -1.3°C and -2.0°C, respectively. Its total error is of 2.0°C RMS, dominated by the bias (2.6°C for the SC method). The improved single-channel method shows slightly more accurate results, and it is used in the following to evaluate the estimation of snow surface temperature by the modelling chain. The modelling chain to estimate T s is evaluated over a diurnal cycle at the FluxAlp station (Fig. 6).    The simulations are colder in general, with a negative bias comprised principally between -4°C and 0°C. Some simulations are however warmer than the satellite observations. The standard deviation of the difference varies mostly between 2 and 4°C. be explained with an early onset of snowmelt (snow patches in the lowest areas), a particular situation that is not yet correctly resolved by the model. Nevertheless, considering all acquisitions, these differences don't seem to be clearly related to the time of the year (i.e. early or late winter), so strong conclusions cannot be drawn.
This results show the performance of the RoughSEB model to simulate the snow surface temperature and the surface energy budget in complex terrain within a reasonable accuracy. Their temporal and spatial variations are also well represented in the 345 study area. To understand these variations, the role of the topographic effects that govern them is addressed in the following section.

Relative importance of the topographic effects
In order to quantify the relative importance of the topographic effects, we have run a series of simulations where every effect considered in this work is disabled at a time. Figure 9 shows the impact on the surface temperature distribution across the area   When no topography is taken into account (i.e., a perfectly flat surface instead of the actual topography), the snow surface temperature is uniform (-6.9°C) across the scene. This value overestimates the mean temperature of the reference simulation 360 (-7.9°C) showing that not only the variations of temperature are not represented but even on average a simulation on a flat terrain is not equivalent to the mean temperature of the area with topography. This basic simulation highlights the considerable effect of the topography on the surface temperature and the importance to take it into account even for large scale simulations.
Considering all the 20 dates, neglecting the topography results in an overall overestimation of 1.0°C (median value) and the standard variations with respect to the reference simulation are high with a median value of 2.5°C and a maximum of 5°C.

365
The simulation with single-scattering only in the SW, (only the first bounce of the photons is considered by RSRT, multiple scattering is neglected) shows small differences with respect to the reference simulation. It is only slightly colder, in particular the coldest pixels. This is mainly due to neglecting the re-illumination caused by multiple bouncing in the shadowed (and cold) areas. Nevertheless, the impact is not significant, the mean difference being of barely different from zero (median of -0.2°C, standard variations of 0.1°C). The smallest impact is observed in December which could be explained by a dependence on the 370 solar zenith angle, and on how the topography modulates the received solar irradiance. Overall, this result means that, at least in our study area, the absorption enhancement caused by multiple bouncing is almost negligible.
When neglecting the spectral dependence of snow albedo and thus considering only the broadband albedo (≈ 0.87 on 18 February 2018), the simulated surface temperature is slightly warmer than the reference with overall differences of 0.4°C (median of the means) and of 0.6°C (median of standard deviations) and a maximum beyond 1°C. This effect involving the 375 coupling between spectral dependence of snow albedo and topographic effect is important to take into account, as it plays a role that needs to be considered.
The effect of the thermal emission by surrounding terrain in the LW is obtained by estimating the surface temperature as if the terrain was flat (first step as detailed in Sect. 2.1.2). The peak of the distribution is less marked and widened in the range (-5°C --8°C), and the impact is a systematic underestimation of T s . This is true for the warmest areas of the scene (higher 380 than -5°C) where the underestimation is of -0.7°C, and in particularly true for the coldest areas (lower than -12°C), with an underestimation of -1.2°C. In the latter, direct radiation is absent and the radiative budget is dominated by diffuse and weak illumination. The thermal emission by surrounding faces warms these cold areas as a function of their sky-view factor. The mean difference goes down to -0.6°C (standard variations of 0.4°C) when considering all the scenes.
With respect to the altitudinal variations of air temperature (lapse rate effect), the distribution shape squeezes and even 385 becomes bimodal, with two marked peaks at -5°C and -10°C. The difference to the reference simulation is significant when considering all the dates, with a median value of 0.9°C and a median standard deviation of 1.3°C. Introducing the lapse rate effect tends to warm the air on the lower parts of the scene (usually the warmer) and the opposite on the higher and colder facets. Since the FluxAlp measurement station (where the reference air temperature is taken) is in the lower range of altitudes of the study area the overall impact neglecting the lapse rate effect is an overestimation (0.9°C). This result is specific to 390 our setting and using a different reference air temperature would change this result. In principle it is possible to choose the reference at the mean altitude of the area which leads to a null bias. Nevertheless, the standard variations are also significant and this is not specific to our setting. It effectively shows that neglecting the altitudinal variations of air temperature results in an overall significant error of 1.3°C (median) in surface temperature over the study area.
Theses results gives the relative importance of the topographic effects investigated here. Neglecting the topography (i.e. a 395 flat surface) is, as expected, the most important source of error when simulating snow surface temperature. Neglecting the altitudinal variations of air temperature (lapse rate effect) is the second effect in terms of importance. It is followed by the thermal emission by surrounding terrain in the LW and the spectral dependence of snow albedo, the latter being slightly less important. Finally, the absorption enhancement caused by multiple bouncing is almost negligible.

Retrieval of surface temperature from satellite observations
The assessment of our modelling chain was performed against in situ and satellite observations, the latter being crucial for the spatial variations. However, these depend on the choice of the method for the atmospheric correction. Here we implemented two single-channel methods for Landsat-8 thermal imagery: the SC method (Jiménez-Muñoz and Sobrino, 2003), and the iSC method (Cristóbal et al., 2018) and we compared them to the recent Surface Temperature product. The iSC method is 405 the most accurate (within ≈ 1°C) at the FluxAlp measurement station. The Collection 2 Surface Temperature product is the less accurate (within 3.5°C). However, this does not preclude of the quality in the whole area. In particular, we assumed that the whole scene is covered by snow, with an emissivity equal to 0.98. However, in alpine areas, vertical mountain ridges or patches without snow on sun-facing slopes are frequent and may results in more variable emissivity over the scene. A possible future improvement would be to consider an emissivity mask, where each pixel would have a particular value depending on 410 the presence of snow, rocks, grass, etc. This is normally achieved by means of NDVI-based classifications (Li et al., 2013).
The Surface Temperature product (which is already generated using this method) presents worse results at our validation point, but overall, the differences across the study area between the iSC method and this official product are of 0.3°C (median of standard deviations) considering the whole dataset, which is considerably lower than the difference at FluxAlp. The satellite thermal observations from Landsat-8 are therefore correct enough to assess the model performance, with an accuracy of the 415 order of 1°C and a precision lower than 0.5°C.

Snow surface temperature estimation
One question that arises is about the performance of the model to estimate the T s and its temporal and spatial variations in complex terrain. The simulations show an overall agreement with measurements and satellite observations (Sect. 3.2). The simulated fluxes to be considered in the surface energy budget and the temporal evolution of T s are well represented for a 420 daily cycle (Fig. 6). The net SW flux is slightly underestimated during most of the day (except in the early morning when it is slightly overestimated), but the impact on surface temperature is balanced by the other terms of the energy budget, in particular due to the slight overestimation of the net LW flux. The turbulent heat fluxes are essential for an accurate calculation of the surface energy balance, in particular during the night, when they balance the LW radiation flux. Their treatment in the model is very simple, the atmosphere is assumed neutral, the aerodynamic roughness length is uniform, the wind is uniform, etc. Some 425 of these simplifications are required to achieve good computational performances and may be challenging to take into account (wind field) while others (e.g. atmospheric stability) could be easily implemented (Essery and Etchevers, 2004).
The ground heat flux was also neglected meaning that the surface temperature reacts immediately to changes of the downwelling radiation. In reality, the snowpack has some thermal inertia which delays and moderates the diurnal amplitude of T s , resulting in an overestimation in the morning and underestimation in the late afternoon, when the cooling of the surface is 430 more pronounced in the simulations as the solar zenith angle increases. Nevertheless, this delay is barely visible in our case ( Fig. 6). Snow is indeed a highly insulating medium and has a small thermal inertia. With a thermal conductivity of around 0.2 Wm −1 K −1 (Sturm et al., 1997), the daily wave penetrates by no more than 20 cm in the snowpack. It means that with a diurnal cycle half-amplitude of ≈ 7°C in surface temperature, the maximum temperature gradient in the upper snowpack is of the order of 35 Km −1 . According to the Fourier law, this implies a maximum ground heat flux G ≈ 7 Wm −2 which is an order 435 of magnitude lower than the radiative and turbulent heat fluxes estimated in our case (Fig.6). This simple calculation confirms that neglecting G is acceptable in a first approximation when snow is relatively fresh and is not melting. In spring, with a denser and less insulating snowpack and with melt, this approximation needs to be reconsidered. However, taking into account the thermal diffusion in the snowpack over millions of facets represent a very significant computational cost and requires approximations.

440
The choice of input parameters such as SSA and aerodynamic roughness length (z 0 ) is critical for the simulations. Local measurements are not always available and may not be representative of the whole area. SSA plays a crucial role on the albedo and so on the short-wave absorption by the snowpack. The sensitivity of snow surface temperature to SSA is shown in Fig.   11 (top), for the diurnal cycle presented in Fig. 6. SSA varies in ± 10 m 2 kg −1 with respect to the reference simulation for the whole time series (20 m 2 kg −1 ). The impact of varying SSA on T s is up to 1°C, and is larger when considering a low 445 SSA than a large SSA. This is mainly due to the fact that the relationship SSA-albedo is not linear, as shown by Domine et al. (2006). In general spatio-temporal variations of SSA are expected in our area because the initial SSA value depends on snowfall temperature mainly and its subsequent evolution is a general decrease depending on the thermal evolution of the snowpack (Domine et al., 2008). SSA is expected to be higher at higher altitudes and in the shadows where the conditions are colder, and so does the short-wave absorption which would tend to increase the cooling in these areas compared to lower SSA 450 areas. This could also explain the differences between the simulated and measured net short-wave flux, as here the value was kept constant for the whole time series.
The aerodynamic roughness length controls the sensible and latent heat fluxes, so it shall also be carefully tuned. In this study, a value of 10 −3 m is assumed, which is standard according to previous works (e.g. Brock et al. (2006)). The sensitivity of T s to aerodynamic roughness length is shown in Fig. 11 (bottom). z 0 varies by a factor of 5 with respect to the reference 455 simulation. The choice of z 0 has a more significant impact on the simulated surface temperature, in particular during the night (up to 2°C), when the SW radiation is absent. Martin and Lejeune (1998) found similar variations using the snow model Crocus, with a slightly different approach that considered changes of atmosphere stability, as well as Kuipers Munneke et al.
(2009) with even another approach for the turbulent heat fluxes. According to Brock et al. (2006), snow z 0 can vary up to three orders of magnitude depending on the time of the season (i.e. early or mid ablation season) and also on snow type (i.e. fresh 460 snow or melting snow). To estimate its expected spatial variations in a mountainous area of several km 2 , a parametrization of SSA evolution would be possible (Domine et al., 2007) or by using SSA retrieved by satellite (e.g. Kokhanovsky et al. (2019)).
The spatial variations of surface temperature are clearly dependent on the topography and are correctly simulated (Fig.   7). They seem to depend in particular on slope orientation, as shown in Fig. 12. For larger slopes (> 30°), the lack of direct radiation governs the surface temperature in some areas, such as the northern (N-W to N-E) slopes covered in the shadow at the 465 southern part of the scene. The south-facing slopes are more numerous and feature an extended range of surface temperature.
The temperature difference is large between opposed slopes (on the order of 5 to 10°C in a few hundred meters). Overall, the mean T s of the northern facets is -10.5°C (± 1.7°C). They are considerably colder than the southern (S-E to S-W) ones (-6.7°C (± 1.5°C)). These differences are consistent with previous studies (e.g. Fierz et al. (2003)), and highlight the necessity of accounting for spatial variations of surface temperature in mountainous areas, where larger slopes prevail over gentle terrain.
470 Figure 12. Distribution of simulated surface temperature as a function of the aspect for slopes larger than 30°. The radius of each of the 16 sectors corresponds to the normed (displayed in percent) quantity of facets.

Role of the topographic effects
The results presented in this study show that the topography controls the energy budget and the surface temperature to a large extent. However modelling all the processes involved in complex terrain may lead to prohibitive computational cost for most applications. This has an adverse outcome: some of them are usually neglected or approximated. Our modelling chain takes into account a relatively comprehensive set of processes, especially on the radiative aspects. The motivation of this study is graphic effects. This agrees with the results presented by Arnold et al. (2006), showing a similar ratio of importance between the role of the shadows and the LW contribution when considering glacier melt. To take into account this effect, the downwelling LW radiation flux in each facet is updated with the thermal emission of the surrounding facets, which is derived with an average T s for the whole scene. Such an approximation saves a lot of computation time but, as a result, the differences in thermal emission due to spatial variations of surface temperature are masked. The warmest facets will emit more LW radiation 490 than the coldest ones, leading to differences in the modelled T s .
Our results show a significant impact on the simulated T s by disabling the altitudinal variations of air temperature, with a warming effect up to 1°C. The spatial and temporal distribution of air temperature, as well as wind dynamics in mountainous areas, have been widely investigated (Jiménez and Dudhia, 2012;Rotach et al., 2015), and different parameterization and estimation methods exist to overcome its complexity (e.g. Wood et al. (2001)). The approach implemented here is simple, with 495 several approximations and important assumptions as a constant aerodynamic roughness length and wind speed across the study area. Future improvements of our approach could benefit from downscaling methods to better represent wind distribution over complex terrain, as in Helbig et al. (2017). Arnold et al. (2006) also pointed out the role of the anisotropic reflectance of snow and ice at high solar zenith angles. This is found to be essential to correctly estimate the radiative budget. The RSRT model can account for this effect (Larue et al.,500 2020) but we did not consider it here. The reason is that we preferred to account for the spectral dependence of snow albedo, which is an improvement with respect to most of the models that only account for broadband albedo. The impact in T s is however limited, but errors up to ≈ 1°C are possible. As running RSRT at many wavelengths is very expensive, we developed the effective method present in Eqs. 5 and 6. Unfortunately, this method relies on the Lambertian approximation which makes it incompatible with taking into account the reflectance anisotropy.

505
Another topographic effect investigated here is the influence of the multiple bounces of photons between the faces in the short-wave domain, also called re-illumination or absorption enhancement. It is simulated with the RSRT model  and represents the most expensive step of the modelling chain in terms of computational cost, as each photon is tracked until it escapes the scene, after one or several bounces. Our results suggest that the re-illumination can be neglected, with a limited impact on the surface temperature estimations (< 0.3°C). This implies that RSRT is only needed to track the shadows, 510 slope and compute the sky view factor. In principle these effects can be taken into account by other methods (e.g. Dozier et al. (1981)). Nevertheless, other authors have drawn different conclusions regarding the role of re-illumination in the same study area. Lamare et al. (2020) found a significant contribution of multiple reflections on the simulated TOA radiance in rugged terrain. An hypothesis to explain this discrepancy is the dependence on the sun's position and the configuration of the terrain.
Earlier simulations in the winter season (closer to December 1st) show a lower impact than later ones. However, this point 515 requires further investigation to derive stronger conclusions, in particular as atmospheric scattering and absorption effects due to the air are neglected in the present work.

Conclusions
We have investigated the relative importance of several topographic effects in a mountainous terrain in the French Alps.
For this, we have first developed a chain of models to predict surface temperature, by combining existing radiative transfer 520 models (RSRT, SBDART) and a new surface energy budget model (RoughSEB). This chain has been evaluated against in situ measurements and remote sensing thermal observations to account for the spatial variations. The latter have been corrected for atmospheric effects with a single-channel algorithm.
A ≈ 36 h long time series was simulated and an overall agreement is achieved with the in situ measurements, within 0.3°C.
Besides, the bias of the simulations at FluxAlp station for the 20 scenes corresponding to the Landsat-8 acquisitions is only 525 -1.2°C (total error of 1.5°C RMS), which highlights the potential of the chain to simulate surface temperature within a reasonable accuracy and precision. The spatial variations across the 50 km 2 scene are also well represented, with a standard deviation of the differences to the satellite observations comprised between 2 and 3°C, which is small compared to the actual surface temperature variations of 5-10°C between the slopes in the study area.
A few topographic effects that are responsible for such spatio-temporal variations have been investigated and their relative 530 importance is: 1) The modulation of solar irradiance by the terrain slope and aspect (i.e. the presence of topography), 2) the altitudinal changes in air temperature (lapse rate effect), 3) the contribution of thermal radiation coming from surrounding terrain, 4) the spectral dependence of snow albedo, and 5) the absorption enhancement caused by multiple bouncing of photons in the SW domain. Their importance has been quantified and the warming (or cooling) effects are up to 1°C, except for the first one (absence of topography) that lead to errors of several degrees Celsius in surface temperature.

535
The modelling chain shows some limitations justified by the assumptions made on some of the parameters controlling the radiative budget, such as snow SSA or the aerodynamic roughness length z 0 . Nevertheless, it shows good performance to estimate snow surface temperature and its spatial variations, and has several applications. It allows a better understanding of the processes that govern the surface energy budget in snow-covered, mountainous areas. A first extension of the model is to investigate the snowmelt. Another application is the preparation, and in the future the calibration/validation of the thermal 540 infrared TRISHNA satellite mission (Lagouarde et al., 2019) which will provide from 2025 high resolution images (50 m) of the Earth surface temperature in mountainous areas.

2004).
Provided several assumptions about the snowpack (semi-infinite, vertical and horizontal homogeneous layers) and the surface (flat and smooth -facets are small enough to be considered as it), they are expressed as follows: α dir (λ, θ s ) = exp − 12(1 + 2 cos θ s ) 7 where θ s is the solar zenith angle, SSA is the snow specific surface area, ρ ice = 917 kgm −3 is the bulk density of ice at 0°C, 550 γ(λ) is the ice absorption coefficient (from Picard et al. (2016)) and B and g are the shape coefficients of snow, taken from Libois et al. (2014).

A2 Surface exchange coefficient
The surface exchange coefficient involved in the calculation of the turbulent heat fluxes is given by: where z t and z w are the height at which air temperature and wind speed are measured, and z 0 is the aerodynamic roughness length. Their values are provided in Table B1.

A3 General solution to the quartic equation for T s
The simplified surface energy budget equation eventually become a quartic equation for T s of the form: where a, d, e are a constant (for each facet). The general solution can be calculated using: ∆ 0 = 12 a e (A7) A4 Land Surface Temperature retrieval from Landsat 8
The atmospheric functions are statistically fitted from the GAPRI database (Mattar et al., 2015) containing 4714 atmospheric 580 profiles from tropical to arctic atmospheric conditions. The following fit is applied here: ψ i = i w 2 + h T 2 a + g w + f T a + e T 2 a w + d T a w + c T a w 2 + b T 2 a w 2 + a (A16) All the coefficient values (from i to a) are in Cristóbal et al. (2018).