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, Fsky. 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

 
formula

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) = Psky(x) + Pt(x)—and assume that Psky(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.

Fig. 1.

Schematic of the cross section of a valley to sketch the RT between two points x and x′.

Fig. 1.

Schematic of the cross section of a valley to sketch the RT between two points x and x′.

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′)/|xx′|2. Then (1) can be written as

 
formula

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

 
formula

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 R0(x′, x). Furthermore, we assume that all surface sources are isotropically emitting; that is,

 
formula

The set of Eqs. (2), (3), and (4) provides our general starting point for studying the radiation balance in complex terrain. By specifying Psky, R0, Q, and τ, the model can be applied to all wavelengths under the assumptions made above.

b. Radiosity equation for shortwave radiation

In the following, we concentrate on the incident broadband global radiation flux P(x) = Sg(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′) = R0(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 R0(x′) is related to the incident SW flux Sg(x′) at x′ via R0(x′) = α(x′)Sg(x′)/π. Lambertian surfaces are a common assumption in RTMs.

  • (iii) We assume that downwelling radiance from the sky R[x, Ω(x)] consists solely of a diffuse contribution that is isotropic and a direct solar beam contribution involving a delta function that is valid for a point-shaped sun. Hence, the sky flux can be written as 
    formula

Here, Sb denotes the direct solar beam flux and Sd the diffuse flux in W m−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

 
formula

in the integration of the isotropic diffuse component. By defining ϑh[ϕ(x)] as the horizon angle in the sloped coordinate system, it follows that

 
formula

or, equivalently,

 
formula

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 Sg(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)Sg(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

 
formula

The first term on the rhs of (9) is given by E(x) = α(x)Ssky(x) as determined from (5) and the second term is α(x) times the terrain reflected radiation St(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 zI = z(xI) with I = 1, 2 … , N in a domain of size Lx = Ly = L with a given horizontal resolution Δx = Δy. This notation is helpful for the subsequent finite element (FE) method. DHM heights zI are allocated to the center of a surface or patch AI 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 AI and AJ is denoted by rIJ.

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 AI, leading to the linear system

 
formula

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

 
formula

where the view factors FIJ (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 FIJ equals the fraction of the total radiant power (flux) emitted by the surface patch AJ that is received by AI. A derivation can be found in Goral et al. (1984), Cohen and Wallace (1993), or Encarnaçâo et al. (1996), among others.

To solve the linear system of (10) we must specify how FIJ and χIJ are computed from the discretized topographies.

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 AI, AJ in Np and Nq small rectangular subpatches dAP and dAQ, respectively. This can be regarded as a sort of “area sampling.” Summing the individual subplane view factors yields a numerical solution of (11):

 
formula

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 AI 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.

Here, the patches AI and AJ are defined to be mutually visible (i.e., χIJ = 1) if all elevations of 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 AI and the boundary patch AJ, 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.

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 AI. 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 − Fsky,I at each xI; that is, we compute the maximum of . Thus, redundant shootings from mountain peaks where 1 − Fsky,I ≈ 0 but where Ssky,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) = BB(k) at the kth 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

 
formula

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

 
formula

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)1Y‖ɛ*(k)1 < 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.

Fig. 2.

Schematic of a transect of a random topography with the characteristic scales, namely the valley-to-peak elevation difference σ and a typical lateral extension ξ of the valleys or mountains in a system of size L. The mean elevation is characterized by z.

Fig. 2.

Schematic of a transect of a random topography with the characteristic scales, namely the valley-to-peak elevation difference σ and a typical lateral extension ξ of the valleys or mountains in a system of size L. The mean elevation is characterized by z.

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 xI, I = 1, 2, … N, which constitute the (discrete) DHM as described in section 3. Thus, our topographies are random vectors zI = z(xI), I = 1, 2, … N drawn from an N-dimensional Gaussian distribution that is characterized by the covariance matrix

 
formula

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).

Fig. 9.

(a) Example of a Gaussian DHM with a mean slope of 20° with ξ = 500 m, σ = 180 m, L = 2500 m, and a gridcell resolution Δx = 25 m. (b) Difference between anisotropically derived St (radiosity approach) and isotropically derived St [cf. (16)]. An albedo value of 0.8 and a sun elevation of 60° were chosen for this simulation.

Fig. 9.

(a) Example of a Gaussian DHM with a mean slope of 20° with ξ = 500 m, σ = 180 m, L = 2500 m, and a gridcell resolution Δx = 25 m. (b) Difference between anisotropically derived St (radiosity approach) and isotropically derived St [cf. (16)]. An albedo value of 0.8 and a sun elevation of 60° were chosen for this simulation.

5. Results

In the following, we focus on the sensitivity of incident reflected radiation St on geometric parameters of the topography and solar elevation angles. First, however, we present some general information about all involved parameters. The sensitivity of St was investigated with direct radiation of Sb,I = 1000 W m−2 cosθIχsun,I, diffuse sky radiation of Sd,I = 150 W m−2 Fsky,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 St were obtained from an ensemble average over the 10 realizations of the topography. In addition, we also spatially averaged St 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.

Table 1.

Each of the mean slope angles is realized by three different combinations of (σ, ξ). The mean slope is computed from tan−1(σ/ξ); the numbers in the table are the corresponding σ values in meters.

Each of the mean slope angles is realized by three different combinations of (σ, ξ). The mean slope is computed from tan−1(σ/ξ); the numbers in the table are the corresponding σ values in meters.
Each of the mean slope angles is realized by three different combinations of (σ, ξ). The mean slope is computed from tan−1(σ/ξ); the numbers in the table are the corresponding σ values in meters.

a. Influence of mean slopes and solar elevation

The influence of varying terrain slopes and solar elevation angles on St 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.

Fig. 3.

Mean St as a function of sun elevation angle for a domain size of L = 2500 m and albedo of α = 0.3. The legend defines the curves in terms of visible terrain given by L/ξ per combination of (σ, ξ). The solid lines denote a mean slope of 10°, the dashed lines a mean slope of 20°, and the dotted lines a mean slope of 30°. The real topography example with a mean slope of 25° (DHM around Weissfluhjoch study site) is denoted by the dashed line with squares.

Fig. 3.

Mean St as a function of sun elevation angle for a domain size of L = 2500 m and albedo of α = 0.3. The legend defines the curves in terms of visible terrain given by L/ξ per combination of (σ, ξ). The solid lines denote a mean slope of 10°, the dashed lines a mean slope of 20°, and the dotted lines a mean slope of 30°. The real topography example with a mean slope of 25° (DHM around Weissfluhjoch study site) is denoted by the dashed line with squares.

Fig. 4.

As in Fig. 3, but for an albedo of α = 0.8. Note the different y-axis scale compared to Fig. 3.

Fig. 4.

As in Fig. 3, but for an albedo of α = 0.8. Note the different y-axis scale compared to Fig. 3.

The mean values of incident reflected radiation St 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 St in Figs. 3 and 4. The line types in each slope group are arranged according to their ratios L/ξ. With an increasing L/ξ (i.e., more terrain is included in the model domain), the increase of St with solar elevation is more pronounced. More St is obtained with decreasing mountain–valley extension ξ 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 St (i.e., Lξ). Remarkably large mean St values are found for large albedos similar to those albedos of new snow surfaces. As a rule of thumb, mean values of St in Gaussian topographies vary between 5 and 70 W m−2 for realistic conditions typical of snow-covered alpine terrain.

Furthermore, in Figs. 3 and 4 it is shown how the simple characterization of St 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 St 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 St 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 St appears with Δ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.

Fig. 5.

Mean St as a function of ξx for lateral extension ξ = 1000 m, domain size L = 2500 m, and mean slope 20°.

Fig. 5.

Mean St as a function of ξx for lateral extension ξ = 1000 m, domain size L = 2500 m, and mean slope 20°.

c. Influence of system size

The observed influence of different domain sizes L on St 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 St is shown for varying ratios of L/ξ. Larger values of St are obtained with larger α 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 St will be obtained. This is obvious because terrain reflections are a manifest 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.

Fig. 6.

Mean St as a function of L/ξ for grid resolution Δx = 50 m, lateral extension ξ = 1000 m, and mean slope 20°.

Fig. 6.

Mean St as a function of L/ξ for grid resolution Δx = 50 m, lateral extension ξ = 1000 m, and mean slope 20°.

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 − Fsky). We compare mean values of St derived by the radiosity approach (regarded as an anisotropic view factor approach) to those derived by the isotropic view factor approach according to

 
formula

Note that the sky view factor is computed here from the anisotropic (terrain) view factor sum [cf. (11)]. In Fig. 7,St 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.

Fig. 7.

Mean St as a function of sun elevation angle for three mean slopes, albedo α = 0.8, and anisotropic and isotropic view factors. Solid lines denote a mean slope of 10°; dashed lines, 20°; and dotted lines, 30°. A lateral extension ξ = 500 m was used.

Fig. 7.

Mean St as a function of sun elevation angle for three mean slopes, albedo α = 0.8, and anisotropic and isotropic view factors. Solid lines denote a mean slope of 10°; dashed lines, 20°; and dotted lines, 30°. A lateral extension ξ = 500 m was used.

For all tested slopes, isotropic view factors lead to underestimated St 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 St the comparison between the two view factor approaches is more interesting. To highlight the spatial St differences in Fig. 8, the mean of minimum and the mean of maximum differences are compared to mean differences and shown for α = 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 St (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−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 Ssky,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.

Fig. 8.

Mean of minimum difference, mean of maximum difference, and mean difference St between anisotropic and isotropic view factors for albedo α = 0.8 and for the dotted lines in Fig. 6 with a mean slope of 30° and lateral extension ξ = 500 m.

Fig. 8.

Mean of minimum difference, mean of maximum difference, and mean difference St between anisotropic and isotropic view factors for albedo α = 0.8 and for the dotted lines in Fig. 6 with a mean slope of 30° and lateral extension ξ = 500 m.

To visualize the spatial inhomogeneity of St between the two view factor approaches, the St difference is shown in Fig. 9b for one Gaussian DHM (Fig. 9a). More asymmetric St values are obtained with anisotropic view factors (i.e., multiple terrain reflections), leading to larger St values on the east-facing and northwest-facing slopes. Additionally, about 50 W m−2 more St 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 St with the isotropic view factor approach on shaded patches. Note that (16) can also erroneously lead to overestimated St on those patches that receive a lot of direct and diffuse sky radiation but see only little terrain.

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

 
formula

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.

Fig. 10.

Differences between domain-averaged albedo α and effective albedo αeff are shown as a function of sun elevation for mean slope angle of 30° realized by σ = 290 m and ξ = 500 m. The gridcell resolution is Δx = 25 m.

Fig. 10.

Differences between domain-averaged albedo α and effective albedo αeff are shown as a function of sun elevation for mean slope angle of 30° realized by σ = 290 m and ξ = 500 m. The gridcell resolution is Δx = 25 m.

The maximum difference at intermediate values can be explained by an approximation formula for the albedo difference. The radiosity equation for Sg(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

 
formula

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.

Fig. 11.

Differences between domain-averaged albedo α and effective albedo αeff from Fig. 10 are additionally averaged over all four sun elevation angles and shown as a function of albedo for mean slope angle of 30° with σ = 290 m and ξ = 500 m. The gridcell resolution is Δx = 25 m. The dotted line is the approximation formula [cf. (18)].

Fig. 11.

Differences between domain-averaged albedo α and effective albedo αeff from Fig. 10 are additionally averaged over all four sun elevation angles and shown as a function of albedo for mean slope angle of 30° with σ = 290 m and ξ = 500 m. The gridcell resolution is Δx = 25 m. The dotted line is the approximation formula [cf. (18)].

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 St 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 Sb and diffuse sky Sd 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 Sb and Sd, 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 St between 5 and 75 W m−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 St can be significantly higher and a reliable estimate of St always requires an appropriate choice of the additional length scales 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 St 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 St are poor when compared to the anisotropic radiosity approach.

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

REFERENCES
Adler
,
R. J.
,
1981
:
The Geometry of Random Fields.
Wiley, 280 pp
.
Amanatides
,
J.
, and
A.
Woo
,
1987
:
A fast voxel traversal algorithm for ray tracing.
Proc. Eurographics ’87, Amsterdam, Netherlands, European Association for Computer Graphics, 3–10
.
Atkinson
,
K. E.
,
1997
:
The Numerical Solution of Integral Equations of the Second Kind.
Cambridge University Press, 584 pp
.
Bellasio
,
R.
,
G.
Maffeis
,
J. S.
Scire
,
M. G.
Longoni
,
R.
Bianconi
, and
N.
Quaranta
,
2005
:
Algorithms to account for topographic shading effects and surface temperature dependence on terrain elevation in diagnostic meteorological models.
Bound.-Layer Meteor.
,
114
,
595
614
.
Borel
,
J-L.
,
Y.
Durand
,
F.
Sillion
, and
C.
Le Gal
,
1998
:
Computing the radiative balance in mountainous areas: Modelling, simulation and visualization.
ERCIM News, No. 34, European Research Consortium for Informatics and Mathematics, Sophia-Antipolis, France, 24–26. [Available online at http://www.ercim.org/publication/Ercim_News/enw34/borel.html]
.
Bourgeois
,
C. S.
,
P.
Calanca
, and
A.
Ohmura
,
2006
:
A field study of the hemispherical directional reflectance factor and spectral albedo of dry snow.
J. Geophys. Res.
,
111
,
D20108
.
doi:10.1029/2006JD007296
.
Bu
,
J.
, and
E. F.
Deprettere
,
1989
:
A VLSI system architecture for high-speed radiative transfer 3D image synthesis.
Visual Comput.
,
5
,
121
133
.
Cahalan
,
R. F.
, and
Coauthors
,
2005
:
The I3RC: Bringing together the most advanced radiative transfer tools for cloudy atmospheres.
Bull. Amer. Meteor. Soc.
,
86
,
1275
1293
.
Chen
,
Y.
,
A.
Hall
, and
K. N.
Liou
,
2006
:
Application of three-dimensional solar radiative transfer to mountains.
J. Geophys. Res.
,
111
,
D21111
.
doi:10.1029/2006JD007163
.
Cohen
,
M. F.
, and
J. R.
Wallace
,
1993
:
Radiosity and Realistic Image Synthesis.
Academic Press Professional, 381 pp
.
Cohen
,
M. F.
,
S. E.
Chen
,
J. R.
Wallace
, and
D. P.
Greenberg
,
1988
:
A progressive refinement approach to fast radiosity image generation.
Comput. Graphics
,
22
,
75
84
.
Corripio
,
J. G.
,
2002
:
Modelling the energy balance of high altitude glacerised basins in the central Andes. Ph.D. thesis, University of Edinburgh, 151 pp
.
Degünther
,
M.
,
R.
Meerkötter
,
A.
Albold
, and
G.
Seckmeyer
,
1998
:
Case study on the influence of inhomogeneous surface albedo on UV irradiance.
Geophys. Res. Lett.
,
25
,
3587
3590
.
Dietler
,
G.
, and
Y-C.
Zhang
,
1992
:
Fractal aspects of the Swiss landscape.
Physica A
,
191
,
213
219
.
Encarnaçâo
,
J.
,
W.
Strasser
, and
R.
Klein
,
1996
:
Gerätetechnik, Programmierung und Anwendung graphischer Systeme.
Vol. 1, Graphische Datenverarbeitung, Oldenborg, 418 pp
.
Encarnaçâo
,
J.
,
W.
Strasser
, and
R.
Klein
,
1997
:
Modellierung komplexer Objekte und photorealistische Bilderzeugung.
Vol. 2, Graphische Datenverarbeitung, Oldenborg, 427 pp
.
Erbs
,
D. G.
,
S. A.
Klein
, and
J. A.
Duffie
,
1982
:
Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation.
Sol. Energy
,
28
,
293
302
.
Essery
,
R.
, and
D.
Marks
,
2007
:
Scaling and parametrization of clear-sky solar radiation over complex topography.
J. Geophys. Res.
,
112
,
D10122
.
doi:10.1029/2006JD007650
.
Evans
,
K. F.
,
1998
:
The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer.
J. Atmos. Sci.
,
55
,
429
446
.
Fortuniak
,
K.
,
2008
:
Numerical estimation of the effective albedo of an urban canyon.
Theor. Appl. Climatol.
,
91
,
245
258
.
Gastellu-Etchegorry
,
J. P.
,
V.
Demarez
,
V.
Pinel
, and
F.
Zagolski
,
1996
:
Modeling radiative transfer in heterogeneous 3-D vegetation canopies.
Remote Sens. Environ.
,
58
,
131
156
.
Gastellu-Etchegorry
,
J. P.
,
E.
Martin
, and
F.
Gascon
,
2004
:
DART: A 3D model for simulating satellite images and studying surface radiation budget.
Int. J. Remote Sens.
,
25
,
73
96
.
Goral
,
C. M.
,
K. E.
Torrance
,
D. P.
Greenberg
, and
B.
Battaile
,
1984
:
Modeling the interaction of light between diffuse surfaces.
Comput. Graphics
,
18
,
213
222
.
Gortler
,
S.
,
M. F.
Cohen
, and
P.
Slusallek
,
1994
:
Radiosity and relaxation methods.
IEEE Comput. Graph. Appl.
,
14
,
48
58
.
Helbig
,
N.
,
2009
:
Application of the radiosity approach to the radiation balance in complex terrain. Ph.D. thesis, University of Zurich, 213 pp
.
Hock
,
R.
, and
B.
Holmgren
,
2005
:
A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden.
J. Glaciol.
,
51
,
25
36
.
Iqbal
,
M.
,
1983
:
An Introduction to Solar Radiation.
Academic Press, 390 pp
.
Klok
,
E. J.
,
W.
Greuell
, and
J.
Oerlemans
,
2003
:
Temporal and spatial variation of the surface albedo of Morteratschgletscher, Switzerland, as derived from 12 Landsat images.
J. Glaciol.
,
49
,
491
502
.
Krayenhoff
,
E. S.
, and
J. A.
Voogt
,
2007
:
A microscale three-dimensional urban energy balance model for studying surface temperatures.
Bound.-Layer Meteor.
,
123
,
433
461
.
Kylling
,
A.
,
A.
Dahlback
, and
B.
Mayer
,
2000a
:
The effect of clouds and surface albedo on UV irradiances at a high-latitude site.
Geophys. Res. Lett.
,
27
,
1411
1414
.
Kylling
,
A.
,
T.
Persen
,
B.
Mayer
, and
T.
Svenøe
,
2000b
:
Determination of an effective spectral surface albedo from ground-based global and direct UV irradiance measurements.
J. Geophys. Res.
,
105
,
4949
4959
.
Lehning
,
M.
,
I.
Völksch
,
D.
Gustafsson
,
T.
Nguyen
,
M.
Stähli
, and
M.
Zappa
,
2006
:
ALPINE3D: A detailed model of mountain surface processes and its application to snow hydrology.
Hydrol. Processes
,
20
,
2111
2128
.
Lehning
,
M.
,
H.
Löwe
,
M.
Ryser
, and
N.
Raderschall
,
2008
:
Inhomogeneous precipitation distribution and snow transport in steep terrain.
Water Resour. Res.
,
44
,
W07404
.
doi:10.1029/2007WR006545
.
Letsinger
,
S. L.
, and
G. A.
Olyphant
,
2007
:
Distributed energy-balance modeling of snow-cover evolution and melt in rugged terrain: Tobacco Root Mountains, Montana, USA.
J. Hydrol.
,
336
,
48
60
.
Liou
,
K. N.
,
2002
:
An Introduction to Atmospheric Radiation.
Academic Press, 583 pp
.
Liu
,
B. Y. H.
, and
R. C.
Jordan
,
1960
:
The interrelationship and characteristic distribution of direct, diffuse and total solar radiation.
Sol. Energy
,
4
,
1
19
.
Marchuk
,
G.
,
G.
Mikhailov
,
M.
Nazaraliev
,
R.
Darbinjan
,
B.
Kargin
, and
B.
Elepov
,
1980
:
The Monte Carlo Methods in Atmospheric Optics.
Springer-Verlag, 208 pp
.
Mayer
,
B.
,
2000
:
I3RC phase 2 results from the MYSTIC Monte Carlo model.
Abstracts, I3RC workshop, Tucson, AZ, International Radiation Commission. [Available online at http://www.bmayer.de/papers/i3rc2\_1.pdf]
.
Mayer
,
B.
, and
M.
Degünther
,
2000
:
Comment on “Measurements of erythemal irradiance near Davis Station, Antarctica: Effect of inhomogeneous surface albedo”.
Geophys. Res. Lett.
,
27
,
3489
3490
.
McCluney
,
W. R.
,
1994
:
Introduction to Radiometry and Photometry.
Artech House, 424 pp
.
Moore
,
I. D.
,
T. W.
Norton
, and
J. E.
Williams
,
1993
:
Modelling environmental heteorogenity in forested landscapes.
J. Hydrol.
,
150
,
717
747
.
Müller
,
M. D.
, and
D.
Scherer
,
2005
:
A grid- and subgrid-scale radiation parameterization of topographic effects for mesoscale weather forecast models.
Mon. Wea. Rev.
,
133
,
1431
1442
.
Nunez
,
M.
,
1980
:
The calculation of solar and net radiation in mountainous terrain.
J. Biogeogr.
,
7
,
173
186
.
Odermatt
,
D.
,
D.
Schläpfer
,
M.
Lehning
,
M.
Schwikowski
,
M.
Kneubühler
, and
K. I.
Itten
,
2005
:
Seasonal study of directional reflectance properties of snow.
EARSeL eProceedings
,
4
,
203
214
.
Oliphant
,
A. J.
,
R. A.
Spronken-Smith
,
A. P.
Sturman
, and
I. F.
Owens
,
2003
:
Spatial variability of surface radiation fluxes in mountainous terrain.
J. Appl. Meteor.
,
42
,
113
128
.
Olyphant
,
G. A.
,
1986
:
Longwave radiation in mountainous areas and its influence on the energy balance of alpine snowfields.
Water Resour. Res.
,
22
,
62
66
.
Orgill
,
J. F.
, and
K. G.
Hollands
,
1977
:
Correlation equation for hourly diffuse radiation on a horizontal surface.
Sol. Energy
,
19
,
357
359
.
Reindl
,
D. T.
,
W. A.
Beckman
, and
J. A.
Duffie
,
1990
:
Diffuse fraction correlations.
Sol. Energy
,
45
,
1
7
.
Schmucki
,
D.
,
S.
Voigt
,
R.
Phillipona
,
C.
Fröhlich
,
J.
Lenoble
,
A.
Ohmura
, and
C.
Wehrli
,
2001
:
Effective albedo derived from UV measurements in the Swiss Alps.
J. Geophys. Res.
,
106
,
5369
5383
.
Siegel
,
R.
, and
J. R.
Howell
,
1978
:
Thermal Radiation Heat Transfer.
Hemisphere Publishing, 814 pp
.
Siegel
,
R.
,
J. R.
Howell
, and
J.
Lohrengel
,
1991
:
Strahlungsaustausch zwischen Oberflächen und in Umhüllungen.
Vol. 2, Wärmeübertragung durch Strahlung, Springer-Verlag, 353 pp
.
Verseghy
,
D. L.
, and
D. S.
Munro
,
1989
:
Sensitivity studies on the calculation of the radiation balance of urban surfaces: I. Shortwave radiation.
Bound.-Layer Meteor.
,
46
,
309
331
.
Wang
,
Q.
,
J.
Tenhunen
,
M.
Schmidt
, and
O.
Kolcun
,
2006
:
A model to estimate global radiation in complex terrain.
Bound.-Layer Meteor.
,
119
,
409
429
.
Weihs
,
P.
, and
Coauthors
,
2001
:
Modeling the effect of an inhomogeneous surface albedo on incident UV radiation in mountainous terrain: Determination of an effective surface albedo.
Geophys. Res. Lett.
,
28
,
3111
3114
.
Weissel
,
J. K.
,
L. F.
Pratson
, and
A.
Malinverno
,
1994
:
The length-scaling properties of topography.
J. Geophys. Res.
,
99
,
13997
14012
.

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