## Abstract

The influence of topography on the radiation balance in complex terrain has so far been investigated either with very simple or very sophisticated approaches that are limited, respectively, by an uncontrolled spatial representation of radiative fluxes or heavy computational efforts. To bridge this gap in complexity, this paper proposes the radiosity approach, well known in computer graphics, to study anisotropic reflections of radiation in complex terrain. To this end the radiosity equation is rederived in the context of three-dimensional radiative transfer. The discretized equation is solved by means of an adapted version of progressive refinement iteration. To systematically study terrain effects, the geometrical disorder provided by the topography is considered in its simplest approximation by Gaussian random fields. These model topographies capture the most important length scales of complex terrain, namely a typical elevation and a typical valley width via the variance and the correlation length of the field, respectively. The mean reflected radiation is computed as a function of these length scales and sun elevation, thereby explicitly addressing finite system sizes and grid resolutions. A comparison with an isotropic parameterization of terrain reflections reveals that mean values are similar whereas spatial distributions vary remarkably. It is also shown that the mean reflected radiation in real topography is reasonably well characterized by the Gaussian approximation. As a final application of the method, the effective albedo of a topography is shown to vary with sun elevation and domain-averaged albedo, leading to albedo differences up to 0.025.

## 1. Introduction

Downwelling radiation is a crucial component of the energy balance required for meteorological, climate, and land surface models. Even in the absence of topography it is a complicated issue to compute the downwelling shortwave flux with realistic radiative transfer models (RTMs) or with empirical models by means of parameterizations based on, respectively, atmospheric conditions or the clearness index and a point measurement of global radiation [cf. the parameterizations of, e.g., Liu and Jordan (1960), Orgill and Hollands (1977), Erbs et al. (1982), Iqbal (1983), and Reindl et al. (1990)].

In the presence of topography two additional effects have a remarkable influence: (i) shading and (ii) multiple terrain reflections or emissions. Particularly in highly complex, snow-covered alpine terrain, the latter effect plays an important role, and different approaches with different levels of complexity exist to take them into account.

There are basically two state-of-the-art model approaches to compute three-dimensional radiative transfer (RT) in complex terrain: deterministic and stochastic. The spherical harmonics discrete ordinate method (SHDOM; Evans 1998) is an example of the deterministic approach. It is widely applied because it is flexible and publicly available (Cahalan et al. 2005). SHDOM is adopted (e.g., in Degünther et al. 1998) to conduct a case study on the interaction between a cloud free atmosphere and a model domain with inhomogeneous surface albedo for ultraviolet (UV) irradiance. Another deterministic model approach is the discrete anisotropic radiative transfer (DART) model that simulates radiative transfer in heterogeneous three-dimensional landscapes with the exact kernel and discrete ordinate methods (Gastellu-Etchegorry et al. 1996, 2004). The second model approach is the stochastic Monte Carlo (MC) method (Marchuk et al. 1980). MC is also widely applied and flexible for three-dimensional RT calculations. It computes the radiation by tracing the stochastic “trajectory” of a photon through the atmosphere. All interactions of a photon with the atmosphere and the ground as boundary surfaces can be incorporated. Two MC models exemplify this approach. The first example uses a forward photon tracing MC method that calculates three-dimensional radiation fluxes while accounting for inhomogeneous surface albedo values and topography in cloudy atmospheres (Mayer 2000; Kylling et al. 2000a). The second MC model example computes the shortwave (SW) radiant flux for a real three-dimensional topography only underneath a clear-sky atmosphere and for homogeneous surface albedo values (Chen et al. 2006). An intercomparison of complex three-dimensional RTMs is given by Cahalan et al. (2005).

On the other hand, various simple models exist for different purposes to estimate terrain reflected or emitted radiation as an isotropic contribution in terms of the so-called sky view factor, *F*_{sky}. Such an approach can be found in energy or radiation balance models (e.g., Nunez 1980; Moore et al. 1993; Corripio 2002; Oliphant et al. 2003; Hock and Holmgren 2005; Wang et al. 2006), diagnostic meteorological models (e.g., Bellasio et al. 2005), or the nonhydrostatic mesoscale model (NMM) with a grid- and subgrid-scale radiation parameterization scheme (Müller and Scherer 2005). An extension of these approaches that takes anisotropic terrain reflection and emissions into account has been put forward by Olyphant (1986) for longwave (LW) radiation and by Letsinger and Olyphant (2007) for SW and LW radiation.

Presently, the gap in complexity between realistic radiative transfer models and models using the so-called isotropic view factor approach is large. However, surface process models such as Alpine3D with applications in hydrology (Lehning et al. 2006) or models handling preferential deposition of snow (Lehning et al. 2008) usually have their main scope well within this gap. To partly bridge it, the radiosity approach is presented in this paper to compute the three-dimensional shortwave radiation balance in complex terrain. The radiosity approach was first introduced in computer graphics by Goral et al. (1984) from techniques originally developed in thermal engineering to describe the exchanges of radiant energy between surfaces. A detailed historical overview of computation methods for radiative energy exchange is given by Siegel et al. (1991). The radiosity approach can be regarded as a three-dimensional RTM where the topography is treated as a Lambertian emitter/reflector surface and the sky as an emitting surface of zero reflectivity. Borel et al. (1998) considered the radiosity approach as a possible methodology for computing the radiation balance in complex terrain, although a systematic investigation of this methodology has never been performed. Although the radiosity method has been applied in urban meteorology for a long time (Verseghy and Munro 1989; Krayenhoff and Voogt 2007), its intrinsic relation to RTM has never been derived explicitly.

A systematic study of the influence of topography on reflected radiation requires a suitable ensemble of topographies. Because even the simplest case of Gaussian random fields has never been investigated, we are using these as our model topographies. With these ensembles complex terrain is characterized by two fundamental length scales: a typical valley-to-peak elevation difference *σ* (height of mountains) and a typical lateral extension *ξ* (width of mountains).

The paper is organized as follows: In section 2 we rederive the radiosity equation in complex terrain from general RT considerations. The numerical solution of the radiosity equation is presented in section 3. In section 4, we specify our ensemble of model topographies. In section 5, simulations of SW RT in the model topographies using the radiosity model are presented, revealing the sensitivity of model results to geometric parameters of the topography. We compare our results to those obtained with the isotropic view factor approach and we compute an important reflected radiation parameter: the effective terrain albedo. The results are discussed in section 6.

## 2. Radiosity approach in complex terrain

### a. Radiative transfer in complex terrain

To enable the use of the radiosity equation for the surface radiation balance over complex topography, we start from basic radiant quantities to keep track of the assumptions on which the model relies. The broadband radiative flux *P*(*x*) (integrated over a fixed but arbitrary wavelength interval) received at a point *x* measured in W m^{−2} is given by the hemispherical integral

Here *ϑ*(*x*) is the zenith angle in the sloped coordinate system at point *x*, which has its *z* axis in the direction of the surface normal (cf. Fig. 1) and *R*[*x*, Ω(*x*)] is the radiance or specific intensity in W m^{−2} sr^{−1} received at *x* within the solid angle Ω(*x*). Now we decompose the integral over the hemisphere in (1) into directions pointing to the sky and to the terrain—*P*(*x*) = *P*_{sky}(*x*) + *P _{t}*(

*x*)—and assume that

*P*

_{sky}(

*x*) can be considered an independent radiation source for the system. This assumption excludes any feedback on downward radiances

*R*[

*x*, Ω(

*x*)] from the sky, which implies that atmospheric backscattering of radiation to the surface is not considered for clear-sky conditions. We assume this to be a reasonable, although simplified, starting point.

Now we rewrite the integration in (1) over directions pointing to the terrain as an integration over nonoccluded parts of the whole terrain surface *A* by expressing the solid angle at *x* in terms of the surface element at *x*′ via *d*Ω(*x*) = *dA*(*x*′) cos*ϑ*(*x*′)/|*x* − *x*′|^{2}. Then (1) can be written as

Here *R*(*x*, *x*′) is the radiance received at *x* from the differential solid angle around the direction pointing to *x*′ and *χ*(*x*, *x*′) denotes mutual visibility; that is, *χ*(*x*, *x*′) = 1 if *x* “sees” *x*′ and is 0 otherwise. The received radiance *R*(*x*, *x*′) is determined by the usual one-dimensional radiative transfer equation (RTE) within the atmospheric column between *x*′ and *x*. In the absence of scattering within this column, the received radiance is given by Liou (2002) as

in terms of the transmittance function *τ* and assumed isotropic atmospheric volume sources *Q*(*x*″). The boundary condition—that is, the radiance emitted at *x*′ into a differential solid angle around the direction pointing to *x*—is denoted by *R*_{0}(*x*′, *x*). Furthermore, we assume that all surface sources are isotropically emitting; that is,

### b. Radiosity equation for shortwave radiation

In the following, we concentrate on the incident broadband global radiation flux *P*(*x*) = *S _{g}*(

*x*). A detailed discussion for the LW flux is left for future work. Because our focus is a systematic study of the influence of topography on the radiation balance, it is reasonable to start from the simplest version of the model obtained from the following additional simplifications:

(i) We assume, in conformance with other simple models, that the atmosphere contribution to radiance can be neglected in (3) for the radiance leaving

*x*′ and incident on*x*(i.e.,*Q*= 0 and*τ*= 1), which implies*R*(*x*,*x*′) =*R*_{0}(*x*′,*x*). Thus, the radiance remains invariant along the path between*x*and*x*′ and the only contribution comes from the surface itself.(ii) We assume ideally diffuse emitting surface sources (i.e., Lambertian surfaces) in terms of the SW albedo

*α*. Hence, the radiance*R*_{0}(*x*′) is related to the incident SW flux*S*(_{g}*x*′) at*x*′ via*R*_{0}(*x*′) =*α*(*x*′)*S*(_{g}*x*′)/*π*. Lambertian surfaces are a common assumption in RTMs.

Here, *S _{b}* denotes the direct solar beam flux and

*S*the diffuse flux in W m

_{d}^{−2}from a theoretical unobstructed hemisphere. The angle between the surface normal at

*x*and the sun is denoted by

*θ*(

*x*) and the indicator function

*χ*

_{sun}(

*x*) denotes illumination; that is,

*χ*

_{sun}(

*x*) = 1 if

*x*is in the sun and

*χ*

_{sun}(

*x*) = 0 if

*x*is in shadow.

In writing down (5) we have implicitly defined the sky view factor

in the integration of the isotropic diffuse component. By defining *ϑ _{h}*[

*ϕ*(

*x*)] as the horizon angle in the sloped coordinate system, it follows that

or, equivalently,

Note that (8) reveals the right formula for the sky view factor in terms of the horizon angle *ϑ _{h}*, which has to be measured in the

*sloped*coordinate system between the surface and the horizon direction. We stress that the sky view factor (6) is the ratio of the flux received by an inclined surface from the visible part of the sky to the flux from a theoretical unobstructed hemisphere. Note that it is not simply the solid angle of the visible sky and it cannot be trivially evaluated in a horizontal coordinate system. The latter two aspects are often ignored or not explicitly stated in the literature.

Indeed, the assumption of an isotropic sky radiance is questionable because the circumsolar radiation component caused by strong forward (Mie) scattering might give rise to a pronounced anisotropy. However, here we follow simple model approaches that assume isotropy for the sky.

By applying the above simplifications [(i)–(iii)] to the radiation balance [(2)–(4)], we end up with a linear integral equation for the SW flux *S _{g}*(

*x*) that is a well-known starting point for computing global illumination in computer graphics (Cohen and Wallace 1993). There, it is commonly formulated in terms of the radiosity

*B*(

*x*) =

*α*(

*x*)

*S*(

_{g}*x*), which is defined as the hemispherically integrated radiant flux per unit area

*leaving*the surface in W m

^{−2}(Goral et al. 1984). The so-called radiosity equation can thus be written as

The first term on the rhs of (9) is given by *E*(*x*) = *α*(*x*)*S*_{sky}(*x*) as determined from (5) and the second term is *α*(*x*) times the terrain reflected radiation *S _{t}*(

*x*) received at

*x*.

In computer graphics terminology, according to (9) mountainous terrain is considered as part of an enclosure bounded by surrounding topography and by the sky, which is a large surface of zero reflectivity that can act as a radiation source.

In summary, we constructed our model to be applicable for microscale to mesoscale mountainous terrain with large albedos on clear-sky days when reflections and shading of remote parts of the surface dominate the effects of the atmosphere between the surfaces.

## 3. Numerical solution of the radiosity equation

### a. Finite element formulation of the radiosity equation

The numerical solution of (9) requires a discretization of a continuous topographic surface *z*(*x*). The discrete digital height model (DHM) is represented by an array of terrain elevations *z _{I}* =

*z*(

*x*) with

_{I}*I*= 1, 2 … ,

*N*in a domain of size

*L*=

_{x}*L*=

_{y}*L*with a given horizontal resolution Δ

*x*= Δ

*y*. This notation is helpful for the subsequent finite element (FE) method. DHM heights

*z*are allocated to the center of a surface or patch

_{I}*A*into which the topographic surface is decomposed. Flat rectangular planes are assumed with the surface normal vectors extracted from the DHM following Corripio (2002). The Euclidean distance between the centers of two patches

_{I}*A*and

_{I}*A*is denoted by

_{J}*r*

_{IJ}.

To solve (9), which is a Fredholm integral equation of the second kind, an FE method is traditionally applied [for an overview of different numerical methods see, e.g., Atkinson (1997)]. Cohen and Wallace (1993) present all steps necessary to discretize (9) using a low-order FE method with constant basis functions on *A _{I}*, leading to the linear system

for the discretization of Eq. (9) with an obvious notation for discrete quantities. The discretized integral kernel forms the so-called view factor matrix

where the view factors *F _{IJ}* (also known as form factors, configuration factors, etc.) are purely geometric factors. They are essential for the radiosity method and were originally introduced to understand the exchanges of radiant energy between surfaces in thermal engineering (Siegel and Howell 1978). The quantity

*F*equals the fraction of the total radiant power (flux) emitted by the surface patch

_{IJ}*A*that is received by

_{J}*A*. A derivation can be found in Goral et al. (1984), Cohen and Wallace (1993), or Encarnaçâo et al. (1996), among others.

_{I}To solve the linear system of (10) we must specify how *F _{IJ}* and

*χ*are computed from the discretized topographies.

_{IJ}### b. View factor computation

We compute the double area integral in (11) with a straightforward numerical integration scheme with a uniform but adaptive area subdivision (substructuring) of the finite patches *A _{I}*,

*A*in

_{J}*N*and

_{p}*N*small rectangular subpatches

_{q}*dA*and

_{P}*dA*, respectively. This can be regarded as a sort of “area sampling.” Summing the individual subplane view factors yields a numerical solution of (11):

_{Q}A threshold of substructuring for both patches is applied according to the suggestion of McCluney (1994) that a finite source can be approximated by a point source if its maximum lateral dimension is less than 10% of the distance to it.

Some numerical and analytical methods for the computation of the double area integral in (11) are given by Encarnaçâo et al. (1997). Analytical solutions for surfaces of specific shapes and orientations are tabulated by Siegel and Howell (1978).

### c. Mutual visibility/sun-or-shadow detection

The remaining computation of mutual visibility *χ _{IJ}* and sun-or-shadow

*χ*

_{sun,I}for the solution of (11) is accomplished in the same manner. This is possible because a patch

*A*is directly illuminated if the patch and the sun are mutually visible. Sun-or-shadow status changes with the position of the sun whereas mutual visibility is fixed for a given topography.

_{I}Here, the patches *A _{I}* and

*A*are defined to be mutually visible (i.e.,

_{J}*χ*= 1) if all elevations of

_{IJ}*intermediate*patches lie below the connecting line between the two patch centers. The

*intermediate*patch center heights

*z*are taken for those patches that are traversed by the projection of the connecting line onto the horizontal plane. The sequence of patches traversed by the projected line is computed similarly to the fast traversal algorithm for ray tracing of Amanatides and Woo (1987).

Sun or shadow is detected similarly; that is, *χ*_{sun,I} = 1 if all center heights of intermediate patches between *A _{I}* and the boundary patch

*A*, which is located in the direction of azimuth angle to the sun, are below the ray connecting the patch and the sun. This implies that the azimuthal resolution of the horizon line for each patch is given by the number of boundary patches. With this method, partial visibility and partial sun or shadow is not accounted for.

_{J}### d. Iterative solution of the radiosity equation

Despite the analogy between the SW radiation balance and the radiosity problem in computer graphics, the question remains whether common iteration methods also apply directly to the SW radiation balance problem. This is not obvious a priori because the geometries are different for the two systems (closed environments such as rooms versus complex topographies with the sky at the top as an area of zero reflectivity). Because the geometry determines the entries for the matrix of the linear system and thereby the convergence properties of the iteration methods, we have tested to see if common iteration methods can be applied.

Progressive refinement or progressive radiosity (PR), introduced by Cohen et al. (1988), was chosen here to solve the radiosity equation [(10)] for complex terrain. It is a fast iterative method to solve the system of *N* linear equations that does not require the storage of the whole *N* × *N* view factor matrix. The overall advantage of PR is a fast convergence rate produced by rearranging the order of considered patches according to their “relevance” for the iteration. This is done by applying the *shooting* criterion of Cohen et al. (1988), which computes the maximum unshot radiant power per iteration step *k*. The unshot radiosity at iteration start is equal to the emission and increases during an iteration by an amount . The rate of convergence of PR is significantly faster when compared to Gauss–Seidel (GS) or Jacobi (JA) iteration, which was explicitly checked for some test cases (cf. Helbig 2009). Note that in contrast JA and GS iteration are algorithms that *gather* the radiation of all other patches in the environment to a patch *A _{I}*. Therefore, one gathering step (i.e., one iteration step

*k*) updates the radiosity of only one patch. For more details about the performance of JA, GS, and PR for the radiosity problem, see, for example, Gortler et al. (1994).

Two extensions of the PR iteration were made for the SW radiation balance problem. First, we modified the shooting criterion of Cohen et al. (1988) to include the terrain view factor 1 − *F*_{sky,I} at each *x _{I}*; that is, we compute the maximum of . Thus, redundant shootings from mountain peaks where 1 −

*F*

_{sky,I}≈ 0 but where

*S*

_{sky,I}is large because of its exposure and the occurrence of large albedos (snow at the top) are avoided. Second, we developed a new stopping criterion to terminate the iterative solution. In the literature, heuristic stopping criteria are commonly used (e.g., Bu and Deprettere 1989). The drawback of such heuristic criteria is that the error ɛ

^{(k)}=

*B*−

*B*

^{(k)}at the

*k*th iteration step usually cannot be inferred from them. We therefore developed a stopping criterion based on the auxiliary quantity ɛ*

^{(k)}with that can be related to ɛ

^{(k)}by a bound ‖ɛ

^{(k)}‖

_{1}<

*Y*‖ɛ*

^{(k)}‖

_{1}. Here, ‖·‖

_{1}denotes the ℓ

_{1}-vector norm. The geometry-dependent factor

controls the tightness of the bound (cf. Helbig 2009 for more details). We terminate the iteration as soon as

Note that we define a stopping tolerance ɛ or iteration error according to the measurement accuracy of a standard net pyranometer (i.e., ±10%). Here, ɛ is chosen such that the error is 10% of the total reflectable radiation at iteration start Δ*B*^{(0)}, namely, ‖ɛ^{(k)}‖_{1} ≤ *Y*‖ɛ*^{(k)}‖_{1} ≤ *Yɛ* < 0.1‖Δ*B*^{(0)}‖_{1}.

Given this accuracy, the numerical solution of the radiosity Eq. (9) in model domains of mean slope angles from 10° to 30° with 100 × 100 grid cells requires approximately 1 to 3 min on an AMD Dual Core Opteron 270 with 2 GHz and 2 Gb RAM.

## 4. Random model topographies

Radiation transfer in mountainous terrain is a typical example of a physical process in a random environment. As such, it requires a sufficiently large number of study cases with similar terrain characteristics to study the influence of terrain on the SW radiation balance systematically. Instead of selecting many study cases from *real* mountain areas, we primarily investigated the radiosity equation on *simulated* topographies with well-defined statistical properties.

We start from the simplest picture and assume the existence of characteristic scales, namely a typical valley-to-peak elevation difference *σ* (height of mountains) and a typical lateral extension *ξ* (width of mountains) in system sizes *L* (cf. Fig. 2). We note that the statistical properties of real topographies might be more complex. The results in Dietler and Zhang (1992) and Weissel et al. (1994) indicate that DHM correlations of the Swiss Alps, Ethiopia, Saudi Arabia, and Somalia exhibit power law scaling, implying the lack of characteristic length scales.

The continuous radiosity Eq. (9) in turn requires acontinuous representation of model topographies, which we choose as Gaussian random fields with a Gaussian covariance. Our motivation for this choice is twofold: First, Gaussian random fields are probably the best investigated type of random field and as such they provide a common starting point for various problems involving physical processes in random environments [for an introduction to random fields, see e.g., Adler (1981)]. Because such an analysis has never been done before for the radiation problem, it is justified to start from simple Gaussian assumption. Second, from a mathematical point of view the radiosity equation [(9)] requires smooth or at least piecewise differentiable surfaces (which is guaranteed by the Gaussian covariance); thus, the influence of grid resolution Δ*x* can always be minimized with a sufficiently small Δ*x*.

For numerical purposes, the continuous Gaussian random field of elevations *z*(*x*) is then sampled at *N* grid positions *x _{I}*,

*I*= 1, 2, …

*N*, which constitute the (discrete) DHM as described in section 3. Thus, our topographies are random vectors

*z*=

_{I}*z*(

*x*),

_{I}*I*= 1, 2, …

*N*drawn from an

*N*-dimensional Gaussian distribution that is characterized by the covariance matrix

and mean *z*, which is taken to be homogeneous with = 2000 m for all *I* = 1, 2, … *N*. The overbar denotes ensemble averaging. Hence, the typical valley-to-peak elevation difference *σ* is given by the square root of the variance and the typical lateral extension is given by the correlation length *ξ* of the field of elevations, respectively. Note that all length scales are measured in meters. With this choice, the mean slope (tangent of the slope angle) of the DHM is equal to *σ*/*ξ* (Adler 1981). Following standard methods, correlated random vectors are generated from independent Gaussian random variables and a Cholesky factorization of the covariance matrix. (An example of a realization with *σ* = 180 m and *ξ* = 500 m is shown in Fig. 9a).

## 5. Results

In the following, we focus on the sensitivity of incident reflected radiation *S _{t}* on geometric parameters of the topography and solar elevation angles. First, however, we present some general information about all involved parameters. The sensitivity of

*S*was investigated with direct radiation of

_{t}*S*

_{b,I}= 1000 W m

^{−2}cos

*θ*

_{I}

*χ*

_{sun,I}, diffuse sky radiation of

*S*

_{d,I}= 150 W m

^{−2 }

*F*

_{sky,I}, a constant sun azimuth in the south direction, four sun elevation angles of 8°, 30°, 60° and 90°, two homogeneous albedo values

*α*=

_{I}*α*of 0.3 and 0.8, and three different mean slopes (10°, 20°, 30°) of the terrain. Each mean slope was realized by three different combinations of

*σ*and

*ξ*(given in Table 1). With 10 different realizations of a Gaussian topography per set of parameter values, this makes up a total of, for example, 720 simulations for the results in section 5a. Mean values of terrain reflected radiation

*S*were obtained from an ensemble average over the 10 realizations of the topography. In addition, we also spatially averaged

_{t}*S*within each realization. The standard deviations depend on the steepness of the terrain, from about 20% for gentle terrain up to 40% in steep terrain. We consider this as sufficient for this broad overview of parameter space. To investigate further details of the radiosity equation on Gaussian topographies, a greater number of realizations is required.

_{t}### a. Influence of mean slopes and solar elevation

The influence of varying terrain slopes and solar elevation angles on *S _{t}* is presented in Fig. 3 for an albedo of 0.3 and in Fig. 4 for an albedo of 0.8, for a fixed domain size

*L*= 2500 m and grid resolution Δ

*x*= 25 m.

The mean values of incident reflected radiation *S _{t}* increase with increasing sun elevation angles (

*x*axis), increasing mean slopes (line types), and increasing albedos (cf. Figs. 3 and 4). Note that two different

*y*-axis scales are applied to show the different dependencies of

*S*in Figs. 3 and 4. The line types in each slope group are arranged according to their ratios

_{t}*L*/

*ξ*. With an increasing

*L*/

*ξ*(i.e., more terrain is included in the model domain), the increase of

*S*with solar elevation is more pronounced. More

_{t}*S*is obtained with decreasing mountain–valley extension

_{t}*ξ*or with an increasing domain size

*L*. From these results, it is obvious that

*L*has to be much larger than

*ξ*to compute a precise

*S*(i.e.,

_{t}*L*≫

*ξ*). Remarkably large mean

*S*values are found for large albedos similar to those albedos of new snow surfaces. As a rule of thumb,

_{t}*mean*values of

*S*in Gaussian topographies vary between 5 and 70 W m

_{t}^{−2}for realistic conditions typical of snow-covered alpine terrain.

Furthermore, in Figs. 3 and 4 it is shown how the simple characterization of *S _{t}* in terms of two length scales compares to that obtained with a real DHM. The real DHM was taken around the Weissfluhjoch study site in the Eastern Alps, Switzerland (http://www.swisstopo.ch/en/) and has a mean slope of 25° in a domain with

*L*= 2500 m and Δ

*x*= 25 m. For both albedos the orders of magnitude obtained for the real DHM are similar to the Gaussian DHM. The realistic values best match those obtained by the combination of

*σ*= 180 m and

*ξ*= 500 m, corresponding to a mean slope of 20° with

*L*/

*ξ*= 5.00 (cf. dashed line with plus signs). From a fit of the height covariance of the real DHM to a Gaussian the parameters,

*σ*= 145 m and

*ξ*= 401 are obtained. Thus, the agreement between the real DHM and the particular Gaussian curve can be also quantitatively confirmed because it represents the best match of the length scales—

*L*/

*ξ*(6.2 compared to 5.00),

*σ*(145 compared to 180), and

*ξ*(401 compared to 500)—and also the mean slope

*σ*/

*ξ*(0.362 compared to 0.36).

### b. Influence of gridcell resolutions

The influence on mean *S _{t}* of grid resolution Δ

*x*is presented for a sun elevation of 60°, two albedos of 0.3 and 0.8, one combination of (

*σ*,

*ξ*) with

*σ*= 360 m and

*ξ*= 1000 m, and a fixed domain size

*L*= 2500 m. Three different resolutions were investigated: Δ

*x*= 25 m, leading to a 100 × 100 grid; Δ

*x*= 50 m leading to a 50 × 50 grid; and Δ

*x*= 100 m, leading to a 25 × 25 grid.

In Fig. 5 mean incident reflected radiation *S _{t}* is shown for varying ratios of

*ξ*/Δ

*x*. Although only a small range of resolutions Δ

*x*was investigated, it is apparent that almost no dependence of

*S*appears with Δ

_{t}*x*(below numerical accuracy). This supports the initial motivation to use differentiable Gaussian random fields as DHMs to avoid a misleading dependence on Δ

*x*as long as Δ

*x*≪

*ξ*. In general, the latter inequality is not always satisfied and care has to be taken.

### c. Influence of system size

The observed influence of different domain sizes *L* on *S _{t}* is further investigated and presented for one sun elevation of 60°, two albedos of 0.3 and 0.8, one combination of (

*σ*,

*ξ*) with

*σ*= 360 m and

*ξ*= 1000 m, and one fixed grid resolution Δ

*x*= 50 m. Three different domain sizes were investigated:

*L*= 5000 m, leading to a 100 × 100 grid;

*L*= 2500 m, leading to a 50 × 50 grid; and

*L*= 1000 m, leading to a 20 × 20 grid.

In Fig. 6 the mean incident reflected radiation *S _{t}* is shown for varying ratios of

*L*/

*ξ*. Larger values of

*S*are obtained with larger

_{t}*α*and larger

*L*. Although only a small range of domain sizes

*L*was investigated, a remarkable increase is observable with larger

*L*(at maximum 34 W m

^{−2}). Hence, for fixed

*ξ*only at very large

*L*, which is not attained in our case, a constant mean

*S*will be obtained. This is obvious because terrain reflections are a manifest

_{t}*nonlocal*effect. Note that the requirement

*L*≫

*ξ*must be carefully met in simulations; otherwise, the results may still contain significant effects of the finite domain size.

### d. Comparison with a common parameterization: Isotropic view factors

Commonly, terrain reflected or emitted radiation is included in models using a simple isotropic view factor parameterization (e.g., as applied by Müller and Scherer 2005). With this approach, terrain radiation is estimated as a single isotropic terrain contribution in terms of the sky view factor (by 1 − *F*_{sky}). We compare mean values of *S _{t}* derived by the radiosity approach (regarded as an anisotropic view factor approach) to those derived by the isotropic view factor approach according to

Note that the sky view factor is computed here from the anisotropic (terrain) view factor sum [cf. (11)]. In Fig. 7,*S _{t}* values are shown for an albedo of 0.8; one fixed lateral extension

*ξ*= 500 m, and three vertical extensions

*σ*of 90, 180, and 290 m (cf. Table 1); a domain size of

*L*= 2500 m; and a grid resolution of Δ

*x*= 25 m computed with anisotropic or isotropic view factors.

For all tested slopes, isotropic view factors lead to underestimated *S _{t}* because only single terrain reflected radiation is accounted for in the isotropic view factor approach. Hence, the largest differences occur for the steepest DHM and largest sun elevation angle for which multiple terrain reflections are most important.

For the spatial distribution of *S _{t}* the comparison between the two view factor approaches is more interesting. To highlight the spatial

*S*differences in Fig. 8, the mean of minimum and the mean of maximum differences are compared to mean differences and shown for

_{t}*α*= 0.8 and for all steep DHMs with mean slopes of 30° (cf. Table 1). Because there are always surfaces that do not see any other terrain [i.e., do not receive any

*S*(e.g., horizontal mountain peaks)] the mean of minimum differences is zero for all sun elevation angles. Whereas mean differences vary between 12 and 32 W m

_{t}^{−2}, the mean of maximum difference reaches remarkable values between 106 and 242 W m

^{−2}. The largest differences are obtained if the sun elevation angle is equal to the mean slope (i.e., 30°). This can be understood since the local values

*S*

_{sky,I}of sky radiation are used in Eq. (16) to estimate remote reflections for the isotropic view factor approach. In general, the largest differences compared to the anisotropic view factor approach are obtained during partial shading, which is most pronounced if the sun elevation is equal to the mean slope.

To visualize the spatial inhomogeneity of *S _{t}* between the two view factor approaches, the

*S*difference is shown in Fig. 9b for one Gaussian DHM (Fig. 9a). More asymmetric

_{t}*S*values are obtained with anisotropic view factors (i.e., multiple terrain reflections), leading to larger

_{t}*S*values on the east-facing and northwest-facing slopes. Additionally, about 50 W m

_{t}^{−2}more

*S*is obtained for anisotropic view factors in the center of the model domain (yellow/orange-colored range). The relative error between both approaches can easily reach 50%. The domain center is shaded by the small mountain located southward from it, leading to underestimated

_{t}*S*with the isotropic view factor approach on shaded patches. Note that (16) can also erroneously lead to overestimated

_{t}*S*on those patches that receive a lot of direct and diffuse sky radiation but see only little terrain.

_{t}### e. Effective albedo values of complex terrain

Finally, we briefly discuss an important application of reflected radiation, namely the effective albedo *α*_{eff} of complex terrain. It plays an important role in applications such as climate modeling, remote sensing, and the computation of the energy budget for complex sites (Kylling et al. 2000b; Schmucki et al. 2001; Weihs et al. 2001; Klok et al. 2003; Fortuniak 2008).

Here, we define an effective albedo as

which is the ratio of the mean outgoing radiation reflected back to the sky and the mean incident radiation received from the sky. This definition corresponds to a coarse-grained view of the topography in which the trapping of radiation is not recognized and a decreased outgoing radiation is observed. In this case, *α*_{eff} must be smaller than the domain average of *α*(*x*) for arbitrary albedo distributions. We note that this definition differs from those given in the articles cited above because appropriate definitions of *α*_{eff} are dictated by the measurement setup. A comparison of these is not the goal of the present work; we shall focus on the common physical effect inherent to many remote sensing methods.

In the following, we consider again the homogeneous case *α*(*x*) = *α*, which already reveals some interesting effects. In Fig. 10 the albedo difference *α* − *α*_{eff} is shown for a domain size of *L* = 2500 m, a grid resolution of Δ*x* = 25 m, and one combination of (*σ*, *ξ*) with *σ* = 290 m and *ξ* = 500 m (cf. Table 1). For larger sun elevations we find an almost constant or at most a slightly increasing difference. If the sun elevation drops below the mean slope (30°), shading becomes increasingly important and the effective albedo drops. The largest difference of 0.025 is obtained for an intermediate value of *α* ≈ 0.5.

The maximum difference at intermediate values can be explained by an approximation formula for the albedo difference. The radiosity equation for *S _{g}*(

*x*) can be solved exactly for small

*α*in terms of the well-known Neumann series expansion (Atkinson 1997). Retaining only linear terms in

*α*and averaging all factors separately leads to ≈ {1 +

*α*[1 − ]}, which has been well confirmed by our comparison with the isotropic view factor approach [cf. (16)]. By inserting this into (17), the albedo difference becomes

Indeed, the dependence on sun elevation dropped out by averaging all factors independently. However, (18) reveals the interplay between the domain-averaged albedo *α* and the topography via the mean sky view factor . This approximation is plotted in Fig. 11 together with the data from Fig. 10 averaged over all sun elevation angles. Note that no additional fitting is involved. The difference must vanish for all topographies at *α* = 0 and *α* = 1 where either all or no radiation is “trapped,” respectively. Obviously, the vanishing difference at *α* = 1 is not captured by the low albedo expansion (18), but the maximum difference at an intermediate value depending on is well predicted.

## 6. Discussion and conclusions

We have applied the radiosity approach to three-dimensional radiative transfer (RT) in complex terrain to systematically study the influence of terrain effects such as shading and anisotropic multiple terrain reflections on the shortwave radiation balance. The radiosity approach in an enclosure originally introduced in computer graphics or thermal engineering can be generalized to RT in mountains in a straightforward way. The numerical methods can easily be tailored to the open environment, and a progressive refinement iteration serves well to solve the finite element discretization.

In contrast to previous work, we started from the continuum description of RT. On the one hand, this reveals, for instance, that formulas commonly applied for the sky view factor (8) given in the literature can be grossly incorrect. In particular cases we found the error caused by such “approximations” to be comparable in magnitude to the terrain reflected radiation *S _{t}* itself. A systematic investigation of this source of error, however, is missing. On the other hand, the continuum description easily allows us to deduce approximation formulas such as (18).

Without loss of generality we used constant values of direct *S _{b}* and diffuse sky

*S*radiation and constant surface properties for our present study and allowed the incident radiation in (5) to vary only with solar elevation (including shading) and terrain geometry. By including a realistic parameterization for

_{d}*S*and

_{b}*S*, taking into account atmospheric conditions, the model can also be verified against real radiation measurements on clear-sky days (N. Helbig et al. 2009, unpublished manuscript). In general, we found mean values of

_{d}*S*between 5 and 75 W m

_{t}^{−2}on a set of random topographies. The values depend on sun elevation, the combination of (

*σ*,

*ξ*) and

*L*, and a domain-averaged albedo, which were chosen to cover typical summer and winter conditions in alpine terrain. Note that maximum values of

*S*can be significantly higher and a reliable estimate of

_{t}*S*always requires an appropriate choice of the additional length scales

_{t}*L*and Δ

*x*such that Δ

*x*≪

*ξ*≪

*L*. The requirement of large domain sizes

*L*coincides well with the demand for a larger

*L*in cloudy RTMs; backscattering from clouds has an impact from distances more than 10 km away (Degünther et al. 1998; Mayer and Degünther 2000).

The predictions of the widely used isotropic parameterization for *S _{t}* such as in (16) are found to be in agreement (mean differences up to 10 W m

^{−2}) with the radiosity approach as long as they are

*averaged*over the domain. However, because these parameterizations use either local [cf. (16)] or domain-averaged incident radiation values, predictions for spatial patterns for

*S*are poor when compared to the anisotropic radiosity approach.

_{t}As a consequence of terrain reflections, effective albedo differs from its domain-averaged value. We found that this effect is most pronounced at intermediate values of the domain averaged albedo because for high reflectivities the trapping of the radiation by topography is less pronounced. We also found a slight influence of sun elevation. Although the influence of terrain revealed by the albedo difference seems to be small, it is on the same order of magnitude as the terrain effects on an effective albedo found by Klok et al. (2003) and Krayenhoff and Voogt (2007). Furthermore, we expect that the albedo difference will exhibit a more complicated behavior for complex terrain with *heterogeneous* albedo distributions, which can easily be addressed by the radiosity approach. In contrast, Kylling et al. (2000b) found larger albedo differences on the order of 0.2 for Tromsø (Norway) if a large fraction of open water was included in their three-dimensional MC modeling. They note that the decrease is wavelength dependent and that further investigations are required, especially for the influence of topography. Similar results were found for the Davos region by Schmucki et al. (2001) from local UV measurements combined with satellite data.

Our general characterization of topography in terms of Gaussian random fields can indeed be questioned, for example because exponential slope distributions found for real DHMs (Essery and Marks 2007) contradict the Gaussian assumption. However, we found that the terrain reflected radiation of a real, non-Gaussian topography fits nicely into this characterization in terms of the fundamental length scales. This can be partly understood since (Essery and Marks 2007) have shown that typical properties of radiation in complex terrain are predominantly determined by extreme value characteristics of the topography. Thus, a certain degree of universality can be expected, although a more careful study of this issue is necessary. However, if the influence of terrain can be effectively described by Gaussian random fields as disorder for the continuum radiosity equation, this well-defined mathematical problem is most likely amenable to approximations that might replace commonly used ad hoc parameterizations.

Note that certain extensions to the radiosity equation [(9)] can easily be made without significantly increasing the numerical complexity of the model. First, including the transmittance function *τ* [cf. (3)] for a plane-parallel atmosphere leads to at most one additional numerical integration for each pair of patches. This would be an improvement because with increasing distance between two points, the assumption of the negligible atmosphere between them [*Q* = 0, *τ* = 1 in (3)] will break down and the influence of the atmosphere in between will dominate. Second, for the incident sky flux [cf. (5)] the assumption of Lambertian reflections can be abandoned and a full anisotropic bidirectional reflectance distribution function (BRDF) could be incorporated. Nevertheless, Lambertian reflections are a reasonable assumption for snow-covered surfaces as found by measured hemispherical directional reflectance factors (HDRFs) with a goniometer by Odermatt et al. (2005) over snow in Davos (Switzerland) and with a gonio-spectrometer by Bourgeois et al.(2006) in Greenland. Thereafter, the HDRF of snow-covered surfaces at noon is nearly isotropic and anisotropy increases only with increasing solar zenith angles (larger forward scattering) and increasing wavelength. In general, however, larger solar zenith angles lead to increased shading and hence decreased terrain reflections in complex terrain such that the anisotropy can be neglected as a first approximation. Note that the HDRF is a nondimensional equivalent to the BRDF, which considers hemispherically integrated incident and directional outgoing radiance. The BRDF can be obtained from the HDRF if the atmospheric composition is available to correct for the diffuse radiation (Bourgeois et al. 2006).

In conclusion, we suggest that our approach is applicable from microscale to mesoscale mountainous terrain such as typical subgrid scales within, for example, numerical weather prediction models. For such domains the radiosity approach is well suited to investigate spatial distributions of radiation as an effect of heterogeneous albedo distributions and its interplay with topography. We believe that our approach will help to connect simple and sophisticated approaches from a conceptual point of view.

## Acknowledgments

We gratefully acknowledge the anonymous referees for their useful suggestions and careful reading of the manuscript. This work was co-financed by the CTI, which is the Confederation’s innovation promotion agency (Project 7723.1 ESPP-ES) and by the European Commission via AWARE.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Nora Helbig, WSL Institute for Snow and Avalanche Research, SLF, Flüelastr 11, 7260 Davos Dorf, Switzerland. Email: helbig@slf.ch