How Forward-Scattering Snow and Terrain Change the Alpine Radiation Balance With Application to Solar Panels

with a ray-tracing model the solar irradiation of mountainous terrain and differentiated single scattered terrain radiation

latter studies the domain reflectance arising from variable surface albedo and rough terrain is anisotropic in nature, the underlying scattering events at the surface are assumed to be Lambertian.In practise, numerous solar radiation models use strongly simplified methods to calculate inter-terrain scattering based on a homogeneous view factor approach (Šúri & Hofierka, 2004;K. Wang et al., 2005).
Although independently of each other the BRDF of snow and the radiative transfer in mountainous terrain has been studied in depth, we are not aware of any study that applies the BRDF of snow to mountainous terrain.To bridge this gap, we present GROUNDEYE, a new model based on a modified radiosity approach, which for the first time includes the combined effect of multiple scattering with the anisotropic forward scattering of snow.Full model validation has been achieved.We provide insight into the interplay of snow-covered terrain and incident sunlight, which leads to complex radiation patterns with local maxima.While snow ultimately reflects a large part of this terrain radiation due to its high albedo, the effect on the radiation balance of efficient absorbers is significantly larger.We therefore investigate how these local maxima can be exploited for photovoltaic production.
A number of studies address the effect of snow on solar panel productivity, but almost all of them address the problem of snow cover on the panels themselves (Andenaes et al., 2018;Andrews & Pearce, 2013;Pawluk et al., 2019).In contrast, Kahl et al. (2019) explicitly studied the positive effect of snow-covered ground on solar panel productivity, but relied on a simple isotropic model of terrain radiation, like other solar energy potential models before it (Hafez et al., 2017;Stanciu & Stanciu, 2014;Yadav & Chandel, 2013).We thus investigate for the first time in detail the role of terrain-reflected radiation in the radiative balance of photovoltaic panels on snow-covered ground.And we demonstrate that the terrain radiation with forward scattering increases the optimal tilt of solar panels in the Alpine area and leads to significantly higher winter irradiance.
The emphasis in this paper is on radiation coming from the terrain, more specifically: reflected by the terrain.This terrain-reflected or ground-reflected radiation is hereafter simply called terrain radiation and is discussed in units of incident flux density 2 [ ] Wm .While the magnitude of the effects depends on the actual topography, the characteristics may in general apply to all snow-covered terrain.

Definitions
This section gives definitions of quantities referred to in the next sections, such as plane albedo, view factor and forward scattering effect.The bidirectional reflectance distribution function (BRDF) characterizes the reflective property of a given surface.Based on Nicodemus (1977) we define the broadband BRDF  1 [ ] sr  i and  v are the zenith angles of the incident and outgoing directions respectively,  i and  v are the corresponding azimuth angles.The syntax used may suggest that  is a purely geometric quantity.But this is not the case, as  also depends on various parameters of the snow cover like grain size, grain shape, surface roughness and dirt (Warren et al., 1998) and on the wavelength of the incident radiation (Schaepman-Strub et al., 2006).By integrating the BRDF over the entire hemisphere, the broadband plane albedo  is obtained, which is the ratio of reflected to incident radiation energy (Nicodemus, 1977): ( , ) ( , , , )sin cos .
View factors are important parameters for calculating the radiation balance of surfaces.However, there is disagreement in literature as to what should be understood by view factor.Generalizing Manners et al. (2012), we define the view factor f as the following function that assigns to a subset  of the visible hemisphere a value between 0 and 1: where  and  are zenith and azimuth angles in the local reference system.With this definition, the view factor  f quantifies the fraction of radiation emitted or received at a given point on a Lambertian surface that passes through the solid angle .The entire hemisphere has a view factor of one.The sky-view factor sky f follows when  is identified with the visible sky.Since the visible hemisphere consists solely of sky and terrain, the terrain-view factor ter f follows directly from sky f : For a given topography, the terrain-view factor , ter p f can be determined for any point p in the terrain.The surrounding terrain, which contributes to , ter p f , can thus be given a terrain-view factor where  p is the solid angle under which terrain is seen from point p.  p f is as defined in Equation 3. We call this quantity  ter f friends' terrain-view factor, because it's an average over the points that p sees all the time.To quantify the effect of anisotropic scattering on terrain radiation, we define the forward scattering effect (FSE) as ,1 1 .

F
is the incident broadband flux density

2
[ ] Wm of single scattered light based on a realistic (anisotropic) BRDF.
,1 Iso ter F is the counterpart, assuming lambertian reflections.Correctly, this should be called anisotropy effect, because the BRDF of snow includes besides the forward scattering peak also other patterns like a darkening at grazing angles (Dumont et al., 2010).However, because the former is the most prominent, we simply call it forward scattering effect.To quantify the effect of multiple scattering on terrain radiation, we define the multiple scattering effect (MSE) as ,1 1 , , therefore the MSE is always equal to or greater than zero.This definition of the MSE refers to anisotropic multiple scattering.Finally, we define the relative difference between anisotropic multiple scattered terrain radiation , and terrain radiation where the first scattering is anisotropic but all further scattering events are isotropic: This rather cumbersome expression quantifies the importance of forward scattering in multiple scattering.
In order to discuss the effect of multiple forward scattering on the energy balance of the whole terrain, we define the terrain albedo  ter as the ratio of reflected and incident radiation on the terrain.Using the fact that the reflected radiation is equal to the incident minus the absorbed radiation,  ter can be expressed as follows: where

F
in Equation 9,  ter differs from the terrain-averaged plane albedo  : With the help of the following definitions and Equation 10, Equation 9 can be rewritten to: Thus, it can be expressed to what extent the isotropic single scattering (   is ) and, moreover, the property of forward scattering (   fs ) and multiple scattering (   ms ) affect the terrain albedo.

Radiative Transfer Model
Embedded in the surface process model Alpine3D (Lehning et al., 2006), the new terrain radiation model GROUNDEYE receives interpolated real weather data with diffuse and direct broadband shortwave radiation for each pixel as well as a spatially variable plane albedo from the module SNOWPACK (Lehning et al., 2002).Details on the two-dimensional treatment of the atmosphere in Alpine3D including the isotropic approximation for solar diffuse radiation can be found in Helbig et al. (2009).
GROUNDEYE performs the three-dimensional radiative transfer by a discrete radiosity method incorporating the BRDF of snow.The basic principle of the new algorithm is the following: The visible hemisphere of each pixel p of the regular triangular grid is discretized into S solid angles s, as illustrated in Figure 2.This is done in such a way that the view factor s f is the same for all s with  1 s f S .Each solid angle s of a pixel p is identified with either "Sky" or the solid angle  s of another pixel  p during initialization: " " if "no intersection with terrain" ( , ) ( , ) else This is achieved by, as shown by Figure 2, checking the straight line through the center of the solid angle s (under zenith and azimuth angles  , p s and  , p s ) for intersections with other pixels  p .If the straight line does not intersect any other pixel, the solid angle s is identified with "Sky."For the subsequent radiative transfer, it is assumed that pixel  p covers the entire solid angle s, although normally  p covers only a part of it.In this sense,  p is said to represent the solid angle s of pixel p.On the other hand, a pixel  p may represent more than one solid angle s of pixel p, in particular if p and  p are nearby.The set    p is in other words a sample of all pixels that are seen by p.The calculation of the sky view factor sky f is then straight forward: It is the sum of all s f identified with "Sky." The governing equation of the new radiative transfer algorithm is as follows: , , , , , ) where   ( , is the reflected short wave radiance from pixel p in direction   ( , ) in units of is the single scattered radiance, it is discussed below.S is the number of solid angles into which the hemisphere is segmented, that is, the hemispheric resolution.The sum goes over all representative pixels, thus over those solid angles, which are not identified with "Sky" (see Equation 15).The angles    where  p is the plane albedo in p.This equation results from Equation 1 and the assumption that diffuse sunlight is scattered isotropically.It may be worth mentioning that Equation 16 is in fact a system of  S P coupled linear equations that satisfy the rendering equation (Kajiya, 1986), where S is the hemispheric resolution and P the number of pixels in the simulated domain.Unlike the standard radiosity approach, in which radiative transfer takes place between all mutually visible pixels (Helbig et al., 2009), the terms in the sum of Equation 16are limited to S. This allows an efficient incorporation of the BRDF.In addition, the complex sub-structuring of nearby pixels (Helbig et al., 2009) can be omitted, as the possibility of multiple representation ensures that automatically.
The quantity of interest, the global flux density global p F in point p is calculated as: ( , wherein  p I is obtained by iteratively solving the system of Equation 16.This iterative process corresponds to the multiple scattering of terrain radiation.Once the change in terrain radiation from one iteration step to the next is smaller than 1 W in all points, the iteration stops.Note that while Equation 16 determines radiative transfer only between terrain pixels, the subsequent rendering of the global radiance to a flux density in Equation 18 can be made for any points.This opens the possibility of irradiance modeling for solar panels of any position and orientation. The radiation model as described above has certain limitations.The two-dimensional treatment of the atmosphere with statistic decomposition of the measured global radiation in direct and isotropic diffuse radiation is significantly simplified (Helbig et al., 2009).However, the focus in the present study is explicitly on the terrain radiation.A more complex treatment of the atmosphere is likely to increase the anisotropy effect as the circumsolar radiation is attributed to the diffuse radiation but is scattered in reality much more like direct radiation.Furthermore, the interaction of radiation with the atmosphere between two terrain scattering events is neglected, most notably atmospheric extinction.In the investigated alpine valley, light travels an average of less than 300 meters between two terrain scattering events.This is, however, significantly shorter than the characteristic extinction length (the inverse of the extinction coefficient), which is in the range of several kilometers (Mei et al., 2017).For the modeling of very wide valleys, a distance-dependent extinction factor may be inserted into Equations 16 and 17.
It is important to note that in GROUNDEYE the following factors, which according to Warren et al. (1998) affect the BRDF of snow, are neglected: single-scattering phase function, snow grain size and surface roughness.He et al. (2017) found that the impact of snow grain size on anisotropy is very small at wavelengths <  1 m.The assumption of spherical snow grains, on the other hand, tends to overestimate the forward scattering (Dang et al., 2016).The same is true with respect to snow surface roughness.Hudson et al. (2006) carried out BRDF measurements on the East Antarctic Plateau and found, most likely due to the presence of sastrugi on the snow surface, a significantly reduced forward scattering peak compared to DISORT model results.However, due to generally smaller macroscale surface roughness and higher snow temperatures, these results cannot be directly transferred to alpine terrain.The uncertainties are therefore difficult to quantify, but are most likely much smaller than the domain dependence of the effects investigated.To avoid modeling unrealistically flat surfaces, we fixed the BRDF  for    , and likewise for  v .

Broadband BRDF
The broadband (shortwave) BRDF of snow used for this study is based on the numerical model by Mishchenko et al. (1999), which has been shown to reproduce the important properties of measured snow BRDF (Dumont et al., 2010), and on the module SNOWPACK that models the plane albedo using an empirical model with meteorological parameters and parameters of the snow cover (Lehning et al., 2002).The Mishchenko model generated the monochromatic BRDF's assuming a power law distribution of sphere radii with an effective radius of  50 m and embedding the wavelength-dependent complex refractive index from the study by Warren and Brandt (2008).The monochromatic BRDF's       ( , , , , ) were further weighted according to the solar spectrum whereby k a is based on the normal terrestrial spectral irradiance AM1.5 days G-173-03 provided by the American Society for Testing and Materials (ASTM, 2012).For the present study, 111 monochromatic BRDF's spaced 20 nm from 300 nm to 2500 nm were used.In a second step,  M was scaled to match the SNOWPACK albedo  S : , , ).
The use of a broadband BRDF is an approximation as the central Equation 16of the radiative transfer model is correct only for a monospectral BRDF.In other words, the weighting in Equation 19 is correct only for single scattering because the spectrum shifts with each reflection towards wavelengths with comparatively high albedo.Since the forward scattering effect (Equation 6) is defined for single scattering, it is not affected by this spectral shift.However, the multiple scattering effect (Equation 7) is likely to be underestimated by up to 10% because, due to the spectral shift, the albedo increases with each scattering event.8) is overestimated by the broadband treatment because the anisotropy is less pronounced at wavelengths with high albedo.However,  FSE MSE is already of a higher order than MSE, the spectral shift can therefore be neglected.Most significantly, the calculations of absorbed radiation are affected by the use of a broadband BRDF.Because the albedo  of snow is usually much larger than 0.5, in relative terms the absorption   1 changes much more than .For example, if the spectral shift increases the 2-fold scattering albedo by 10 percent from 0.8 to 0.88, the absorption decreases from 0.2 to 0.12, that is, by 40 percent.Apart from underestimated multiple scattering, the broadband method has no effect on the calculations for solar panels, since no absorption is calculated for them.

Modeling Solar Panels
The implementation of solar panels in the model is straightforward: just as Equation 18 is evaluated for terrain pixels, it can be done for rectangular panels of any size and orientation above the surface.The irradiance is calculated only for a point in the center of the panel, while its size determines the shadow cast, which is not covered by the preceding radiative transfer.Well-placed solar modules are indeed very effective at casting shadows.This can have an effect on productivity, especially with bifacial panels, as the back side usually has a shadow in the field of view.A function was therefore added to GROUNDEYE to calculate the shadows of the solar panels.It does not model direct cast of shadows on other panels but the cast of shadows on the terrain.The model considers exclusively shadows that are caused by the sun's direct radiation.Although diffuse radiation is partly shaded by the panels as well, this diffuse shadow was not modeled as it is more complex and less important in most circumstances.The function works as follows: If a certain viewing direction of a panel meets terrain, the function calculates whether the line between this point and the sun intersects the panel itself or another panel.In case of intersection the viewing direction is considered to see shade.Therefore, no direct radiation is scattered by this location.
GROUNDEYE is able to calculate the optimal panel tilt and alignment angle of solar panels.For this purpose the Nelder-Mead blackbox optimization method (Nelder & Mead, 1965) was implemented.

Model Validation
The new model GROUNDEYE has been verified with analytical solutions based on a topography shown in Figure 3b: This so-called hemispherical cavity has the property that at any point on the surface the skyview factor is  1 2 sky f (Treuenfels, 1963;Manners et al., 2012).Thus in case of non-absorbing Lambertian surfaces, half of the radiation is scattered back into the cavity at every scattering event, so the total terrain RÜTTE ET AL.

10.1029/2020JD034333
radiation can be calculated easily with a geometric series.In Figure 3a, model results are compared with these analytical considerations for three different incident types of radiation: Isotropic sky-diffuse radiation, direct radiation at  75 zenith angle and zenith radiation on random pixels.For all configurations, the deviation is well below 1%.
At the Totalp in the Swiss ski resort Davos-Parsenn, as part of a photovoltaic test facility, global (and direct) SW radiation is measured under various angles of incidence.Given the fact that for strongly inclined sensors a substantial portion of measured radiation is terrain radiation, these measurements are very well suited to validate GROUNDEYE.A digital elevation model with a grid size of 10 m was used such that the domain of the simulation contains the horizon seen by the test facility.Furthermore, horizontal global radiation and direct radiation from the test facility as well as snow height, precipitation, temperature and humidity of the nearby Weissfluhjoch are used.For validation, we chose measurements of two south-oriented radiation sensors with an inclination of 30 and 90 degrees respectively and one north-oriented radiation sensor with an inclination of 90 degrees.Figure 4 shows for these sensors a comparison of measurements and simulation both for four clear-sky days from November to February as well as for three whole weeks in March of the winter 2017-2018, in the following called all-sky days.Under clear-sky conditions, the average deviation is less than 6% in each case.For the all-sky days the radiation is underestimated with the exception of the north-facing  90 sensor.This is most likely a consequence of the isotropic sky model, which neglects circumsolar radiation (see Reindl et al. (1990b)).RÜTTE ET AL.The similar performance of the south-facing  30 and  90 sensors, and the fact that terrain radiation is much more important for the latter, indicates that the errors of terrain radiation are small compared to the errors in sky radiation.Furthermore, the agreement is improved by including anisotropic, multiple scattering and self-shading, which justifies the use of the new sophisticated algorithm.Sole exception is the north-facing panel, for which the deviation increases through the inclusion of self-shading and multiple forward scattering.However, the deviations are below 6%.Since direct radiation measurements are only available for the time periods used for validation, a statistical decomposition of global radiation by (Reindl et al., 1990a) has been used for the full year studies.Comparison with measured values at the Totalp shows that the direct solar radiation tends to be overestimated by the statistical decomposition under all-sky conditions.

Alpine Irradiance Simulation
The Alpine valley "Meierhofertäli" in the Swiss ski resort Davos-Parsenn (see Figure 5a) is a typical inner-Alpine valley with a length of 2 km, a width of about 800 m and a depth of roughly 200 m.Situated at an average elevation of 2331 . .m a s l., the valley descends 600 m towards the east (see Figure 5b).As can be seen in Figure 5c, the terrain-view factor in the Meierhofertäli is on average occur on the north-facing slope and in the canyon-like valley outlet in the east.Using meteorological data from a station within the domain and a digital elevation model with 10 m resolution, radiative transfer has been modeled for the period September 2017 -August 2018 at hourly resolution.Irradiance was calculated both for the terrain as well as for different solar panels.Radiative transfer mostly converged to the criterion described in Section 3.1 within 3-6 iterations.The entire simulation took about one day for a hemispherical resolution of S = 1200 on a computer with 20 CPU.
The effect of snow-covered terrain on the irradiance of solar panels is studied for the test site on the Totalp ( 2475. .m a s l., marker pin in Figure 5a).The model allows to break down the simulated total irradiance into direct solar radiation, sky-diffuse radiation and terrain radiation.The latter is further subdivided into radiation from a hypothetically snow-free terrain with a soil albedo of   0.2 soil , the additional radiation due to the much higher albedo of snow (albedo effect), multiple scattering effect (MSE), forward scattering effect (FSE), and finally the loss due to panel-induced shading of the surrounding terrain.

Effect of Multiple Forward Scattering on Terrain Irradiance
Figures 6 and 7 show the spatial distribution of different simulated radiative quantities within the studied valley (see section 3.5), on February 19, 2018 respectively on May 6, 2018 at 12:00 p.m. under clear-sky conditions.By using the two different dates in February and May, the terrain radiation can be studied on a clear winter day as well as during snow melt.The total irradiance in Figures 6 and 7a is the sum of direct, sky-diffuse and terrain radiation in the pixel plane.Clearly visible on February 19 is the shadow cast on the north-facing slope as well as patterns due to different slope inclinations.The large spatial variability from Figures 6 and 7b show the incident terrain radiation as a percentage of the total irradiance.This contribution is particularly substantial for shady slopes on February 19: widely more than 40%, locally well over 60%.In absolute terms, the terrain radiation is on average 2 46 Wm and locally reaches up to 2 185 Wm .Such remarkable values result from the coincidence of high albedo and high complexity of the snow-covered terrain.Over the entire domain, terrain radiation accounts for 8.3% of the total irradiance.On May 6, this value is 5.4% ( 2 49 Wm ).Since the north-facing slope receives direct solar radiation on May 6, terrain radiation only locally accounts for more than 20 % of the total radiation.Averaged over the period September 2017-August 2018, terrain radiation accounts for 5.8% of the total irradiance or To what extent is the terrain radiation affected by the forward scattering of light by snow?Figures 6 and 7d show the percentage difference between single scattered terrain radiation based on anisotropic scattering and isotropic scattering respectively, that is, the forward scattering effect (FSE) defined in Equation 6.It RÜTTE ET AL.    6; (e) multiple scattering effect in percent as defined in Equation 7; (f) Amplification of multiple scattering through forward scattering in percent as defined in Equation 8.  6; (e) multiple scattering effect in percent as defined in Equation 7; (f) Amplification of multiple scattering through forward scattering in percent as defined in Equation 8. can be seen that this effect is very pronounced on February 19: More than two thirds of the pixels show a difference of more than 10%, for some pixels the effect exceeds 2 40 Wm .It has been found that the direct irradiance correlates moderately with the FSE (  0.61 r ).This indicates that sun-facing slopes benefit from the effect, while slopes directed away from the sun get less than in the isotropic case.A probable explanation can be found in the distinct forward scattering peak of snow for small incident elevation angles (which is often the case in winter) and the fact that slopes facing the sun tend to get forward scattered light.Further, the forward scattering is peaked at low elevation angles of outgoing light.Since terrain is mostly seen under small elevation angles, the forward scattering not only leads to a spatial redistribution but also to a positive average effect of (4.7%) in terrain radiation ( 2 1.8 Wm ).The FSE is somewhat less pronounced on May 6, when more than half of the pixels show a difference of more than 10%.In contrast to February 19, the FSE on the south-facing slope is widely negative, correlation with direct irradiance is weaker (  0.41 r ).The reason for this is probably the darkening at grazing angles of the BRDF for the pixels of the flat valley floor (Dumont et al., 2010), which are hit by the high May sun at a large angle.This leads to a terrain averaged FSE of 1.9% ( 2 0.8 Wm ).Over the whole day, the FSE is positive with 5.2%, probably due to the lower elevation of the sun outside noontime.On an annual average, this effect is 4.3% ( 2 0.3 Wm ).
Figures 6 and 7e show the percentage increase of terrain radiation due to the process of multiple scattering, that is, the multiple scattering effect (MSE) defined in Equation 7. As can be seen in Figure 6e, the MSE ranges on February 19 widely between 8% and 20%, locally it's over 40%.The absolute values are largely between 2 2 Wm and 2 6 Wm , locally over 2 20 Wm .On May 6, the MSE ranges widely between 5% and 12%, the terrain average of 9.1% is lower than in February (13.1%) due to the lower albedo.On an annual average, the MSE over the entire terrain is 11% ( 2 0.9 Wm ), locally up to 2 3.9 Wm .The MSE correlates only weakly with the terrain-view factor ter f (  0.24 r , February 19), but more strongly with the average terrain-view of the visible terrain (  0.55 r ).This quantity, we called it friends' terrain-view factor  ter f , is shown for the Meierhofertäli in Figure 5d.It is pronounced in narrow canyons or cavities, but much less on prominent slopes.On the flat valley floor, on the other hand,  ter f is usually higher than ter f .Averaged over the entire terrain,  ter f is 18 percent higher than ter f .This apparent asymmetry is generally known as friendship paradox (Eom & Jo, 2014).
To keep things simple, the forward scattering effect FSE refers to single scattered terrain radiation.However, the forward scattering also takes place in multiple scattering.To estimate this coupling, Equation 8 defines  FSE MSE, that is, the difference between anisotropic and isotropic multiple scattering effect.Figures 6f  and 7a show that the anisotropic multiple scattering effect is widely enhanced by 5%-25% when compared to the isotropic multiple scattering effect, locally by over 40%.On February 19 and May 6, the effect is negative in 12% and 4% of the pixels, respectively, the terrain average is 18% and 21% percent, respectively.For the whole domain, the effect is on an annual average 19%.Thus, forward scattering significantly amplifies multiple scattering.This is probably because incident terrain radiation usually hits the surface at glancing incidence angles whereby the forward scattering peak is very pronounced.Nevertheless, it must be noted that the consideration of anisotropy in multiple scattering has clearly a higher order effect compared to the anisotropy of single scattering or isotropic multiple scattering.

Comparison With Simple Terrain Radiation Algorithms
Many radiation models use highly simplified local approximations of the terrain radiation based on the so-called homogeneous view factor approach (Ruiz-Arias et al., 2009;Šúri & Hofierka, 2004;K. Wang et al., 2005).The terrain radiation in point p is thereby estimated as horizontal direct and diffuse radiation Three different variations of this highly efficient method are compared in the following to the present complex model for the the Meierhofertäli (see Figure 5a) on February 19, 2018 at 12:00 p.m. under clear-sky conditions.The first version by Mészároš and Miklánek (2006) or by Šúri and Hofierka (2004) uses an approximation of the terrain-view factor ter f based only on the local inclination The performance of this method is shown in Figure 8a.Deviations are significant and systematic.Especially in the floor of the valley, terrain radiation is massively underestimated.50% of the pixels deviate by more than 50%, the correlation of the approximation with the complex model is moderate (  0.61 r ).The reason for this is primarily the calculation of the terrain-view factor, which is based exclusively on the local slope inclination and is therefore considerably underestimated in the floor of the valley.Because the terrain radiation in the floor of the valley is small compared to the slopes, the underestimation in the terrain average with 28% is less than what might be expected from Figure 8a.In the daily mean, terrain radiation is underestimated by 26%, in the annual mean by 18%.
The second variation of the simple method of Equation 22 uses correct terrain-view factors, calculated on the basis of the entire digital elevation model.The spatial patterns of its deviation from the complex model are shown in Figure 8b.Compared to Figure 8a, the deviations are smaller and less systematic.They result from neglecting forward and multiple scattering, and from the use of a local albedo and homogeneous terrain irradiance, the latter being probably the most important.The largest deviations occur in a band near the north-facing slope, as well as on the sun facing side of the valley outlet in the east.These are locations that see shady terrain with significantly less irradiance than the horizontal global radiation used in Equation 22.In the domain average, terrain radiation is overestimated by 17.6%.The correlation with the complex model is strong (  0.81 r ), as ter f largely governs terrain radiation.
Figure 8c shows the comparison of this method for the annual irradiance from September 2017 to August 2018.The correlation with the complex model is very strong (  0.95 r ).The deviation is strongest at RÜTTE ET AL.

10.1029/2020JD034333
14 of 20 locations that see a lot of shaded terrain, most pronounced in the valley outlet in the east.These locations have a high friends' terrain-view factor as can be seen in comparison with Figure 5c.In fact, the deviation correlates much more strongly with the friends' terrain-view factor (  0.68 r ) than with the terrain-view factor (  0.27 r ).Not only does this method overestimate the annual terrain radiation at these specific points, the terrain average too is overestimated by 21.3%.The reason for this is probably the following: The simple models of Equation 22 use the horizontal global radiation as a basis for the isotropic terrain radiation.Due to the non-planar topography, however, this radiation distributed over more area, and the average global radiation is therefore smaller than the horizontal global radiation.The area of Meierhofertäli is 15% larger than the horizontal projection, so the horizontal global radiation should be reduced by this value.This would leave an annual average deviation of 5.5%.
To take into account multiple scattering, Sirguey (2009) proposes to extend ter f by the sum of a geometric sequence with common ratio  f ter : where  is the average albedo and f ter is the average terrain-view factor of the valley.This third variation of the simple method increases the overestimation of the annual terrain radiation in the domain average to 29.3%, while correlation with the complex model is unaffected since the innovation of Equation 24 is a constant factor.For the same reason, spatial patterns are the same as those of the second variation in Figure 8b and are therefore not shown separately.
It can be noted that in the spatial mean the simple local approximations agree well with the elaborate model under the condition of correct calculation of the terrain-view factors.Locally, the deviations are are often substantial however with more than 2 80 Wm .However, the inclusion of the correction factor for multiple scattering is not an improvement, because using the horizontal global radiation as a basis for the isotropic terrain radiation already leads to a systematic positive deviation larger than the positive correction factor.

Effect of Multiple Forward Scattering on Terrain Energy Balance
In the previous sections, incident terrain radiation has been studied in detail.However, it has not been studied what part of the incident terrain radiation is absorbed by the snow cover, or to what extent the terrain albedo is changed.In Section 4, it was shown that for the simulated valley, both forward scattering and multiple scattering increase the average terrain radiation.It is therefore to be expected that the radiation absorbed by the terrain will increase and thus the terrain albedo  ter , as defined in Equation 9, will be reduced by multiple forward scattering.
Figure 9a shows the terrain averaged plane albedo  as defined in Equation 10.Clearly visible are the first snowfalls in autumn as well as the gradual snow melt from April on.From December to March the albedo is mostly between   0.8 and   0.9.Figure 9b shows the components of the terrain albedo according to Equation 14.It can be seen that isotropic single scattering reduces the terrain albedo (   is ) by about 1-2 percentage points.The effect is most pronounced after the first snowfalls in autumn and especially during the snow melt from April to July, when it amounts to    2.3%Using Figure 9a, it can be seen that   is is largest for terrain averaged plane albedos between   30% and   70%.This is due to the following reason: As shown in Equation 11, as well as to the terrain radiation, which itself is proportional to the plane albedo .The product ) , is largest for   50%., respectively.Forward scattering reduces the terrain albedo most when the forward scattered radiation hits snow-free terrain and the backscattered radiation hits snow.This is the case in autumn and spring, when south-facing slopes are snow free; in winter, the effects of forward scattering on the terrain albedo and thus the energy balance of the entire terrain is marginal.The reduction of the terrain albedo by multiple scattering   ms shows a similar seasonal pattern as that of the isotropic single scattering   is , but decreases to lower values in summer than in winter.The reason for this is that the former is proportional to ) and is therefore more strongly suppressed by low values of  in summer.
As mentioned in Section 3.2, the use of a broadband BRDF is unsuitable for precise calculation of absorbed terrain radiation.However, the main point remains valid: the terrain albedo is only marginally reduced by forward scattering and multiple scattering, because the terrain radiation is increased by these effects in times when the albedo is high.

Effect of Multiple Forward Scattering on Irradiance of Solar Panels
Solar panels are very efficient radiation sinks, it thus appears worthwhile to install them in highly reflective terrain.The effect of snow-covered terrain on the irradiance of solar panels is studied for the test site on the Totalp ( 2475. .m a s l., marker pin in Figure 5a).Figure 10a shows the different radiative components for the south-facing front side (  90 S ) and the north-facing back side (  90 N ) of a vertical bifacial panel, both for a clear-sky day at the end of February and for a overcast day at the beginning of March 2018.At both times, the valley was completely covered by snow.
On the clear-sky day,  90 S receives two thirds direct radiation, about 30% comes from the terrain.With a share of 6.8% of the total radiation, the FSE is more than twice the size of the MSE.On the back side  90 N , together the negative FSE and the shadow cast by the panel reduce the total irradiance by more than 25%.The exposure on  S on the overcast day is one-quarter that of the clear-sky day.On the back side however, the exposure remains nearly the same, as the lower terrain radiation is compensated by higher sky-diffuse radiation and the absence of both direct shadow cast and negative FSE.
Figure 10b shows average irradiance values for the entire period September 2017-August 2018 as well as for the months November-January which normally have the lowest production of solar electricity in northern latitudes.In addition to the two sides of the vertical panel, corresponding results are shown for a south-facing panel with an inclination of 30 degrees (   14) by isotropic single scattering   is (Equation 11), forward scattering   fs (Equation 12), and multiple scattering   ms (Equation 13). 90 S is 47% higher than on  30 S . The steep inclination combined with the high snow albedo leads to a significant shift of exposure from summer to winter as stated by Kahl et al. (2019)   it can be noted that the exposure from November to January can be increased by 55% without lowering the annual average with regard to terrain-neglecting calculations ( 1Opt ).In this context, the forward scattering plays an important role, accounting for 27% of the terrain radiation and 8.1% of the total radiation for  90 S in winter.On an annual average, FSE and MSE for south-facing panels are less important with about 0.8% of total irradiance each for  30 S and 3% for  90 S .

Conclusions
With the help of the new model GROUNDEYE, the radiative transfer of sunlight in mountainous snow-covered terrain was studied in great detail.Shadow casting, variable surface albedo, forward scattering and multiple scattering lead to a very complex spatial radiation distribution with local maxima.For a typical Alpine valley near Davos, Switzerland, representing high-mountain valleys in mid latitudes, we showed that on winter days, the share of ground-reflected radiation in the radiative balance of the terrain is locally above 60%, with absolute values exceeding 2 180 Wm .The forward-scattering property of snow leads to a significant redistribution of radiation within the terrain in favor of sun-facing slopes in the order of a few Watts per square meter (locally up to 2 50 Wm ) on clear-sky days in winter.However, the forward scattering does not lead to a zero-sum redistribution; it increases the terrain radiation on average and slightly lowers the terrain-average albedo.In particular, we were able to show that the multiple scattering of radiation is amplified by about 20% through forward scattering in the studied domain.As most of the additional radiation is ultimately reflected from the terrain due to the high albedo of snow, the effect of forward and multiple scattering on the area-wide albedo is small.We showed that in the studied Alpine valley multiple forward scattering decreases the terrain albedo by less than 0.2 percentage points on annual average, while isotropic single scattering reduces the terrain albedo by 1.45 percentage points.
Simple methods for the calculation of terrain radiation based on the so-called homogeneous view factor approach were evaluated.The comparison with the complex model shows that the large-scale averages are well approximated if the terrain-view factors are calculated accurately.Since basing the isotropic terrain radiation on the horizontal global radiation already leads to a systematic positive deviation that is, larger than the effect of multiple scattering, we do not recommend the use of an extra factor accounting for multiple scattering.An estimation of the terrain-view factor exclusively with the local inclination is strongly discouraged, because it leads to a strong underestimation of the terrain radiation especially for the valley floor.At small scales, the accuracy of the simple methods is limited, even if the terrain-view factors are calculated accurately.
The complex radiative transfer in snow-covered terrain with its high surface albedo has a direct impact on the productivity of photovoltaic installations in Alpine regions.We showed that for a south oriented, vertical panel at the study site, the terrain radiation accounts for 30% of the total irradiance from November to January.7.6% alone are owed to the forward scattering.In sunny locations, where the optimal panel orientation is controlled by direct sunlight, the careful modeling of the terrain radiation is important for an efficient increase of winter radiation in the following way: Reflected radiation from the ground outweighs diffuse radiation from the sky, leading to high optimal panel tilts.This, in combination with the sun's low elevation in the Alpine winter, shifts solar production potential from summer towards winter.We were able to show that the winter yield can be increased by more than 50% without reducing the annual yield with respect to terrain-neglecting calculations.This increase in winter production is of crucial importance for a future carbon-free energy supply in the mid-latitudes.

Figure 1 .
Figure 1.Visualization of radiative transfer in mountainous, snow-covered terrain: The large arrows illustrate the trajectory of a light beam, the small arrows symbolize anisotropic scattering.(Ernst Ludwig Kirchner, Davos in Winter.Davos in Snow, in: Kunstmuseum Basel, Sammlung Online; Modified).

,
defined in Figure2.Further,  p is the BRDF defined in Equation 1 in p.While the sum on the right hand side of Equation 16 controls the multiple scattering in the terrain, single p I corresponds to the single scattered radiance.The latter is reflected radiance resulting from direct and diffuse solar irradiance

Figure 2 .
Figure 2. Sketch of the geometric principle of GROUNDEYE: The hemisphere of the triangular pixel p is discretized.The straight line through each solid angle s of it is checked for intersection with other triangular pixels  p .In case of intersection, the solid angle s of pixel p is identified with the solid angle  s of pixel  p (right).If the line does not intersect any other pixel, it points to the sky (blue solid angle).

Figure 3 .
Figure 3. (a) Comparison of theoretical and simulated average terrain radiation within the hemispherical cavity as a function of the iteration steps for different incident configurations.The final convergence accuracy is listed in the lower right corner.(b) The hemispherical cavity.

Figure 4 .
Figure 4. Average deviation between model and measurement of SW flux density for all-sky and clear-sky days at different inclinations at the test facility Totalp.Hourly standard deviation between measurement and the full model is indicated by the gray ribbon.Anisotropic and multiple scattering as well as selfshading are omitted for the reduced model (dashed line).The average daily deviation  and its standard deviation  is shown for each case.

270
Wm to more than 2 1100 Wm is a feature of the mountainous topography.In comparison, the spatial variability is less pronounced on May 6, while the domain average is with

28. 8
Figures 6 and 7c show the modeled plane albedo.While it varies only minimally around   0.83 on Feb- ruary 19, the differences are larger on May 6: Parts of the valley exit are already free of snow with   0.2.Large parts of the remaining valley have values of   0.6, in some places up to   0.8.The annual average albedo is   0.67.

Figure 6 .
Figure 6.Results of the Meierhofertäli Simulation for February 19, 2018 at 12 p.m., at a solar elevation angle of 31.4 degrees: (a) Total irradiance in 2 Wm ; (b) Terrain radiation as percentage of total irradiance; (c) Plane albedo; (d) Forward scattering effect in percent as defined in Equation6; (e) multiple scattering effect in percent as defined in Equation7; (f) Amplification of multiple scattering through forward scattering in percent as defined in Equation8.

Figure 7 .
Figure 7. Results of the Meierhofertäli Simulation for May 6, 2018 at 12 p.m., at a solar elevation angle of 56.5 degrees: (a) Total irradiance in 2 Wm ; (b) Terrain radiation as percentage of total irradiance; (c) Plane albedo; (d) Forward scattering effect in percent as defined in Equation6; (e) multiple scattering effect in percent as defined in Equation7; (f) Amplification of multiple scattering through forward scattering in percent as defined in Equation8.
local or semi-local albedo  p times the terrain-view factor,

Figure 8 .
Figure 8. Evaluation of simple approximations of terrain radiation based on the homogeneous view factor approach for the Meierhofertäli.Shown is the relative deviation of the approximation from the complex model.(a) Method with roughly estimated terrain-view factor (Equation 23) for February 19, 2018 at 12:00 p.m.(b) Method with exact terrainview factor for February 19, 2018 at 12:00 p.m. (c) Method with exact terrain-view factor, annual mean.

Figure 9 .
Figure 9. (a) Terrain averaged plane albedo  as defined in Equation10for the totalp valley (see section 3.5) from september 2017 to august 2018.(b) Change in terrain albedo  ter (Equation14) by isotropic single scattering 
. Wm which corresponds to 40% of the radiation received by  90 S. From November to January it is only 28%.The optimal inclination with respect to the annual exposure at the test site was calculated both comprehensively ( Opt are likewise shown in Figure10b.It turns out: By considering the terrain radiation, the annual average increases by 12%, the winter radiation by 37%.Since 279.3