Surface spectral irradiance and irradiance partitioning in a complex mountain environment: understanding location-dependent topographic effects in satellite imagery

Abstract

Remote sensing of mountain environments is extremely difficult due to atmosphere-topography-landcover couplings that govern the irradiant fluxes. Researchers have not adequately accounted for topographic effects in terms of irradiance components over the landscape. Consequently, we evaluated surface spectral-irradiance components and irradiance partitioning with respect to anisotropic-reflectance correction (ARC). Our work focuses on addressing issues of scale, parameter sensitivity analysis and irradiance partitioning to provide new insights into understanding the complex nature of topographic modulation of the radiation-transfer cascade. We accomplish this by characterizing irradiance components in numerical simulations that account for topographic effects and couplings. Our results reveal atmosphere-topographic couplings associated with irradiance components that are not spatially coincident. Furthermore, we found that commonly utilized adjacent-terrain irradiance parameterization schemes produce different results compared to a more comprehensive parameterization scheme. Parameter sensitivity analysis revealed that the exclusion of topographic effects also produces different results.

Keywords:

• Adjacent-terrain irradiance
• anisotropic-reflectance correction
• atmosphere-topography-landcover coupling
• radiation-transfer cascade
• surface irradiance

1. Introduction

It is well known that multi-scale topographic effects govern the short-wave irradiance fluxes that govern the magnitude of surface biophysical parameters and processes (Gratton et al. 1994; Huang et al. 2019; Chu et al. 2021; Bishop et al. 2022). The topic of multi-scale topographic effects is complex and encapsulates spatial concepts that formalize the operational scale-dependencies of radiation-transfer processes. Local, directional, distance and area dependence of matter-energy interactions involving terrain, atmosphere and land cover variations must be accounted for to accurately estimate irradiant fluxes. The concept also refers to variations in scale-dependence of processes that uniquely governs each irradiance component and the coupling of numerous radiation-transfer parameters. Furthermore, research has demonstrated that spatially-dependent coupling is not uniform for different irradiance components (Bishop et al. 2022). Consequently, it is difficult to accurately predict surface irradiance and surface biophysical parameters and processes, especially in mountain environments (Chi et al. 2022; Ma et al. 2021; Jasrotia et al. 2022). Unfortunately, we have not adequately accounted for scale-dependent integrative coupling of radiation-transfer and topographic parameters, given complex topographic effects and the anisotropic coupling of atmosphere-topography-landcover conditions (Oliphant et al. 2003; Bishop et al. 2022). Furthermore, there are numerous issues associated with numerical modeling and anisotropic-reflectance correction (ARC) approaches that have not been adequately addressed, which require advances in numerical modeling efforts (Bishop et al. 2019, 2022).

From a numerical modeling perspective, spatio-temporal dependencies associated with radiation-transfer parameters need to be accounted for to generate accurate estimates of surface-irradiance components. Investigators have recognized the need for: 1) characterizing spatially complex topography (Wang et al. 2005); 2) understanding the significance of various topographic effects (Li et al. 2016); 3) assessing the role of topographic effects and atmosphere-topography couplings (Schulmann et al. 2015; Bishop et al. 2019, 2022); 4) development and evaluation of improved irradiance parameterization schemes (Lai et al. 2010; Chu et al. 2021); 5) characterizing uncertainty in numerical modeling efforts (Pinty and Verstraete 1991; Wu et al. 2018); and 6) understanding the nature of irradiance partitioning from environmental, temporal, spatial and spectral perspectives (Bishop et al. 2019). Collectively, these issues bring about significant uncertainty in characterizing topographic effects that govern irradiance and reflectance variations that are ultimately recorded by satellite sensors, as common assumptions about topographic effects, significance and nature of radiation-transfer parameters (RTP), and irradiance parameterization schemes can generate highly generalized results.

Although the direct irradiance (E_b) is perhaps the easiest component to estimate, time-dependent parameters are frequently not accounted for (Bishop et al. 2019, 2022). Similarly, assumptions about the diffuse-skylight irradiance (E_d), involving its anisotropic nature and the anisotropic nature of topographic shielding, do not adequately account for atmosphere-topography coupling (Bishop et al. 2022). Similarly, it is generally assumed that the adjacent-terrain irradiance (E_t) does not contribute significant energy, and therefore, may not be accounted for given the computational complexity associated with its estimation. Nevertheless, Chu et al. (2021) and Gao et al. (2016) indicate that there is scientific consensus that E_t is important and is non-negligible. This may especially be the case given specific environmental settings (i.e. mountains) that exhibit a range of highly reflective features, extreme relief, and a strong terrain-orientation fabric, such that at various times of the day, numerous portions of the landscape may be dominated by E_d and/or E_t.

Many of the aforementioned issues also govern effective ARC of satellite imagery. Satellite imagery acquired over mountain environments must be radiometrically calibrated to remove coupled atmosphere-topography effects, so that image spectral variation can be effectively used for quantitative biophysical estimation of surface parameters and thematic-mapping (Richter et al. 2009; Vanonckelen et al. 2013; Li et al. 2016; Bishop et al. 2022). Removing topographic effects from imagery is notoriously difficult, and must account for surface irradiance variations and the Bidirectional Reflectance Distribution Function (BRDF) that encapsulates topographic effects and land-cover structure across a range of distances and directions (Pinty and Verstraete 1991; Schaaf et al. 1994; Wen et al. 2009; Li et al. 2012; Bishop et al. 2019, 2022). This may be attempted through numerical modeling efforts that account for these factors, or by using empirical ARC methods that attempt to characterize the spatial complexity of topography. Unfortunately, both approaches may not adequately characterize various topographic effects and couplings, given inadequate parameterization schemes that do not account for important RTP. Such is the case for the more complicated E_d and E_t irradiance components, as spatial-dependent interactions must be taken into consideration (Bishop et al. 2022). Perhaps one of the most difficult parameters to account for is E_t, as atmosphere-topography-landcover couplings need to be formalized.

The objectives of this research, therefore, are to assess topographic effects, couplings and irradiance partitioning over the Nanga Parbat Massif in Pakistan. The Nanga Parbat Massif is an ideal location for evaluating surface irradiance and topographic effects in imagery because the landscape exhibits relatively high spatial variability of topographic properties and land cover variations due to complex climate, surface processes and tectonic interactions. Our work represents an extension of our previous modeling efforts over this area that accounts for atmosphere-topographic coupling that is not usually accounted for in ARC investigations. In this work, we focus on the development and evaluation of the E_t component and surface-irradiance partitioning to better understand if location- and wavelength-dependent topographic effects should be accounted for in ARC methods. More specifically, our objectives are as follows:

• Estimate and characterize the spatial variability of E_t over the landscape.

• Evaluate the significance of E_t parameterization schemes and various topographic effects.

• Based upon irradiance simulations of E_b, E_d and E_t, characterize the irradiance partitioning patterns over the landscape to identify the location- and wavelength-dependent nature of different topographic effects.

It is important to recognize that given the complexity of the problem and the inherent computational issues associated with estimating E_b, E_d and E_t, it is not possible to address all of the topographic effects that govern surface irradiance and reflectance (e.g. BRDF), such that our results represent a reasonable first-order approximation towards estimating the spectral-irradiance components, given prescribed atmospheric conditions. Furthermore, it is not the purpose of this study to conduct a comprehensive sensitivity analysis that formalizes the influence of numerous RTP parameters and their synergistic impact throughout the radiative-transfer cascade, or the cumulative accounting of energy across the short-wave region of the spectrum. Rather, we generate wavelength-dependent estimates, and compare them to other parameterization schemes to provide insights into the importance of topographic effects for estimating the surface irradiance and how irradiance partitioning could impact ARC investigations.

2. Background

The surface irradiance (E) is characterized as: (1) $$E(\lambda )=E_b(\lambda )+E_d(\lambda )+E_t(\lambda ).$$(1)

Topographic effects govern the magnitude and spatial variability of each component, such that there are spatial, spectral and temporal dependencies associated with sub-parameters and radiation-transfer processes. Unfortunately, there are a multitude of parameterization schemes for the various components and their sub-parameters, which vary in their degree of comprehensiveness, that are the result of specific assumptions, the relative importance of RTP, the degree to which a scheme accounts for various dependencies, and the computational approach that may be method-specific. Unfortunately, the use of different assumptions and parameterizations do not necessarily generate similar magnitudes or distributions, as clearly demonstrated by Bishop et al. (2022), given their sensitivity analysis results regarding E_b and E_d.

2.1. Direct irradiance

This component is modulated by Earth-Sun orbital dynamics, atmospheric processes that govern solar geometry and atmospheric attenuation, local-scale topographic effects, and meso-scale topographic effects that govern cast-shadows (Li et al. 2016; Bishop et al. 2019, 2022). Issues associated with accurately characterizing this component include accounting for time-dependent orbital dynamics, not using the small angle-approximation due to the spatial dependence of the solar geometry, accounting for atmospheric attenuation given the coupling of the atmosphere and topography, accurate estimation of local slope angle and slope azimuth conditions, and accounting for cast shadows, given the terrain relief structure in relation to orbital dynamics and solar geometry.

In empirical ARC investigations, the orbital dynamics, topographic effects governing solar geometry, atmosphere-topography coupling, and terrain-relief structure are not usually accounted for to reduce spectral variation in imagery. It is also important to note that the umbra and penumbra regions of cast shadows can increase the spatial variability in irradiance depending on the extent of shadowing as a function of time, such that the penumbera region of a cast shadow can exhibit surface irradiance from the direct-irradiance component (Cameron and Kumar 2018; Chu et al. 2021; Bishop et al. 2022). Common ARC methods effectively account for the local topographic effects, although other sub-parameters of direct irradiance are not usually utilized. Various researchers, however, have noted that atmosphere-topography coupling and cast shadows need to be accounted for (Li et al. 2016; Bishop et al. 2019).

2.2. Diffuse-skylight irradiance

The E_d component is known to account for Rayleigh and aerosol scattering, as well as secondary ground-atmospheric scattering that are related to atmospheric optical depth, ground-reflectance and sky-reflectance conditions (Proy et al. 1989; Gueymard 1995; Zhang et al. 2015). Atmospheric scattering is modulated by local topographic effects and meso-scale topographic shielding of scattered irradiance across the sky hemisphere, such that circumsolar and horizontal brightening govern the anisotropic-atmospheric scattering distribution over the sky hemisphere (Perez et al. 1986; Proy et al. 1989; Perez et al. 1990). Characterizing the anisotropic nature of this component is challenging given the spatio-temporal variability in atmospheric constituents (i.e. aerosol type and concentrations, and various degrees of clear skies given cloud conditions), as well as the computational complexity associated with coupling anisotropic scattering with anisotropic topography effects. Such interactions produces highly variable coupling conditions that can vary significantly over mountain landscapes (Bishop et al. 2022).

Investigators have attempted to characterize E_d assuming isotropic scattering conditions and have used a numerical scaling coefficient to represent the anisotropic nature of topographic shielding, which is commonly referred to as the skyview factor (e.g. Wen et al. 2009; Helbig et al. 2010; Wu et al. 2018). Other investigators use proxy parameters to account for meso-scale terrain anisotropy that are based upon local topographic properties (e.g. Temps and Coulson 1977; Perez et al. 1990). Such approaches to E_d estimation have been demonstrated to generate different magnitudes, variance structures and different spatial distribution patterns compared to relatively more complicated parameterization schemes that formally account for the anisotropic coupling of atmospheric scattering and topographic shielding (Bishop et al. 2022). Furthermore, many schemes do not account for anisotropic coupling of scattering due to local topographic effects.

Results from Bishop et al. (2022) suggest that perhaps the greatest challenge in estimating E_d will be to accurately characterize atmospheric conditions that govern the nature of anisotropic scattering. Atmospheric remote-sensing estimates regarding aerosols and cloud conditions can greatly facilitate numerical modeling efforts. Nevertheless, this irradiance component has not been adequately characterized or evaluated in common empirical ARC methods, where empirical data-based proxy parameters are used as a substitute, which cannot account for anisotropic couplings (Bishop et al. 2022). Consequently, the magnitude of this component is not known very well over complex mountain landscapes, and its overall surface energy contribution may be relatively high in cast-shadow regions and where terrain self-shadowing occurs due to steep slopes at lower altitudes. Not accounting for this parameter in ARC methods is thought to be the primary reason for ARC over-correction errors, although the E_t component may also dominate at high altitudes in cast shadow regions.

2.3. Adjacent-terrain irradiance

The surface energy contribution from this component is known to be important in surface energy-balance studies (Chu et al. 2021). Its importance in ARC, however, is not known with any certainty, and this component is not accounted for in most empirical ARC parameterization schemes (Gao et al. 2016). Accurate numerical estimates are difficult to generate due to the complexity of characterizing the anisotropic nature of terrain shielding, local topographic effects, atmospheric attenuation and the BRDF. Another complicating issue is characterizing the spatially-dependent influence of the adjacent terrain for a particular target location on the landscape, as contributions from the adjacent terrain are synergistically governed by the coupling of shielding, terrain self-shadowing, atmospheric attenuation and land-cover conditions. To our knowledge, the nature of these spatial-dependencies have not been adequately characterized or mapped, and research into this aspect of determining the impact of the adjacent terrain on this irradiance component is sorely needed. We speculate that the spatial properties of source areas is highly variable in a mountain landscape given spatial variations in topographic properties and land cover conditions, as mountain geodynamics can cause rapid environmental change given climate, surface processes and tectonic interactions (Shroder and Bishop 2000; Bishop et al. 2003).

Specific complications involve accurate characterization of the BRDF for adjacent-terrain source locations, and accounting for the collective topographic effects that govern E. The BRDF is governed by land cover structure and topographic effects, and the nature of anisotropic reflectance is not known with any degree of certainty in complex mountain environments (Zhang et al. 2015). Furthermore, the BRDF is governed by anisotropic adjacent-terrain reflectance, which partially governs the anisotropic nature of the BRDF. Bishop et al. (2022) refers to this issue as the ‘chicken or the egg’ problem, as E_t estimation is dependent upon knowing adjacent-terrain BRDFs, but a BRDF is governed by E that includes E_t. More research regarding numerical characterization of terrain modulated BRDFs in mountain environments is sorely needed.

Investigators have attempted to address these complexities in various ways by developing and evaluating various parameterization-schemes (e.g. ray tracing, accumulation and terrain-configuration factor approaches). Common limitations in generating estimates include not accounting for anisotropic surface reflectance, anisotropic atmospheric attenuation and anisotropic topographic effects. The Lambertian assumption is often necessary, as the variability of the adjacent-terrain BRDFs is not known, and we lack fundamental data to validate BRDF modeling efforts, although investigators have used satellite imagery to facilitate estimating numerical scaling coefficients for semi-empirical BRDF models. Atmospheric effects are usually considered to be insignificant, although extreme relief conditions and wavelenth-dependent atmospheric scattering can potentially cause attenuation at shorter wavelengths. Finally, investigators attempt to reduce the computational complexity by not accounting for anisotropic-terrain effects and generalize adjacent-terrain contributions using proxy parameters (i.e. assumed source-target distance or terrain-configuration factor) that can be calculated in a variety of ways. Consequently, we should not expect various parameterization schemes to generate accurate or precise results, and more research is required to formally evaluate topographic effects, and compare results from different approaches for estimating E_t (Chu et al. 2021).

2.4. Surface-irradiance partitioning

The reporting of surface-irradiance partitioning statistics from complex mountain locations are available for a handful of geographic locations at specific times (e.g. Dozier 1980; Gratton et al. 1994), however, accurately estimating the impact of E_d and E_t can be challenging given the spatio-temporal and spectral dependencies associated with topographic effects and land-cover conditions. It is generally thought that there is a systematic decrease in the magnitude of the irradiance components over the entire landscape, although this is most-likely not true in complex mountain environments where landscape complexity may govern component magnitudes and distribution patterns. Nevertheless, we do not have reasonable estimates of the spatial variability of partitioning magnitudes and patterns. Such information provides insights into the location-dependent topographic effects that need to be removed from satellite imagery.

It is known that E_b will not be the dominant irradiance component in the umbra portion of a cast-shadow region, where E_d or E_t could dominate depending upon atmospheric, terrain or land cover conditions. The magnitude and partitioning patterns in the penumbra region are very uncertain depending upon it’s length, such that E_d or E_t will dominate close to the umbra boundary, and E_b will systematically increase in dominance towards the furthest boundary of the cast shadow. Given that cast-shadow extent is dependent upon the coupled influence of geographic location and time (i.e. diurnal and annual), there could be significant spatial variations in irradiance partitioning patterns over a cast-shadow region. Such partitioning information over cast-shadow regions could be very helpful in explaining why most empirical ARC methods fail in these areas, and could help facilitate our understanding of how ARC methods should account for various types of topographic effects that dominate in various locations.

Similarly, the concept of spatial continuity and heterogeneity (i.e. spatial complexity) in surface biophysical properties can also modulate the magnitudes of E_d and E_t, and there may be spatial conditions that support relatively high or equivalent magnitudes of E_t compared to E_d. The spatial properties in surface reflectance, however, are closely related to the operational scale-dependencies of surface processes and topography, such that significant variations related to rock mass uplift, mass movements, fluvial erosion, glacierization and plant growth, have the potential to generate irradiance-partitioning patterns that are highly variable over a landscape. Much more research is required to assess the nature of irradiance component magnitudes and partitioning patterns over space and time.

Finally, the spatial complexity of terrain conditions related to topographic properties and terrain-orientation structure will also have an impact on magnitude and partition patterns, as topographic complexity is highly variable in mountain environments. Unfortunately, parameterization schemes and computational approaches may not account for a multitude of couplings involving interacting sub-parameters. Our research will formally address this important issue in our simulations to better estimate the distribution of partition patterns, and provide important information for addressing the limitations associated with empirical ARC methods.

3. Materials and methods

3.1. Data

We used datasets and simulation results from Bishop et al. (2022) to generate E_t estimates and irradiance partitioning statistics for the Nanga Parbat Massif. The irradiance components were simulated using an ASTER Global Digital Elevation Model Version 3 (ASTER GDEM; NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team 2009) to account for topographic effects. We accounted for atmospheric transmittance in all irradiance simulations based upon the U.S. standard atmosphere model and used mean exo-atmospheric irradiance values, atmospheric absorption coefficients for atmospheric constituents, and atmospheric scattering coefficients from Gueymard (1995) in simulating atmospheric conditions. See Bishop et al. (2022) for details on atmospheric parameterization schemes and constants used.

We simulated surface reflectance based on the spatial structure of generalized land cover conditions as described by Bishop et al. (2022). We used Landsat 8 OLI imagery captured in 2018 on August 4th (path 149, rows 35-36) and September 9th (path 150, rows 35-36) to do this. It is important to note that the use of the imagery was not to produce an accurate map of all land cover classes or biophysical variations. Nor did we make use of image solar geometry. We did use spectral indices to produce a first-order thematic map of the fundamental spatial distributions of water, snow, vegetation and non-vegetation, which served as a basis for establishing limits on what spectral signatures were used to characterize the biophysical variation within those land cover themes. In this way our simulations account for prescribed land cover and biophysical variation to estimate irradiance components, based upon prescribed surface reflectance conditions. We simulated three wavelengths for surface reflectance and all surface-irradiance components that correspond to the central spectral band wavelengths for OLI imagery in the green, red and near-infrared regions of the spectrum (${\lambda }_1=0.56141$ µm, ${\lambda }_2=0.65459$ µm, ${\lambda }_3=0.86467$ µm). The date and time of each simulation was 9-15-2022 at 10:00AM. Consequently, all irradiance images are in W m⁻²µm^–1.

3.2. Direct irradiance

We used E_b simulation results from Bishop et al. (2022). The parameterization scheme accounted for orbital dynamics, atmospheric-topography interactions, and local and meso-scale topographic effects, and was computed as: (2) $$E_b(\lambda )=E^0(\lambda )T^{\downarrow }(\lambda )\hspace{0.17em}\cos \hspace{0.17em}iS,$$(2) where E⁰ is the exoatmospheric irradiance, $T^{\downarrow }$ is the total downward atmospheric transmittance, i is the incidence angle that accounts for local topographic effects, and S represents the cast-shadow coefficient ($0.0-1.0$) that accounts for the meso-scale relief structure of the terrain. It should be noted that S encapsulates the umbra and penumbra regions of a cast shadow in our simulation.

3.3. Diffuse-skylight irradiance

We used E_d simulation results from Bishop et al. (2022). We computed E_d for every pixel and formally accounted for anisotropic atmosphere-topography coupling such that: (3) $$E_d(\lambda )={\int }_{{\varphi }_i=0}^{2\pi }{\int }_{{\theta }_i=0}^{\pi /2}{L}^{\downarrow }\left({\theta }_i,{\varphi }_i\right)\hspace{0.17em}\cos \hspace{0.17em}I{S}_t({\theta }_i,{\varphi }_i)\mathit{\text{dθdϕ}},$$(3) where $L^{\downarrow }$ is the downward radiance, θ_i and ${\varphi }_i$ are the zenith and azimuth angles of the incident energy from the sky hemisphere, I is the incidence angle of the direction defined by sky hemisphere and terrain geometry, and S_t is a binary coefficient which is 0.0 given terrain shielding (i.e. higher hemispherical zenith angles) and 1.0, given relatively lower hemispherical zenith angles.

It is important to note that we did not use the skyview factor in our computation of E_d, as our previous research revealed that its use does not generate similar magnitude and spatial variation results in E_d, compared to a full anisotropic computation of the hemispherical zenith angle (i.e. maximum relief angle that limits the directional hemisphere contribution) that characterizes topographic shielding. Our computation evaluated each direction to compute an integrated parameter that accounts for the anisotropic coupling of diffuse irradiance and topographic shielding. This also couples the local terrain first-order derivatives in relation to the hemisphere geometry to characterize the incidence angle.

3.4. Adjacent-terrain irradiance

We provide a first-order estimate of E_t by accounting for terrain blockage, atmospheric attenuation, local topographic conditions, and the inherent locational dependencies associated with adjacent source pixels that contribute reflected energy to a target pixel. We note that E_t is governed by the BRDFs from adjacent-terrain locations, however, we utilize the assumption of Lambertian surface reflectance, as the BRDF is also influenced by topographic effects (Wen et al. 2018; Bishop et al. 2019), although the nature and variability of BRDFs in complex mountainous terrain is difficult to model, as atmosphere, land cover structure and topographic effects have to be accounted for (Wen et al. 2018). Our parameterization scheme is as follows: (4) $$E_t(\lambda )={\int }_{{\varphi }_i=0}^{2\pi }{\int }_{{\theta }_i=0}^{\pi }L\left({\theta }_i,{\varphi }_i,\lambda \right)T_t^{\downarrow \uparrow }(\Delta z,{\theta }_v,\lambda )\hspace{0.17em}\cos \hspace{0.17em}{I}_t{S}_b\mathit{\text{dθdϕ}}$$(4) where ${\varphi }_i$ is the hemispherical incident azimuth angle, θ_i is the incident vertical hemispherical zenith angle, L represents incoming reflected energy that is direction dependent ($L=\rho (E_b+E_d)/\pi ,$ where ρ is the surface reflectance), $T_t^{\downarrow \uparrow }$ is the atmospheric transmittance given the optical depth of the atmosphere due to relief ($\Delta z$) and propagation zenith angle (θ_v) between two locations, I_t is the terrain incidence angle given the influence of the local terrain geometry in relation to the incident directional geometry, and S_b represents the potential of terrain blocking of the incident surface radiance between any two points of the surrounding terrain, as topographic eminences (i.e. crests and ridges) can extend above the trajectory line between two pixels that causes terrain blockage. For the terrain incidence angle, we specifically utilize the local terrain first-order derivatives at the target pixel in relation to the first-order derivatives of all contributing source pixels.

The computation of E_t is highly dependent on the number and spatial distribution of adjacent source pixels, as the topographic complexity caused by erosion and uplift dynamics at Nanga Parbat alters the relief structure and landscape dissection orientation of the terrain, such that the spatial distribution of adjacent pixels that contributes energy to a target pixel (source pixels) should be highly variable with respect to distance and direction. To our knowledge, the spatial-distributional nature of adjacent-terrain source locations have not been characterized or mapped over Nanga Parbat. It has been indicated that the computation of E_t should account for terrain conditions out to approximately 5 km from a pixel location (Proy et al. 1989), but this approach does not account for the terrain relief/orientation structure and the coupling of atmospheric and surface reflectance conditions. Consequently, we simulate E_t for a sample of locations and produce E_t maps for a relatively small area, given the significant amount of computation time associated with large-area estimation.

3.4.1. Scale-dependent analysis

Characterizing and mapping the scale-dependent nature of E_t was accomplished by designing our software so that we can identify adjacent-terrain source pixels that contribute energy to a target pixel over a particular computational radius, so that we can determine the functional relationship between the planimetric area of source pixels and the magnitude of E_t, and the variability of environmental factors governing E_t. We initially utilized a stratified-random sampling design to identify 500 locations over Nanga Parbat (Figure 1). We then computed E_t for each location over six computational radius limits (2.5, 5, 7.5, 10, 12.5 and 15 km) and extracted various parameters associated with the distance-dependent distribution of source pixels. Specifically, we examined the relationship between E_t and planimetric area (A_p), surface area (A_s), mean altitude ($\overline{z}$), maximum relief ($\Delta z_{\max }$), mean slope angle ($\overline{{\theta }_t}$) and average surface reflectance ($\overline{\rho }(\lambda )$).

Figure 1. Stratified-random sampling of locations over the Nanga Parbat Massif in Pakistan. Sample locations (yellow points) are depicted over a false-color composite of Landsat 8 imagery (R = NIR, G = red, B = green). the Landsat imagery is a mosaic of scenes that were selected for minimal cloud cover and temporal proximity: paths 149-150, rows 35-36, acquired on August 4th and September 12th, 2018.

We might expect that the magnitude of E_t will increase with larger computational radius distances, however, we might also expect that this trend may not be consistently uniform over complex mountain topography, as the scale-dependencies associated with terrain relief/orientation structure is highly variable over this region. Therefore, we map the adjacent-terrain source pixels using the largest radius by displaying the adjacent-terrain reflected radiance for source pixels that contribute to a target pixel. In this way, we can examine the spatial patterns and the nature of spatial continuity and heterogeneity associated with terrain and land cover complexity.

3.4.2. Parameter sensitivity analysis

We also compared our parameterization scheme to that of Proy et al. ($E_{t-\mathit{Proy}}$: 1989), which has been commonly referenced by investigators (e.g. Gratton et al. 1993; Richter 1997; Shepherd & Dymond 2003; Gao et al. 2016) such that: (5a) $$L_p(\lambda )=\rho (\lambda )\left[\frac{E_b(\lambda )+E_d(\lambda )}{\pi }\right],$$(5a) (5b) $$L(\lambda )=\frac{(L_p{A}_{\mathit{src}}\hspace{0.17em}\cos \hspace{0.17em}{\theta }_{\mathit{src}}{A}_{\mathit{tar}}\hspace{0.17em}\cos \hspace{0.17em}{\theta }_{\mathit{tar}})}{d^2},$$(5b) (5c) $$E_{t-\mathit{Proy}}={\int }_i^N{L}_{i}di,$$(5c) where L_p is the reflected radiance from a source pixel, L is the total irradiance from a source pixel, A_src and A_tar are the area of the source and target pixels, respectively, θ_src and θ_tar are the angles between the trajectory line from target to source pixel and the surface normal vector for the source and target pixel respectively, and d is the distance between the source and target pixel.

We also evaluated other commonly used parameterization schemes that are based upon a terrain-configuration factor (Chu et al. 2021). This approach is widely utilized (e.g. Kondratyev 1977; Dozier and Frew 1990; Shepherd & Dymond 2003; Wang et al. 2006), however, there are a multitude of parameterization schemes, and there does not appear to be consistent implementation related to the computation of the terrain-configuration factor (TCF) over the computational scale. Consequently, we evaluated three schemes for computing E_t based upon published work. Scheme one is computed as follows: (6) $$E_{t-\mathit{tcf}1}(\lambda )=F\left(\overline{\rho }(\lambda ){\overline{E}}^{\downarrow }(\lambda )\right),$$(6) where F is the terrain configuration factor that represents the fraction of adjacent source pixels found within the computational area that contributes energy to the target pixel, $\overline{\rho }$ is the average reflectance of the adjacent source pixels that contributes energy to the target pixel, and ${\overline{E}}^{\downarrow }$ is the average downward irradiance from the adjacent source pixels that contribute energy to the target pixel (i.e. $E^{\downarrow }=E_b+E_d$).

Scheme two is based upon the TCF of Kondratyev (1977) that is computed as: (7) $$E_{t-\mathit{tcf}2}(\lambda )=\left[\frac{1-\hspace{0.17em}\cos \hspace{0.17em}{\theta }_t}{2}\right]\overline{\rho }(\lambda ){\overline{E}}^{\downarrow }(\lambda ),$$(7) where θ_t is the slope angle at the target pixel.

Finally, we evaluate a TCF that is based upon the skyview factor (V_sky) that is used by many investigators (e.g. Helbig and Lowe 2012; Ma et al. 2016). It is computed as: (8) $$E_{t-\mathit{tcf}3}(\lambda )=\left(1.0-V_{\mathit{Sky}}\right)\overline{\rho }(\lambda ){\overline{E}}^{\downarrow }(\lambda ),$$(8) where V_sky was computed as: (9) $$V_{\mathit{sky}}=\sum _{\varphi =0.0}^{360}{\hspace{0.17em}\cos \hspace{0.17em}}^2{\theta }_{\max }(\varphi ,d)\frac{\Delta \varphi }{360}.$$(9)

It should be noted that we restricted the extent of the algorithm to find the maximum angles within the computational scale limit. Other parameterization schemes can also be used to account for the TCF using a multi-reflection term (Chu et al. 2021).

We then conducted a numerical sensitivity analysis similar to Bishop et al. (2022) to determine if a particular RTP, or a different parameterization scheme, significantly governs the magnitude and variance structure of E_t compared to our more comprehensive parameterization scheme (i.e. accounts for more topographic effects). This was done to identify important RTP that may need to be accounted for in ARC methods. We included or excluded specific RTP into our simulations that represent a series of control parameter scenarios, so that we can examine statistics and compare control parameter scenarios to the more comprehensive parameterization scheme. Specifically, we compute the minimum, maximum, mean and standard deviation for the control scenarios, and then compare results using the root-mean-squared error metric, the structural similarity index (SSI; Wang et al. 2005), Students t-test analysis for unequal variances, and an F-test to determine the significance of a parameter with respect to this component. See Table 1 for simulation and scenario descriptions.

Table 1. Simulation and control scenario (CS) descriptions for generating estimates of the adjacent-terrain irradiance.

CS	${\mathbf{P}}_{\mathbf{c}}$	Description
1	$E_{t-\mathit{Proy}}^i$	Proy parameterization scheme.
2	$E_{t-\mathit{tcf}1}^i$	Parameterization scheme based upon using a terrain configuration factor.
3	$E_{t-\mathit{tcf}2}^i$	Parameterization scheme based upon using a terrain configuration factor.
4	$E_{t-\mathit{tcf}3}^i$	Parameterization scheme based upon using a terrain configuration factor.
5	$S_b^x$	Terrain shielding.
6	TT^x	Total atmospheric transmittance.
7	${(\hspace{0.17em}\cos \hspace{0.17em}I)}^x$	Cosine of the terrain incidence angle.
8	${\mathbf{E}}_{\mathbf{t}}$	Base for comparisons using equation 4(4) $$E_t(\lambda )={\int }_{{\varphi }_i=0}^{2\pi }{\int }_{{\theta }_i=0}^{\pi }L\left({\theta }_i,{\varphi }_i,\lambda \right)T_t^{\downarrow \uparrow }(\Delta z,{\theta }_v,\lambda )\hspace{0.17em}\cos \hspace{0.17em}{I}_t{S}_b\mathit{\text{dθdϕ}}$$(4) .

Note: Control parameters (P_c) for a scenario were included (ⁱ) or excluded (^x) for comparison to a more fully comprehensive parameterization scheme.

3.5. Irradiance partitioning

Irradiance partitioning distribution patterns that reflect the dominance of the irradiance components are not generally known over complex mountain terrain. Such information needs to be systematically assessed with respect to space and time to determine which topographic effects govern irradiance and reflectance over the landscape. Furthermore, it is possible that the dominance of a particular irradiance component (i.e. highest percentage of energy) is a function of wavelength given atmosphere-topography-landcover coupling. Such information is essential to provide insights into effectiveness of various topographic correction methods.

To accomplish this, we used our E_b, E_d and E_t simulations over the central portion of the Nanga Parbat Massif to estimate the total surface irradiance (as a function of wavelength) and computed partitioning percentage images that depict the energy dominance of the components. The irradiance partitioning patterns (IPP) were classified for each location on the landscape such that: (10) $$\mathit{IPP}(x,y)=\{\begin{array}{cc}1\hfill & \text{if}\hspace{1em}{E}_b>E_d>E_t\hfill \\ 2\hfill & \text{if}\hspace{1em}{E}_b>E_t>E_d\hfill \\ 3\hfill & \text{if}\hspace{1em}{E}_d>E_t>E_b\hfill \\ 4\hfill & \text{if}\hspace{1em}{E}_d>E_b>E_t\hfill \\ 5\hfill & \text{if}\hspace{1em}{E}_t>E_b>E_d\hfill \\ 6\hfill & \text{if}\hspace{1em}{E}_t>E_d>E_b\hfill \end{array}$$(10)

The IPP locations were then statistically analyzed to characterize the surface-irradiance components. In this way, we generate information that can facilitate our understanding of the dominant types of topographic effects that need to be addressed by ARC methods at various terrain locations.

4 Results

4.1. Scale dependencies (E_t)

Computational scale-dependent E_t estimates for the sampled terrain locations are presented in Table 2. In general, larger computation scales are associated with a decrease in mean irradiance and standard deviation for all wavelengths. This trend may potentially be explained by atmospheric attenuation, as greater distances that account for more relief can decrease the surface radiance from adjacent-source pixels. As expected, when we compared smaller scales to our baseline scale (15 km), RMSE and SSI values systematically decreased and increased, respectively, for all wavelengths. At smaller scales, E_t estimates are not similar in magnitude, variance or correlation structure. This is highlighted by our inferential statistical tests where E_t magnitudes and variance at the 2.5 and 5.0 km scales were significantly different from the baseline scale.

Table 2. Adjacent-terrain irradiance estimates for sampled locations (Figure 1) based upon computational radius.

Radius [km]	Min [W m^−2µm^– 1]	Max [W m^−2µm^– 1]	$\overline{\mu }$ [W m^−2µm^– 1]	σ	RMSE [W m^−2µm^– 1]	SSI	t	F
$\lambda =0.56141$ [µm]
2.5	11.54	284.13	68.96	34.66	27.89	0.3857	−4.06*	0.53*
5.0	7.70	188.10	65.20	28.90	17.54	0.6467	−2.35*	0.77*
7.5	17.01	158.19	63.59	26.50	11.51	0.8197	−1.48	0.91
10.0	16.04	150.14	62.47	25.44	7.36	0.9196	−0.81	0.99
12.5	15.90	168.62	61.73	25.42	3.74	0.9788	−0.36	0.99
15.0	15.84	195.35	61.16	25.32	0.00	1.0000	0.00	1.00
$\lambda =0.65459$ [µm]
2.5	10.06	233.45	60.57	29.69	24.02	0.3837	−4.01*	0.55*
5.0	6.78	160.09	57.42	25.13	15.22	0.6481	−2.33*	0.77*
7.5	15.04	137.04	56.04	23.12	9.96	0.8222	−1.47	0.91
10.0	14.22	126.99	55.07	22.18	6.30	0.9224	−0.81	0.99
12.5	14.16	146.73	54.43	22.14	3.18	0.9798	−0.35	0.99
15.0	14.14	167.70	53.94	22.05	0.00	1.0000	0.00	1.00
$\lambda =0.86467$ [µm]
2.5	6.01	126.84	38.60	17.45	14.29	0.3900	−4.39*	0.59*
5.0	4.21	93.35	36.56	15.18	9.22	0.6515	−2.52*	0.78*
7.5	9.38	85.09	35.63	14.07	5.93	0.8312	−1.55	0.91
10.0	8.98	78.06	35.02	13.56	3.71	0.9283	−0.88	0.98
12.5	9.02	86.87	34.62	13.52	1.87	0.9815	−0.41	0.99
15.0	8.62	96.97	34.28	13.44	0.00	1.0000	0.00	1.00

Notes: All results were compared to a baseline representing results computed at the 15 km radius. Statistical parameter symbols include the minimum E_t (Min), maximum E_t (Max), mean E_t ($\overline{\mu }$), standard deviation (σ), root-mean-squared error (RMSE), structural similarity index (SSI), t-statistic (t), and the variance ratio (F). Inferential statistical values noted with a superscript symbol (*) represent significantly different magnitudes and variance. These t-statistics and F-ratios were significant at the $p\le 0.01$ level indicating that the magnitude and variance structure of E_t is significantly different due to a relatively small computational radius.

Nevertheless, E_t estimates at the 7.5, 10 and 12.5 km radius scales were not found to be significantly different in terms of magnitude and variance. We suspect that these results are specifically related to the relatively small sample size (n = 500) that does not capture the true magnitude and variance of the spatial variation of E_t with respect to the topographic spatial complexity at Nanga Parbat. A larger sample size would account for more topographic variation and increase the variability of E_t estimates. Consequently, we might expect that E_t estimates at the 7.5 and 10 km radius scales might then be significantly different. Regardless, our statistical results support the notion of E_t being highly dependent on computational scale (i.e. spatial extent used for performing calculations).

Additional computational evidence in support of this comes from computed planimetric and surface areas associated with adjacent-terrain source pixels that contribute reflected energy to sampled target pixels. Based upon our sampling, the adjacent-terrain source planimetric and surface areas ranged from $1-124$ and $1-150$ km², respectively, given a radius limit of 15 km. There was extensive variability in the source area from sampled target pixels that was governed by the terrain complexity related to relief and valley-orientation structure. We might expect a relationship between the magnitude of E_t and the adjacent-terrain source area size, but our analysis revealed no significant linear or non-linear relationships. We also examined the relationships between E_t and mean altitude, mean slope, mean reflectance and mean surface irradiance from the terrain source areas, although no significant linear or non-linear relationships were found. This is most likely the result of numerous factors that synergistically govern the magnitude of E_t including local topographic effects, atmospheric attenuation, variation in surface reflectance, and variations in surface irradiance given variations in the terrain source locations for a target pixel. Furthermore, the dominance of these factors is likely to be spatially variable across the landscape, thereby decreasing the causative role of source area size as being a primary governing factor. Other spatial concepts of scale may also be responsible for the scale-dependent nature of E_t.

We therefore mapped the spatial patterns of relatively large (60-200 km²) and intermediate (50-60 km²) scale-dependent terrain source areas surrounding 12 different target locations (Figures 2 and 3). Examination of both graphics clearly reveals numerous spatial properties associated with adjacent-source areas that characterize the nature of scale-dependence:

Figure 2. Relatively large terrain source areas for 6 target pixel locations over the Nanga Parbat Massif. The yellow dot in the center of each location is the target pixel. Reflected surface radiance values from the green portion of the spectrum ($\lambda =0.56141 \mu$m) that contribute to the target pixel defines the geographic distribution of terrain source areas.

Figure 3. Moderate sized ($50-60$ km²) terrain source areas for 6 target pixel locations over the Nanga Parbat Massif. The yellow dot in the center of each location is the target pixel. Reflected surface radiance values from the green portion of the spectrum ($\lambda =0.56141 \mu$m) that contribute to the target pixel defines the geographic distribution of terrain source areas.

1. Size variation. The contribution areas can vary significantly in their size (i.e. number of source pixels) from one target pixel to the next. This was described in our statistical results and can be visualized when comparing contribution areas in Figure 2 to Figure 3.

2. Shape variations. Spatially contiguous contribution areas can have a variety of shapes depending upon the nature of the topographic structure. Figures 2 and 3 both depict simple shapes and complex shapes, as defined by the perimeter to planimetric area ratio, and elongated shapes given valley orientation patterns.

3. Anisotropic variation. Adjacent-source areas can vary in their primary directional orientation from the target pixel. This can be characterized as an azimuthal range that depends upon the source area distribution. Target pixels could also have a relatively high isotropic source area distribution at the base of basin headwall regions, where slopes generally face towards the center of the valley floor.

4. Spatial continuity. Contribution areas can exhibit a range of spatial patterns, where contributions areas are not spatially contiguous. This is primarily governed by the terrain orientation structure in relation to the altitude of target pixels. Orientation structure governs terrain blockage thereby generating multiple spatially discontinuous source areas. The discontinuous patterns can be seen in both figures.

Collectively, our spatial maps reveal that the scale-dependent nature of E_t must be formally evaluated, and that any computation of parameters over the entire computation radius, for estimating E_t, may not adequately account for complex source area distribution patterns. Examination of Figures 2 and 3 also reveal that our 15 km computation radius does not account for other adjacent-source pixels outside our radius, as the maps clearly reveal source pixels at the very edge of our radius for some target pixels. Given the terrain complexity at Nanga Parbat, a larger radius is warranted to produce better estimates.

Finally, it should be noted that the adjacent-source area for a target pixel does not change as a function of wavelength (Figure 4). It is primarily governed by the topographic structure. The magnitude and spatial patterns of the reflected radiance distribution for the source area, however, does vary with wavelength, as the surface irradiance varies over the source areas given terrain variations, and reflectance values vary depending upon surface matter type and biophysical property variations.

Figure 4. Terrain source areas for two target pixel locations over the Nanga Parbat Massif. The yellow dot in the center of each location is the target pixel. A. Location 1. Reflected surface radiance values from the green portion of the spectrum ($\lambda =0.56141 \mu$m) that contribute to the target pixel. B. Location 1. Reflected surface radiance values from the red portion of the spectrum ($\lambda =0.65459 \mu$m) that contribute to the target pixel. C. Location 1. Reflected surface radiance values from the near-infrared portion of the spectrum ($\lambda =0.86467 \mu$m) that contribute to the target pixel. D. Location 2. Reflected surface radiance values from the green portion of the spectrum ($\lambda =0.56141 \mu$m) that contribute to the target pixel. E. Location 2. Reflected surface radiance values from the red portion of the spectrum ($\lambda =0.65459 \mu$m) that contribute to the target pixel. F. Location 2. Reflected surface radiance values from the near-infrared portion of the spectrum ($\lambda =0.86467 \mu$m) that contribute to the target pixel.

4.2. Parameter sensitivity analysis (E_t)

Control-scenario sensitivity analysis results are presented in Table 3. Parameterization scheme results (CS 1-4) demonstrate that commonly utilized parameterization schemes that are based upon Proy et al. (1989) and the terrain-configuration factor, do not produce E_t estimates that are comparable to our more comprehensive parameterization scheme. Root-mean-squared error values were relatively high and SSI values are relatively low indicating that there are significant differences in the magnitude, variance, and correlation structure of the estimates. Inferential statistical results also support this conclusion.

Table 3. Control scenario (CS) statistical results for adjacent-terrain irradiance.

CS	Pc	Min	Max	$\overline{\mu }$	σ	RMSE	SSI	t	F
$\lambda =0.56141$ [µm]
CS1	$E_{t-\mathit{Proy}}^i$	0.64	1009.40	32.68	54.48	66.73	0.0001	10.60	0.22
CS2	$E_{t-\mathit{tcf}1}^i$	0.21	71.91	18.14	13.63	49.35	0.0306	33.46	3.45
CS3	$E_{t-\mathit{tcf}2}^i$	0.01	247.68	34.41	25.27	30.14	0.5267	16.72	1.00*
CS4	$E_{t-\mathit{tcf}3}^i$	2.27	265.57	80.47	41.71	41.49	0.1964	−8.85	0.37
CS5	$S_b^x$	10.12	110.16	47.68	16.37	21.00	0.5294	10.00	2.39
CS6	TT^x	17.07	205.07	63.83	26.54	3.14	0.9953	−1.62*	0.91*
CS7	${(\hspace{0.17em}\cos \hspace{0.17em}I)}^x$	67.36	288.27	137.84	35.40	82.12	0.1717	−39.40	0.51
CS8	${\mathbf{E}}_{\mathbf{t}}$	15.84	195.35	61.16	25.32	0.00	1.0000	0.00*	1.00*
$\lambda =0.65459$ [µm]
CS1	$E_{t-\mathit{proy}}^i$	0.56	912.37	28.68	48.94	59.70	0.0003	10.52	0.20
CS2	$E_{t-\mathit{tcf}1}^i$	0.19	64.02	15.78	11.92	43.59	0.0301	34.05	3.42
CS3	$E_{t-\mathit{tcf}2}^i$	0.01	209.26	29.89	21.71	26.87	0.5190	17.38	1.03*
CS4	$E_{t-\mathit{tcf}3}^i$	2.00	237.69	69.66	35.41	35.43	0.1844	−8.43	0.39
CS5	$S_b^x$	8.75	96.22	42.43	14.41	18.15	0.5319	9.77	2.34
CS6	TT^x	14.66	172.73	55.52	22.75	1.89	0.9976	−1.12*	0.94*
CS7	${(\hspace{0.17em}\cos \hspace{0.17em}I)}^x$	57.54	242.69	121.52	30.35	72.20	0.1690	−40.29	0.53
CS8	${\mathbf{E}}_{\mathbf{t}}$	14.14	167.70	53.94	22.05	0.00	1.0000	0.00*	1.00*
$\lambda =0.86467$ [µm]
CS1	$E_{t-\mathit{proy}}^i$	0.34	566.80	18.17	30.73	37.79	0.0016	10.74	0.19
CS2	$E_{t-\mathit{tcf}1}^i$	0.12	41.75	9.83	7.46	27.63	0.0305	35.56	3.25
CS3	$E_{t-\mathit{tcf}2}^i$	0.01	119.09	18.59	13.15	17.28	0.5092	18.65	1.05*
CS4	$E_{t-\mathit{tcf}3}^i$	1.29	151.88	43.14	21.10	21.55	0.1544	−7.92	0.41
CS5	$S_b^x$	5.46	57.66	27.77	9.06	10.78	0.5397	8.97	2.20
CS6	TT^x	8.70	97.80	34.54	13.56	0.31	0.9998	−0.30*	0.98*
CS7	${(\hspace{0.17em}\cos \hspace{0.17em}I)}^x$	39.48	136.26	77.00	17.37	45.30	0.1597	−43.49	0.60
CS8	${\mathbf{E}}_{\mathbf{t}}$	8.62	96.97	34.28	13.44	0.00	1.0000	0.00*	1.00*

Control parameters (P_c) for a scenario were included (ⁱ) or excluded (^x) for comparison to a more fully comprehensive parameterization scheme (${\mathbf{E}}_{\mathbf{t}}$). statistical parameter symbols include the minimum E_t (Min), maximum E_t (Max), mean E_t ($\overline{\mu }$), standard deviation (σ), root-mean-squared error (RMSE), structural similarity index (SSI), t-statistic (t), and the variance ratio (F). control parameters symbols are defined in the methodology and symbol notation table in the appendix (Table A1). Inferential statistical values noted with a superscript symbol (*) represent no significant difference in the mean or variance. All other t-statistics and F-ratios were significant at the $p\le 0.00001$ level for all regions of the spectrum, indicating that the magnitude and variance structure of E_t is significantly different when parameters characterizing topographic effects are not accounted for, or when a different parameterization scheme is used.

We did not expect the results to be so different, as the implementation of these parameterization schemes were all based upon our analysis, such that the computation accounted for source pixel distributions. In the case of the Proy et al. (1989) algorithm, atmospheric attenuation was not accounted for, local topographic effects are characterized differently, and a target pixel to source-pixel distance parameter governs the magnitude of the estimate. We suspect that the coupling of local topographic effects and atmospheric attenuation generates significant spatial variation that cannot be accounted for in this algorithm, and the distance parameter does not accurately characterize a radiation-transfer process.

Similarly, the terrain-configuration factor algorithms utilize parameters that do not accurately characterize surface radiance, irradiance and terrain variation within adjacent-source areas. Local topographic conditions and topographic-atmospheric couplings are not accounted for, and these parameters can be significant given the extreme relief and relatively high anisotropic terrain conditions. The basic statistical approach does not formalize radiation-transfer processes, and effectively relies on a terrain-scaling factor to obtain estimates within a range. This, however, cannot account for the spatial variation in irradiance, as demonstrated in our results.

Sensitivity analysis of the parameters in our E_t parameterization scheme are presented in Table 3 (CS 5-7). Local and meso-scale topographic effects were the most significant factors. Significantly different results were obtained when not accounting for terrain blockage and the terrain incident angle for all wavelengths. The most significant factor was local topographic effects, as the $\hspace{0.17em}\cos \hspace{0.17em}{I}_t$ parameter exhibited the lowest SSI value. Given the extreme relief at Nanga Parbat, we expected the total transmittance parameter to significantly govern E_t variation, especially in the green region of the spectrum. We found, however, no significant difference in E_t magnitude or variation when we did not account for atmospheric attenuation. These results, however, are based upon the clear atmosphere assumption that was used in our simulations. Increasing the atmospheric aerosol content and accounting for a larger area in estimating E_t for a target pixel could account for more spatial variation in this irradiance component.

4.3. Irradiance components

Simulations of the direct irradiance (Figure 5; A,B,C) clearly depict significant spatial variation caused by cast shadowing and local topographic effects. The overall magnitude is greatest in the green portion of the spectrum and systematically decreases with increasing wavelength given the exoatmospheric irradiance from the sun. The magnitude of the spatial variation in E_b is greatest at shorter wavelengths given the coupled influence of atmosphere-topography coupling that is dominant at shorter wavelengths. Nevertheless, it is important to note that that local topographic and cast shadowing effects are present irrespective of wavelength.

Figure 5. Surface irradiance components over the Nanga Parbat Massif. A. Direct irradiance (green; 0.56141 µm). B. Direct irradiance (RED; 0.65459 µm). C. Direct irradiance (NIR; 0.86467 µm). D. Diffuse-skylight irradiance (green; 0.56141 µm). E. Diffuse-skylight irradiance (RED; 0.65459 µm). F. Diffuse-skylight irradiance (NIR; 0.86467 µm). G. Adjacent-terrain irradiance (green; 0.56141 µm). H. Adjacent-terrain irradiance (RED; 0.65459 µm). I. Adjacent-terrain irradiance (NIR; 0.86467 µm).

Simulations of diffuse-skylight irradiance (Figure 5; D,E,F) clearly depict spatial variation caused by anisotropic diffuse irradiance, local topography, and coupled atmosphere-topography coupling, such that the magnitude of E_d is generally lower on slopes that face away from the sun. The relief also accounts for altitudinal variation in irradiance, with higher irradiance at lower altitudes, although this variation is difficult to visually see given more dominant variations in local topographic effects and hemispherical topographic shielding. Similarly, the general magnitude of irradiance decreases with increasing wavelength, and topographic effects and patterns persist irrespective of wavelength.

Simulations of adjacent-terrain irradiance (Figure 5; G,H,I) clearly depict significant spatial variation caused by land cover, local topography, terrain blockage and atmosphere-topography coupling. The terrain orientation fabric strongly governs the distribution of irradiance, and more extensive hillslopes that are adjacent to target pixels that face towards them govern higher E_t values. Similarly, the general magnitude of irradiance decreases with increasing wavelength, and the topographic effects and patterns persist irrespective of wavelength.

It is important to note that the spatial variance structure and patterns of irradiance are different for each of the irradiance components. Cast shadows and local topography in relation to solar geometry most strongly governs E_b, while anisotropic atmosphere-topography coupling dominates E_d variation. Terrain relief and orientation fabric strongly controls E_t variation and produces different irradiance patterns. Consequently, knowledge of irradiance partitioning patterns represents an important aspect of ARC, as different topographic effects dominant in different locations on the landscape. This strongly suggests that ARC methods and modeling efforts may need to account for all three, or a subset of irradiance components for some geographic regions.

4.4. Irradiance partitioning

Irradiance partitioning images are presented in Figure 6. In general, the surface irradiance in the green region of the spectrum (panels A,D,G) is dominated by E_b, sequentially followed by E_d and E_t. This same general pattern is not found for the red (panels B,E,H) and near-infrared (panels C,F,I) regions of the spectrum, where E_t percentages are relatively large compared to E_d. The first partitioning partitioning pattern (IPP 1) is generally thought to be consistent over the landscape and wavelength regions, although our simulations demonstrate that this may not be the case, especially at longer wavelengths and in areas of cast shadows. Within cast shadow regions, the surface irradiance is thought to be dominated by E_d, although our results clearly reveal that E_t contributes more energy, with its spatial dominance increasing with an increase in wavelength.

Figure 6. Spectral irradiance partitioning patterns of overall surface irradiance. A. Direct irradiance contribution (green; 0.56141 µm). B. Direct irradiance contribution (RED; 0.65459 µm). C. Direct irradiance contribution (NIR; 0.86467 µm). D. Diffuse-skylight irradiance contribution (green; 0.56141 µm). E. Diffuse-skylight irradiance contribution (RED; 0.65459 µm). F. Diffuse-skylight irradiance contribution (NIR; 0.86467 µm). G. Adjacent-terrain irradiance contribution (green; 0.56141 µm). H. Adjacent-terrain irradiance contribution (RED; 0.65459 µm). I. Adjacent-terrain irradiance contribution (NIR; 0.86467 µm).

Wavelength dependence is also noted for E_b and E_d, where the general dominance of E_b increases and E_d decreases in dominance with an increase in wavelength. This is primarily caused by the decreasing influence of atmospheric scattering which is highly wavelength dependent. For E_t the wavelength dependence is thought to be related to decreased scattering and an increase in surface reflectance associated with minerals/rocks and vegetation.

The irradiance partitioning patterns for the green, red and NIR regions of the spectrum are displayed in Figure 7, while the statistics for surface irradiance components are presented in Table 4. Cast shadow regions (IPP 6) found at relatively high altitudes exhibit surface irradiance dominated by E_t and E_d components. In general, the surface irradiance decreases with increasing wavelength, and the areal percentage of this pattern did not change significantly as a function of wavelength.

Figure 7. Spectral irradiance partitioning patterns over the Nanga Parbat landscape. A. Pattern distribution in the green spectral region (0.56141 µm). B. Pattern distribution in the red spectral region (0.65459 µm). C. Pattern distribution in the NIR spectral region (0.86467 µm).

Table 4. Surface irradiance statistical results for spectral irradiance partitioning patterns (IPP) over the Nanga Parbat landscape.

IPP	%	Eb				Ed				Et				E
		Min	Max	$\overline{\mu }$	σ	Min	Max	$\overline{\mu }$	σ	Min	Max	$\overline{\mu }$	σ	Min	Max	$\overline{\mu }$	Σ
$\lambda =0.56141$ [µm]
1	84.29	78.87	1188.26	810.00	256.55	63.65	189.26	140.40	22.00	3.85	150.52	62.78	22.92	217.18	1476.79	1013.18	269.94
2	9.90	84.44	1055.66	342.50	137.44	67.91	150.69	100.08	8.90	78.23	248.90	126.70	18.71	247.71	1287.15	569.28	146.73
3	0.56	0.00	99.58	8.63	22.28	41.11	118.63	96.37	7.10	3.05	111.02	83.37	12.83	44.65	303.11	188.37	28.47
4	0.02	61.16	103.90	88.24	8.48	78.07	118.01	95.06	6.85	29.87	100.77	79.12	12.08	203.06	305.90	262.42	22.53
5	0.44	74.57	239.24	118.03	21.60	68.82	103.62	87.52	6.12	83.59	239.56	142.73	24.45	238.13	559.93	348.28	39.89
6	4.79	0.00	99.97	6.46	19.66	40.97	112.82	85.87	8.10	67.20	274.52	131.41	26.63	132.90	373.74	223.74	28.69
$\lambda =0.65459$ [µm]
1	64.94	75.03	1038.50	785.16	172.65	42.53	102.23	80.51	10.14	3.41	97.32	46.97	15.69	166.73	1225.51	912.63	182.82
2	29.27	53.47	1035.66	398.30	176.38	34.66	97.68	60.55	8.17	35.03	205.73	89.18	18.56	151.59	1225.12	548.03	181.18
3	0.04	0.00	55.67	0.48	4.51	22.67	64.22	55.49	4.67	1.77	61.18	45.54	8.74	24.58	177.12	101.51	12.92
4	0.00	54.89	58.95	57.55	2.30	57.28	62.56	60.28	2.71	25.90	57.90	43.89	16.37	145.84	179.27	161.71	16.78
5	0.77	38.76	169.41	85.88	22.70	37.38	60.55	48.98	3.49	55.69	188.55	112.87	21.34	149.54	407.02	247.73	37.72
6	4.97	0.00	56.00	1.81	7.69	22.59	64.61	47.23	4.66	43.99	229.22	111.44	23.86	84.56	251.81	160.48	22.06
$\lambda =0.86467$ [µm]
1	19.62	36.86	648.95	522.71	72.35	14.50	28.61	23.06	2.27	2.11	28.42	17.65	4.36	69.31	703.90	563.42	76.59
2	74.64	33.43	649.14	387.41	159.97	9.62	28.62	20.07	3.76	14.29	112.53	41.33	13.09	78.55	730.25	448.81	156.54
3	0.00	0.00	0.00	0.00	0.03	6.38	16.66	13.58	3.23	0.49	16.49	9.20	6.56	6.87	33.16	22.78	9.61
4	0.00	–	–	–	–	–	–	–	–	–	–	–	–	–	–	–	–
5	0.94	10.78	93.61	45.41	16.25	10.32	17.15	13.77	0.99	29.35	103.89	66.42	10.16	57.69	202.94	125.59	22.06
6	4.81	0.00	16.00	0.20	1.38	6.37	18.10	13.13	1.32	15.79	125.99	67.73	12.86	27.73	132.35	81.07	12.16

Statistical parameter symbols include the areal percentage of a sequence (%), minimum irradiance (Min), maximum irradiance (Max), mean irradiance ($\overline{\mu }$) and standard deviation (σ).

Our results demonstrate that IPP 2 ($E_b>E_t>E_d$) was the second most dominate irradiance pattern over the landscape, although the dominance and spatial variability of this pattern is highly variable. These areas receive a considerable amount of surface irradiance, and we note that the spatial percentage of this pattern significantly increases with an increase in wavelength, such that in the NIR, 74.64% of the region was dominated by this irradiance pattern. This result was unexpected and demonstrates the significance of the E_t component in accounting for surface irradiance variation. All the other irradiance patterns did not have a significant spatial presence on the landscape, given the date and time associated with simulations.

Collectively, these results have significant implications for ARC of satellite imagery in topographically complex mountain environments. They clearly demonstrate that different topographic effects dominate over different geographic locations and that the contributions of various irradiance components exhibit wavelength dependence. It is also imperative to note that irradiance variations are spatio-temporal scale dependent with respect to coupled atmosphere-topography-landcover variations that govern the BRDF and the anisotropic surface reflectance variations that are recorded by platform sensors.

5. Discussion

5.1. Scale dependencies

In the remote sensing community, a variety of spatial concepts have not been accounted for to address the irradiance variations needed for adequate ARC. Instead, traditional empirical ARC methods rely upon using a numerical scaling-factor coefficient in topographic-correction algorithms in an attempt to reduce topographic effects in imagery (Bishop et al. 2022). It is important to note, that the spatial properties that regulate irradiance variations are unique for all of the irradiance components that are a function of time and location. Bishop et al. (2022) demonstrate this for the E_b and E_d components, especially for cast shadows and the spatial coupling of anisotropic diffuse irradiance and topographic shielding, respectively.

Our results also clearly depict the nature of scale dependencies associated with the E_t component at the pixel level. This component is spatially variable depending upon land-cover distribution, terrain blockage from adjacent terrain to a target pixel, and atmospheric and local topographic effects. We have demonstrated that the landscape spatial influence on irradiance can occur over larger distances than what has been described in the literature, and that the nature of scale dependence is spatially complex, in terms of direction, distance, areal size and spatial continuity. Such complex spatial characteristics are governed by the relief structure, terrain-orientation fabric, and distribution patterns of landcover. At Nanga Parbat, we found that we needed to evaluate the terrain past 15 km from a target pixel for some locations, in order to accurately account for source pixels. This would increase the sampling of vegetation and snow source pixels, which could increase the magnitude of E_t values at target pixels in the visible spectrum (i.e. snow) and NIR region (i.e. vegetation).

Our results, when compared with Bishop et al. (2022) demonstrate that the unique scale dependencies of the irradiance components may not be spatially coincident, which governs the spatial variability in surface irradiance and dictates irradiance patterns over the landscape. Ultimately, the nature of irradiance partitioning is governed by the degree of environmental spatial and temporal landscape complexity. Consequently, effective ARC requires knowledge of the dominant scale-dependent topographic effects that govern surface irradiance variations at a particular location.

5.2. Et parameterization schemes

Our results clearly demonstrate that the use of different parameterization schemes will potentially have a significant impact on the magnitude and spatial variability of adjacent-terrain irradiance. Similar results were obtained from (Bishop et al. 2022) for the E_b and E_d components. Given the importance of the E_t component in energy-balance studies and our partitioning results, integrative coupling appears to be required in order to account for scale-dependencies. This parameter is rarely included in empirical ARC parameterization schemes. When evaluated, researchers commonly use the Proy et al. (1989) approach or schemes that are based upon a terrain configuration factor, which effectively represents numerical scaling. In our work, we noted that Proy’s scheme is heavily dependent upon the distance parameter with a 2D characterization producing lower irradiance values compared to a 3D characterization, all other factors remaining constant. Furthermore, we noted that irradiance values are highly dependent on the computational scale, which governs the distance parameter. Consequently, small scales of analysis produce more irradiance while larger scales of analysis produce less irradiance. This represents a disconnect with radiation-transfer processes such as atmospheric attenuation. Given the complexity of location-dependent adjacent-source areas and their inherent spatial variability, a fixed computation scale warrants significant error and uncertainty regarding E_t estimates.

Results from terrain-configuration schemes also produced significantly different results from our more comprehensive scheme. This was anticipated, as we noted that the computation of the terrain scaling coefficients are highly ambiguous with respect to characterizing topographic effects, and that generalized parameters cannot adequately account for variations in adjacent-terrain reflectance and irradiance. It should be noted that the first problem is selecting an adequate computational scale, which is effectively unknown for any particular location, given our analysis results or variations in topographic complexity. One must then decide if summary statistics from within the computation scale (i.e. window) should be used or be restricted to just those source areas. Finally, the percentage of the source pixels within the window, which is not related to any processes, systematically changes as a function of computation scale. We note that data outside of the source areas should not be used, and that the terrain numerical scaling factor has no physical meaning.

Similarly, computing the terrain-configuration factor based upon the skyview coefficient does not have any physical meaning, as one could compute this parameter over the computational scale, over the adjacent sources areas, or for the basin that the target pixel is found within. The notion here is that the adjacent-terrain influence numerically scales with skyview, although our results demonstrate that not all directions of the landscape have contributing source areas, that not all altitudes are represented in sources areas, and that source areas can be spatial discontinuous. Consequently, the relief-based parameter cannot adequately numerically scale the results in a meaningful way. We suspect that this approach is a simplistic way to numerically scale or calibrate irradiance in less complex environments given field reference data.

Finally, we note that excluding terrain blockage or local topographic effects generates results that are significantly different compared to our scheme that incorporates these parameters. Given the 15 km radius distance of our analysis, we found that excluding the atmospheric transmittance parameter does not produce significantly different results. We were surprised by this finding, as the extreme relief and altitudinal land cover gradients at Nanga Parbat at larger scales would incorporate more highly reflective features (i.e. eroded hillslopes, vegetation, glaciers, snow) and increase atmospheric optical depth. Analysis at a larger scale is required to determine the degree to which more sources areas contribute to the irradiance component for those locations exhibiting relatively high relief. Furthermore, we note a limitation to our analysis, in that we do not account for the surface BRDF in computing E_d and E_t estimates. We also do not iteratively or recursively account for E_t in estimating E_t. More research on accounting for the anisotropic nature of surface reflectance is sorely needed, as the BRDF is governed by topographic effects and the land-cover structure.

Accurately estimating the surface irradiance is therefore dependent upon adequate parameterization schemes that can uniquely characterize the spatial properties associated with topographic effects for all the irradiance components. Adequate characterization accounts for temporal and spatial scale dependencies for individual sub-parameters, which collectively governs the dominance of different terrain effects on the landscape.

5.3. ARC

Our irradiance partitioning results provide new insights into the effective utility of ARC methods in complex mountain environments. In our simulations, we found that irradiance partitioning patterns vary over the landscape. The dominant irradiance component was found to significantly change over the different wavelengths that we examined, and exemplifies the differential influence of scale-dependent topographic effects. Furthermore, given contributions of all irradiance components to surface irradiance over our region, it is clear that numerous topographic effects are generating variations in surface irradiance that ultimately govern the nature of the BRDF along with intimate and areal matter-property distributions and spatial structure. Consequently, ARC methods need to account for the integrative coupling of topographic effects at sub-pixel, local and regional scales. More research relating the irradiance components to the nature of the BRDF is sorely needed to account for the complex nature of the BRDF.

Ultimately, the degree of environmental complexity regulates the nature of the spatial properties that govern topographic effects and irradiance partitioning. The concept of complexity can be viewed from a structural perspective that accounts for spatial variation in topographic and land-cover properties and patterns. Complexity can also be viewed from a systems perspective that accounts for the spatio-temporal variability in governing factors that modulate radiation-transfer processes that ultimately control spectral variation in satellite imagery. Given the dynamic nature of the problem with respect to topographic effects, it seems clear that ARC parameterization schemes need to account for parameters that can represent different irradiance components and the spatial variability of irradiance depending upon the degree of environmental/topographic complexity.

With regard to complexity, many may assume that our results may not be applicable to other less complex mountain environments, as Nanga Parbat exhibits extreme relief and therefore it is assumed to exhibit relatively high topographic complexity. The topic of topographic complexity has been recognized to be important in numerous Earth science and remote sensing investigations. Researchers have attempted to characterize this concept based upon on local statistical characterization, spatial/structural characterizations and more rigorous system approaches. Topographic complexity is a multi-faceted concept and we argue that the topographic complexity at Nanga Parbat is relatively moderate (i.e. based upon local variance of topographic properties). The nature of the geodynamics over time (i.e. rapid uplift and erosion) has generated relatively low to moderate spatial frequencies of topographic property variation (i.e. topographic properties do not significantly vary over short distances). Many other mountain systems actually exhibit greater topographic complexity (i.e. higher spatial frequency variations of topographic properties), although the relief is much less. Temporal frequency variations in irradiance and partitioning should also be taken into consideration, as space-time effects are coupled.

Our simulations effectively represent a systems approach to evaluating topographic complexity, as we account for a multitude of topographic effects and the spatial variation in irradiance components and partitioning. Although we do not evaluate the temporal dependence of RTP over this area, there is considerable spatial variation in irradiance and partitioning patterns. In a more complex mountain environment exhibiting higher spatial frequency variations in topographic and land cover conditions, the irradiance variations and partitioning sequences could vary even more. We are aware of many other mountain systems that exhibit such high frequency variations in topographic properties. Consequently, we argue that our results have potential applicability to many other mountain systems where ARC is also required.

Regarding the influence of atmospheric constituents on E_t, while it is most likely true that higher aerosol content could decrease terrain irradiance in the shorter wavelengths, our results demonstrate the importance of this irradiance component at longer wavelengths. We would argue that land cover and topographic complexity are more dominant factors than aerosol content/transmittance, especially given higher topographic complexity, as shorter distances between source and target pixels and less relief results in less aerosol mass.

Unfortunately, empirical ARC methods do not account for several E_b sub-parameters, E_d, E_t and BRDF sources of variation. Bishop et al. (2022) discussed the problems associated with ARC parameterization schemes and noted the strong possibility that such scaling methods compress information and may potentially cause over-correction and/or the destruction of information. Furthermore, the degree of over-correction or information compression due to numerical scaling has not been spatially evaluated by the research community. ARC methods should be able to account for the inherent spatial variations associated with topographic effects and couplings in order to facilitate biophysical assessment and mapping in mountain environments.

6. Conclusions

Biophysical remote sensing and thematic mapping of mountain environments is extremely difficult due to atmosphere-topography-landcover couplings that govern the irradiant fluxes and the BRDF. Researchers in the Earth sciences have noted that all of the irradiance components are important, although they have not been adequately accounted for in terms of operational scale dependencies. Most empirical anisotropic-reflectance correction (ARC) approaches are based upon a numerical scaling coefficient that does not adequately account for various topographic effects. Consequently, we evaluated the spatial variation of surface spectral irradiance components with respect to parameterization schemes, computational scales, parameter sensitivity analysis and irradiance-partitioning patterns to enable new insights into the complex nature of topographic modulation of the RTC.

Our results clearly depict spatially-dependent adjacent-source areas that govern E_t. We note that adjacent source areas can be found at larger distances than what has been described in the literature, and that the spatial nature of topographic effects is complex, in terms of direction, distance, areal size and spatial continuity. Such complex characteristics were found to be governed by the relief structure, terrain-orientation fabric, and distribution patterns of landcover. Furthermore, compared to our previous work, we found that the operational scale dependencies associated with various irradiance components may not be spatially coincident, which governs the spatial variability in surface irradiance and dictates the dominance and distribution of irradiance patterns. Consequently, effective ARC requires knowledge of the dominant topographic effects that govern reflectance at a particular location.

Many irradiance parameterization schemes do not adequately account for multi-scale topographic effects. Given parameter sensitivity analysis, we found this to be the case for all the irradiance components for the Nanga Parbat region. Most schemes do not account for the fundamental anisotropic couplings between the atmosphere and topography that govern different patterns of E_d and E_t. We found that Proy’s scheme is heavily dependent upon the distance parameter, therefore dictating the selection of an adequate computational radius, although results were found to be significantly different compared to a more comprehensive scheme. Results from terrain-configuration schemes also produced significantly different results, as generalized parameters cannot adequately account for spatial variations in adjacent-terrain reflectance and irradiance.

Finally, it is important to note that dominance-based irradiance patterns can vary over the landscape. The irradiance patterns were also found to significantly change over the different wavelengths that we examined, and exemplifies the differential influence of atmosphere-topography-landcover coupling. Furthermore, given variations in irradiance partitioning patterns and topographic anisotropy, it is clear that numerous topographic effects are potentially governing the anisotropic nature of the BRDF. Consequently, ARC methods should attempt to account for multi-scale topographic effects that incorporate sub-pixel, local and regional anisotropic couplings.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

We used Landsat 8 OLI imagery captured on 2018 on August 4th (path 149, rows 35-36) and September 9th (path 150, rows 35-36) and the ASTER GDEM version 3. These can be accessed online (e.g. through Earth Data; https://www.earthdata.nasa.gov/). This research and C++ implementations of equations presented in this paper and referenced literature represents ongoing work and serves as the basis for future research. Simulation results are available upon reasonable request.

Appendix A.

 Mathematical symbol notation

Table A1. Mathematical symbols used in radiation-transfer cascade calculations.

Symbol	Units	Description
$\Delta z$	[m]	Relief.
λ	[µm]	Wavelength of light.
ϕ	[radians]	Azimuth angle.
ϕi	[radians]	Azimuth angle of hemispherical incident energy.
ϕs	[radians]	Solar azimuth angle.
ρ	[kg m^−3]	Surface reflectance.
$\overline{\rho }$	dimensionless	Average reflectance.
θmax	[radians]	Maximum relief angle to the horizon.
θi	[radians]	Zenith angle of hemispherical incident energy.
θs	[radians]	Solar zenith angle.
θmax	[radians]	Maximum horizon angle.
θt	[radians]	Slope angle.
$\overline{{\theta }_t}$	[radians]	Average slope angle.
θsrc	[radians]	Slope angle of adjacent source pixel that reflects light onto the target pixel.
θtar	[radians]	Slope angle of the target pixel.
Ap	[m^2]	Planimetric area.
As	[m^2]	Surface area.
Asrc	[m^2]	Area of the adjacent source pixel that reflects light onto the target pixel.
Atar	[m^2]	Area of the target pixel.
d	[km]	Distance between two pixels.
E	[W m^−2µm^–1]	Surface irradiance.
E^0	[W m^−2µm^–1]	Exoatmospheric irradiance.
Eb	[W m^−2µm^–1]	Direct-beam irradiance from the Sun.
Ed	[W m^−2µm^–1]	Diffuse-skylight irradiance.
Et	[W m^−2µm^–1]	Adjacent-terrain irradiance.
${\overline{E}}^{\downarrow }$	[W m^−2µm^–1]	Average downward irradiance from adjacent terrain that contributes energy to the target location (i.e. $E^{\downarrow }=E_b+E_d$).
Et– Proy	[W m^−2µm^–1]	Adjacent-terrain irradiance, after Proy et al. (1989).
$E_{t-\mathit{tcf}1}$	[W m^−2µm^–1]	Adjacent-terrain irradiance based upon using a terrain configuration factor.
$E_{t-\mathit{tcf}2}$	[W m^−2µm^–1]	Adjacent-terrain irradiance based upon using a terrain configuration factor, after Kondratyev (1977).
$E_{t-\mathit{tcf}3}$	[W m^−2µm^–1]	Adjacent-terrain irradiance based upon upon Vsky.
i	[radians]	Solar-terrain incidence angle.
I	[radians]	Hemispherical incidence angle.
It	[radians]	Adjacent terrain incidence angle.
L	[W m^−2sr${}^{-1}$ μm^–1]	Incident surface-reflected radiance.
$L^{\downarrow }$	[W m^−2sr${}^{-1}$ µm^–1]	Hemispherical downward diffuse radiance.
Lp	[m]	Length of the penumbra.
$\Delta z_{\max }$	[m]	Maximum relief.
S	dimensionless	Coefficient for cast shadows, fraction of Eb incident on the landscape.
Sb	dimensionless	Coefficient for terrain blockage for Et.
St	dimensionless	Coefficient for terrain shielding for Ed.
$T^{\downarrow }$	dimensionless	Total downward transmittance.
$T_t^{\downarrow \uparrow }$	dimensionless	Total atmospheric transmittance due to terrain relief.
Vsky	dimensionless	Skyview-factor coefficient.
$\overline{z}$	[m]	Average altitude.