Numerical Modeling and Parameter Sensitivity Analysis for Understanding Scale-Dependent Topographic Effects Governing Anisotropic Reflectance Correction of Satellite Imagery

Anisotropic reflectance correction (ARC) of satellite imagery is required to remove multiscale topographic effects in imagery. Commonly utilized ARC approaches have not effectively accounted for atmosphere-topographic coupling. Furthermore, it is not clear which topographic effects need to be formally accounted for. Consequently, we simulate the direct and diffuse-skylight irradiance components and formally account for multi-scale topographic effects. A sensitivity analysis was used to determine if characterization schemes can account for a collective treatment of effects, using our parameterization scheme as a basis for comparison. We found that commonly used assumptions could not account for topographic modulation in our simulations. We also found that the use of isotropic diffuse irradiance and a topographic shielding parameter also failed to characterize topographic modulation. Our results reveal that topographic effects govern irradiance variations in a synergistic way, and that issues of ARC need to be formally addressed given atmospheretopography coupling. Collectively, our results suggest that empirical ARC methods cannot be used to effectively address topographic effects, given inadequate parameterization schemes. Characterizing and removing spectral variation from multispectral imagery will most likely require numerical modeling efforts. More research is warranted to develop/evaluate parameterization schemes that better characterize the anisotropic nature of atmosphere-topography coupling.


Introduction
Systematic trends and variability in atmospheric conditions are causing significant changes to the Earth's climate, cryospheric, ecological, hydrological and geomorphological systems [1][2][3][4]. Mountain environments are very susceptible to rapid environmental change caused by the coupling of these changing systems, which also governs geohazard potential [5]. Unfortunately, due to practical and scientific issues, we do not fully understand the dynamic nature of climate forcing and coupled systems that governs mountain geodynamics and regulates environmental change [6,7]. We do know, however, that progress will require numerical modeling of key surface parameters and coupled systems, and quantitative and thematic information derived from remote-sensing observations.
Numerous surface processes are regulated by the surface irradiance, including evapotranspiration, surface temperature, and ablation. Given the paucity of suitable geospatial irradiance data at the local scale, numerical modeling efforts must be used to account for the scale dependent topgographic effects that govern radiation-transfer processes that dictate the magnitude and variability of surface irradiance. Furthermore, there are numerous parameterization schemes that can be used to simulate radiation-transfer parameters (RTP) and the irradiance components, although various parameterization schemes do not consistently make use of the same parameters and/or make specific assumptions (e.g., isotropic versus anisotropic; parameter insignificance, lack of formalization). It is important to note that many RTP are governed by multi-scale topographic effects that represent first-, second-, or high-order controls on the magnitude and variability of surface irradiance. Unfortunately, there is a paucity of research that evaluates variations in the magnitude of influence of various RTP on the magnitude and variance of irradiance components. Ultimately, it is essential that parameterization schemes can inherently account for the unique scale-dependencies and couplings associated with specific radiation-transfer processes and the radiation-transfer cascade (RTC).
In remote sensing, it is widely known that extracting accurate information from satellite imagery acquired over mountain terrain is notoriously difficult because of multi-scale topographic effects [8][9][10]. Anisotropic reflectance correction (ARC) is required, and the use of semi-empirical ARC methods have been routinely applied [9]. Bishop et al. [11] demonstrate and discuss the significant limiting issues associated with these empirical approaches, as parameterization schemes do not account for important RTP or the scale dependencies associated with various processes of the RTC. Other investigators have also recognized this and have developed and evaluated new parameterization schemes in at attempt to better account for cast shadows, diffuse-skylight irradiance, the Bidirectional Reflectance Distribution Function (BRDF), and the adjacent-terrain irradiance (e.g., [12,13]). Unfortunately, we do not adequately understand and know with certainty which RTP must be accounted for in ARC parameterization schemes to effectively remove multi-scale topographic effects and atmosphere-topography coupling from multispectral imagery, given geographical variations in the topographic complexity that governs a particular mountain environment.
The objective of this research is therefore to more formally characterize multi-scale topographic effects and evaluate various RTP that regulate the magnitude and spatial variation in surface irradiance. Our sensitivity analysis of control parameters should provide insights into which topographic effects may be important to account for in ARC parameterization schemes. We conduct our research over the Nanga Parbat Massif of northern Pakistan, that exhibits extreme topographic complexity in terms of local and mesoscale topographic properties and land cover conditions. More specifically, our research objectives are to: 1.
Simulate the direct irradiance and conduct a sensitivity analysis based upon the inclusion/exclusion of numerous RTP that are assumed to be significant/insignificant. We evaluate solar geometry, atmospheric attenuation, local topography, and cast shadows, where parameterization schemes are compared based upon different assumptions and parameter dependencies.

2.
Simulate the diffuse-skylight irradiance and conduct a sensitivity analysis based upon the inclusion/exclusion of numerous RTP that are assumed to be significant/insignificant. We evaluate the isotropic versus anisotropic assumption, the influence of secondary ground reflectance, local topography, and meso-scale relief, where parameterization schemes are compared based upon assumptions and parameter dependencies.
This research addresses a fundamental issue in ARC where the synergistic influence of topographic effects associated with the RTC have not been effectively characterized and evaluated. Consequently, our results should provide new insights into the nature of scale dependencies, atmosphere-topographic coupling and RTP that should be considered in ARC parameterization schemes.
It is important to realize that the inherent complexity associated with addressing such a large number of RTP, and the need to develop and evaluate new parameterization schemes, dictates that we limit the scope of this investigation. The nature of this research is multifaceted and we do not address issues associated with the temporal dimension, the BRDF, simulating the entire RTC, or sensor-system modulation of topographic effects. Nor is it the objective of this research to evaluate ARC methods or to develop new ARC methods. Rather, this research attempts to serve as a foundational basis for better understanding topographic effects of the dominant irradiance components and their potential for characterization in ARC methods.

Background
The surface irradiance (E) is strongly controlled by multi-scale topographic effects [14] and they partially regulate the the magnitude and spatial distribution of numerous surface processes. These processes govern landscape erosion and deposition dynamics, thereby altering topographic properties which then modulate the magnitude of surface irradiance, in complex radiation-topography-process feedbacks and system couplings that govern mountain geodynamics [15,16]. Therefore, radiation, precipitation, topographic and tectonic forcings can facilitate rapid environmental change in mountain environments. Consequently, the magnitude and partitioning of irradiance components and associated scale dependencies needs to be better investigated given the inherent limitations associated with using regional and/or global climate models [7].
The partitioning of these components with respect to their contribution to E are generally thought to systematically decrease in significance, although the spatio-temporal scale dependencies related to RT processes and various atmosphere-topography couplings have not been adequately accounted for. It is well known that: where E b is the direct irradiance from the solar beam, E d is the diffuse-skylight irradiance due to atmospheric scattering, E t is the adjacent-terrain irradiance, and λ is wavelength. Many Earth-science studies do not adequately account for these components at local scales many RT parameters are not formalized at the appropriate scale and do not account for important topographic effects [7]. Similarly, remote-sensing investigations involving empirical ARC methods involving the use of a scaling factor for topographic correction do not usually account for numerous topographic effects that govern E b , E d and E t , as the focus is on accounting for local topographic variations associated with E b (e.g., [9,[17][18][19][20]). Consequently, research involving the development and evaluation of new RTP schemes for these methods is sorely needed.

Direct Irradiance
The direct irradiance is the dominant irradiance component that is governed by orbital dynamics, atmospheric-topography interactions, and local and meso-scale topographic effects. It is computed as: where E 0 is the exoatmospheric irradiance, T ↓ is the total downward atmospheric transmittance that accounts for atmospheric attenuation, cos i is the cosine of the incidence angle (i) that accounts for solar and terrain geometry relationships that governs terrain self shadowing, and S is the shadow coefficient (0.0-1.0) that should account for cast shadow regions caused by the meso-scale relief structure of the terrain. It should be noted that direct shadowing is accounted for by the factor S. Remote-sensing investigations involving empirical ARC approaches, that are based upon the use of a scaling factor, do not usually account for the orbital dynamics that govern temporal changes in the Earth-Sun distance and variations in the solar zenith angle spatially (i.e., small-angle approximation; [11]). Similarly, many atmospheric processes may not be considered such as optical-depth variations, atmospheric refraction, and absorption and scattering of radiation for various atmospheric constituents, that represent an atmosphere-topography coupling effect. Such topographic effects must be accounted for, although we do not know with certainty the significance of various effects given environmental variations [11].
Most ARC methods utilize the cos i parameter, as variations in local topographic conditions can vary significantly, and this parameter is strongly correlated with surface reflectance variations in imagery. Nevertheless, the solar zenith angle at the center of a scene is most typically utilized to compute this parameter, although the solar geometry changes for every pixel across the scene, and it is not known with any certainty how the variance structure might be different given a relatively large scene, although we intuitively know that irradiance variation should increase as the size of the geographic area increases.
Ray tracing can be used to characterize cast shadows and investigators typically assume a point-source of light from the sun. Most ARC parameterization schemes do not account for cast shadows, although researchers are beginning to demonstrate its significance for accounting for irradiance and spectral variation in imagery (e.g., [10,21]). Unfortunately, the point-source assumption is inadequate to account for irradiance variations in the penumbra region, and numerical models and ARC methods need to account for such variations generated by the solar disk, as the umbra does not receive any E b (S = 0.0), while the penumbra region receives irradiance (S > 0.0; S = 0.0-1.0; [16,21]).To our knowledge, investigators have not evaluated the difference between these two formalizations for cast shadows, and whether the RTP should be used to account for meso-scale topographic effects. Collectively, these RTP are modulated by location, time and topography, and their significance for accurately estimating E b needs to be evaluated.

Diffuse-Skylight Irradiance
Atmospheric scattering will produce hemispherical anisotropic irradiance [22][23][24]. A two-stream approximation model that is widely utilized for computing the skylight irradiance for a horizontal surface (E dh ) assumes that: where E dr is the Rayleigh scattering component, E da is the aerosol scattering component, and E dg is the ground/sky backscattering component caused by multiple interactions between the ground surface and the atmosphere [23,25,26]. One can account for single scattering or a multiple-scattering Rayleigh atmosphere. The single-scattering albedo is required to account for aerosol scattering which is a function of both wavelength and humidity. Aerosol type and models such as rural, urban, maritime or tropospheric can also be accounted for using parameterization schemes or recent radiation-transfer models including MODTRAN and DISORT [25]. Finally, the ground backscatter component can be modeled by accounting for the zonal ground reflectance from the direct irradiance, the reflectance from the diffuse irradiance, and the overall sky reflectance [25]. Equation (3) does not account for anisotropic skylight irradiance or local and terrainshielding topographic effects. To represent these factors, E d must be computed for every pixel to account for scale-dependencies at each location on the landscape. We provide a computation solution such that: where L ↓ is the downward radiance from incident directions θ i and φ i , θ i and φ i are the zenith and azimuth angles of the incident energy from the hemisphere, I is the incidence angle of the direction defined by hemisphere and terrain geometry (similar to the computation of cos i), and S t is a binary terrain-shielding coefficient which is 0.0 given terrain shielding, or 1.0 at relative high hemispherical zenith angles. The E d component is generally not accounted for in ARC investigations, and proxy parameters may be utilized to ineffectively represent it (e.g., [18,20]). It has been accounted for by some investigators (e.g., [27][28][29]), although they utilize a parameterization scheme that makes use of Equation (3) and the skyview coefficient, which attempts to represent hemispherical topographic shielding. Local topographic effects and the anisotropic nature of the hemispherical irradiance distribution are not usually accounted for.
The anisotropic nature of E d has been the largest source of error associated with estimating this component, as it is governed by circumsolar brightening due to forward scattering of aerosols, and horizontal brightening due to multiple Rayleigh scattering [22,24,25]. An exact computation is computationally expensive given that E d is wavelength dependent, and the hemispherical anisotropic coupling of irradiance and local and meso-scale topographic shielding must be accounted for. Sensitivity analysis associated with using various parameterization schemes and combinations of RTP is sorely needed to provide insights into the use of this parameter for investigating topographic effects.

Adjacent-Terrain Irradiance
The aforementioned irradiance components and surface BRDF interact with the surrounding terrain geometry to generate E t . Numerous researchers have indicated that it should be accounted for in mountains due to snow, glaciers, vegetation and extreme relief over short distances [14,23]. It is an extremely complicated parameter to compute, as it is governed by atmospheric attenuation and numerous multi-scale topographic effects and surface anisotropic reflectance conditions. A first-order approximation to E t was formulated by Proy et al. [23], assuming Lambertian reflectance and is still used in investigations that attempt to account for E t . Nevertheless, it is rarely accounted for in empirical ARC studies.
It can be computed as: where the superscript e represents the effective zenith and azimuth angles that account for the influence of the terrain slope and slope azimuth angles on the incident (θ e i , φ e i ) and viewing geometry (θ e v , φ e v ) of the surface BRDF compared to a horizontal surface. L represents the surface reflected radiance coming from the effective incident direction, ρ brd f is the surface BRDF, φ i is the hemispherical incident azimuth angle, θ i is the incident vertical hemispherical zenith angle to account for terrain radiance above and below a pixel location, T ↓↑ t is the atmospheric transmittance given the optical depth of the atmosphere due to relief and propagation zenith angle through the atmosphere (θ v ) between two locations, cos I t represents the terrain incidence angle given the influence of the local terrain geometry in relation to the incident directional geometry, and S b represents the blockage of the adjacentterrain reflectance between two points that occurs when crests and/or ridge altitudes extend above the altitude of the 3-D trajectory connecting two points. Potential terrain blockage is direction and distance dependent and is different from the terrain shielding effects of S t . It is generally thought that the computation of E t should account for the terrain conditions extending out to approximately 5 km from each pixel location [23], although the scale dependency of contributing terrain varies with the complexity of the topography and relief, and most likely is highly variable across the landscape. Assessing the scale dependencies for each pixel location has not been adequately investigated. Another complicating factor is that each pixel exhibits a unique BRDF, as topography and varying land cover structure and biophysical properties govern the BRDF [11,30]. Consequently, accurate estimates of the magnitude and the spatial distibution of E t need to be assessed in the field and through numerical modeling efforts.
Most empirical ARC investigations do not adequately account for Et, as they assume Lambertian reflectance, insignificant atmospheric attenuation, and may not account for anisotropic topographic effects. Such assumptions are not valid in many mountain environments, as steep slopes and highly variable biophysical properties caused by high-magnitude surface processes, create anisotropic reflectance conditions and extreme relief [17,31,32]. Furthermore, topographic effects on the surface BRDF are not fully understood, although we know that the viewing geometry can change for every pixel given variations in terrain slope and slope-azimuth angles [13,26,33]. Furthermore, accurately estimating E t requires accurate characterization of the anisotropic irradiance conditions of E b and E t which govern the anisotropic nature of the BRDF. The problem represents the issue of the "chicken or the egg", as E t is dependent upon knowing the BRDF, but the BRDF is governed by irradiance components, as partially described in Equation (5a). More research regarding assessment of the significance of these RTP are needed to determine the most important parameters to be included in ARC methods or models.

Study Area
We simulate a first-order approximation of E t and E d over the Nanga Parbat massif, in northern Pakistan (see study area figure in Bishop et al. [11]). This geographic region is excellent for evaluating multi-scale topographic effects on the irradiance components, as it exhibits the full range of topographic properties due to atmosphere, surface process and tectonic interactions that are responsible for high magnitude erosion, uplift and relief production [34][35][36]. We specifically simulate irradiance over a 60 × 60 km area around Nanga Parbat to account for relief structure, meso-scale topographic shielding and surface conditions.
The study area has been extensively described in terms of climate, topographic, land cover, landforms, and geological conditions. Landscape conditions are spatially complex and dictate anisotropic reflectance conditions. Unfortunately, surface irradiance conditions are not well known for this area, as the multi-scale topographic effects have not been adequately accounted for. For excellent study area details see Shroder and Bishop [34], Zeitler et al. [35], Bishop et al. [36], Schneider et al. [37], Owen et al. [38].

Data
The ASTER Global Digital Elevation Model Version 3 (ASTER GDEM; [39]) was utilized to represent the topography and account for various multi-scale topographic effects governing surface irradiance components. Other data used for simulating E b and E d include spectra and constants from published work. Mean Exo-atmospheric irradiance values, atmospheric absorption coefficients for atmospheric constituents, and atmospheric scattering coefficients from Gueymard [25] were used in simulating atmospheric conditions and irradiant components.
We simulated surface reflectance based on the spatial structure of generalized land cover conditions. We used Landsat 8 OLI imagery captured in 2018 on 4 August (path 149, rows 35-36) and 9 September (path 150, rows [35][36], minimizing cloud cover and temporal differences between scenes. To assist in classification efforts, we computed the normalized difference vegetation index (NDVI), the normalized difference snow index (NDSI), modified normalized difference water index (mNDWI) as: where Green, Red, N IR and SW IR represent the green, red, near-infrared and shortwaveinfrared region of the spectrum using appropriate Landsat 8 spectral bands. We chose four basic land-cover classes for our analysis: 1. Snow, where NDSI > 0.75 and altitude > 3000 m.

2.
Vegetation, where there was no snow and NDVI > 0.1.

3.
Water, where there was neither snow nor vegetation and NDSI > 0.79.

4.
Rock and sediment in the remaining space.
To achieve a more homogeneous result, we removed all solitary water pixels and iteratively applied a majority filter to the snow and vegetation classification maps until no reclassifications occurred. We then simulated surface reflectance (i.e., assuming Lambertian reflectance) for these classes using the same linear spectral-mixing method and spectral libraries as Bishop et al. [11].

Orbital Dynamics
Earth-Sun orbital parameters control the magnitude of the exoatmospheric irradiance (E 0 ), therefore, we account for eccentricity (e), obliquity ( ) and other parameters such as the longitude of perhelion [40]. We used Berger [40] amplitudes, rates and phases of trigonometric expansions to predict these orbital parameters to compute the Earth-Sun distance (d es ), the distance-correction factor ( f r ) and solar geometry, such that: whereĒ 0 represents the solar irradiance spectra from Gueymard [25] at 1 AU and θ s is the solar zenith angle. Orbital parameters are also required to compute the spatial variation in the solar geometry such that: where δ is the solar declination, and λ l is the true longitude of the Earth relative to the vernal equinox. The solar zenith angle and solar azimuth angle φ s vary for each pixel and govern a multitude of parameters. The geocentric solar zenith angle is: where ϕ is latitude, and H is the hour angle of the sun. We modified this angle using a parallax correction. Parallax correction accounts for the radius of the Earth at a particular latitude, the height relative to the ellipsoid, and the distance from the sun. Atmospheric refraction also modulates θ s and we evaluated several refraction algorithms given that there are numerous parameterization schemes that account for different parameters (Table 1). Specifically, we evaluated the tan 5 formula, as presented in Corbard et al. [41], and the parameterization scheme of United Kingdom Hydrographic Office [42], and compared refraction estimates to measured data and full integration results produced by Corbard et al. [41]. In our simulations, we used the tan 5 formula, as it produced better results at higher observed zenith angles. The solar azimuth angle for each pixel was computed as: where the sine and cosine components are computed as: and where α s is the solar elevation angle. The solar azimuth angle was then corrected for grid convergence which is a function of latitude and the longitude of a pixel location with respect to the central meridian of the projection used (i.e., Transverse Mercator).  7 , and h = 1323 m, where λ represents wavelength, T is temperature, p is pressure, RH is relative humidity, ϕ is latitude and h is altitude. Data based upon Corbard et al. [41]. tan 5 represents the truncation of the expansion in the odd power-of-tangent refraction model that accounts for topographic, atmospheric and wavelength parameters, and AA represents the refraction formulation in United Kingdom Hydrographic Office [42] that only accounts for temperature and pressure. The parameters γ 1 and γ 2 represent tan 5 − FI and AA − FI, respectively.

Atmospheric Transmittance
We compute the downward total transmittance (T ↓ ) based upon primary atmospheric constituents and processes. These include Rayleigh scattering transmittance (T r ), aerosol scattering transmittance (T a ) using a rural clear-sky aerosol model, and atmospheric absorption transmittance of water (T w ), ozone (T o ) and primary gases (T g ). The parameterization scheme for ozone was based upon Bird and Riordan [43], while all the others were based upon Gueymard [25]. The atmospheric parameters used in these parameterization schemes to simulate atmospheric transmittance and compute irradiance components are based upon the U.S. standard atmospheric model for our location (Table 2). Finally, the total downward atmospheric transmittance was computed as: For more details regarding atmosphere transmittance parameters, see Bishop et al. [11].

Direct Irradiance
We computed E b using Equation (2). We do not formally define the first three parameters, as they have already been defined or are widely know and described in the literature (e.g., cos i). We focus here on the cast shadow parameter, as most studies utilize the solar point-source assumption and ray tracing to compute a binary coefficient (e.g., [44][45][46]). This approach generates an estimate of the distribution of cast shadows, but it spatially over-estimates the distribution and fails to accurately characterize the magnitude of E b in the cast shadow region, as the penumbra region of the cast shadow is not accounted for.
Cast shadows are governed by the size of the solar disk in relation to the Earth-Sun distance. More specifically, the solar angular width (α(t), [radians]) in relation to topographic relief, will govern the overall length of a cast shadow in the direction of φ s . It can be computed as: where R n is the nominal solar radius (R n = 695,700 km, [47]) and d es is the Earth-Sun distance (1AU = 149,597,870.7 km) that is a function of time (t). Given the aforementioned values, the solar angular width is ≈0.53290 • . This angle and other cast shadow parameters are depicted in Figure 1. The shadow parameters include the length of the penumbra (L p , [m]) and the maximum potential planimetric length from the point of shielding for the umbra (L u , [m]) and the entire cast-shadow region (L s , [m]). They can be computed as [21]: Finally, the fraction of E b that is incident on the landscape in the penumbra region can be computed as: The umbra is the region over the landscape where E b is totally obstructed by the topography. Conversely, the penumbra region receives a fraction of E b because a portion of the solar disk is obstructed by the topography. As depicted in Figure 1, the umbra region would exhibit shadow coefficients of 0.0, while the penumbra region would exit a gradient of coefficient values ranging from 0.0-1.0, depending upon the positional distance within this zone. Consequently, cast shadows also include irradiance from E b [21]. To our knowledge, this characterization of cast shadows has not been evaluated compared to a point-source characterization and the degree to which this RTP may influence ARC.

Diffuse-Skylight Irradiance
We computed E d using a variety of parameterization schemes. Given that E dh from Equation (3) represents the isotropic background diffuse irradiance component for numerous parameterization schemes, we compute the governing equations for the components of E dh such that [25,43]: where F r is the downward fraction of scattered radiation for a single scattering atmosphere, T aa is the transmittance of the aerosol absorption process, F a is the downward fraction of scattered flux, T as is the transmittance for aerosol scattering, E bn is the direct normal irradiance, ρ sky is the sky reflectance and ρ g is the ground reflectance averaged over a radius of 5 km. A commonly used parameterization scheme for estimating this component that attempts to account for topographic shielding (E dt ), involves the calculation of the skyviewfactor coefficient (V sky ) such that: where φ represents the hemispherical (i.e., skydome) azimuth angle, θ max represents the maximum horizon angle in the direction of φ and over a distance of d. This scheme assumes isotropic irradiance conditions and does not account for local topographic effects. We used a distance of 20 km around each pixel and used a φ interval of 5 • to compute the coefficient.
In this way, we can compare this parameterization scheme to a more comprehensive parameterization scheme (Equation (4)) that better accounts for topographic effects.
We also wanted to examine a widely used anisotropic diffuse parameterization scheme (E dp by Perez et al. [48]) that has been evaluated in the literature (e.g., [49,50]). It is computed as: where E dh represents the isotropic background, a 1 /a 2 is the solid angle for the circumsolar region, F 1 is the horizontal brightness, and F 2 is the circumsolar coefficient, where isotropic conditions represent F 1 = F 2 = 1.0. This scheme accounts for circumsolar and horizon scattering as well as local topographic conditions. It does not, however, formally account for hemispherical terrain shielding, although an approximation term ( 1+cos θ t 2 ) is used to represent this meso-scale topographic effect.
Our final and more comprehensive parameterization scheme for E d that we compute is based upon Equation (4). Specifically, we perform sky luminance modeling as described in Perez et al. [24], to estimate the anisotropic nature of the diffuse irradiance and account for horizon scattering and circumsolar scattering. We assume clear-sky conditions and specifically utilize the CIE clear-sky formulation [51] to compute the hemispherical sky illuminance. We then generate normalized luminance value for all the elements of the sky dome, and use the Perez et al. [24] diffuse luminous efficacy model to convert values to irradiance. We then integrate over the entire skydome accounting for local topographic effects that are characterized by skydome and terrain gemometry, and limit the integration based upon the θ max angle for a particular hemispherical azimuth direction, thereby accounting for hemispherical topographic shielding. In this way, we formally account for anisotropy and the scale dependent nature of topographic shielding.

Sensitivity Analysis
For each irradiance component, we conducted a numerical sensitivity analysis to determine if a particular RTP, or a different parameterization scheme, significantly governs the magnitude and variance structure of a particular irradiance component compared to a more comprehensive parameterization scheme that is deemed to be more accurate (i.e., accounts for more topographic effects). This was done to identify significant RTP that may need to be accounted for in ARC methods or models. Furthermore, we simulated results for 3 different spectral wavelengths based upon the central spectral bandwidth wavelength of the Landsat 8 Operational Land Imager sensor (i.e., λ 1 = 0.56141 µm, λ 2 = 0.65459 µm, λ 3 = 0.86467 µm).
It is important to note that we evaluate three specific wavelengths and characterize specific atmosphere-topography couplings. Spectral integration would have generalized our results, and our focus is to evaluate parameterization schemes to see how well they can characterize topographic effects. It is well known that the spatial, spectral and radiometric resolution of sensor systems will modulate the magnitude of the recorded sensor signal from the radiation-transfer cascade (RTC), such that environmental information (i.e., topographic variations) will be modified, based upon RTP dominance in relation to the degree of generalization cause by various resolutions, which alters spectral variation in images. We do not yet have enough knowledge about how sensor system characteristics and spectral integration will influence topographic variation in imagery, so we did not want to complicate the problem of generalizing our results over a portion of the spectrum. We primarily focus on terrain-irradiance variations and we do not compare irradiance images to spectral imagery.
For all irradiance components, we included or excluded specific RTP into our simulations and analyses that represent a series of control parameter scenarios, so that we can examine statistics and compare control parameter scenarios to a more comprehensive parameterization scheme over the Nanga Parbat Massif. 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 globally computed structural similarity index (SSI) [52], Students t-test analysis for unequal variances, and an F-Test to determine the significance of a parameter with respect to the irradiance component. See Table 3 for simulation and scenario descriptions. Table 3. Simulation and control scenario (SS) descriptions for generating estimates of irradiance components. Control parameters (P c ) for a scenario were included ( i ) or excluded ( x ) for comparison to a more fully comprehensive parameterization scheme.

SS P c Description
Solar zenith angle only at scene center. S1S2 P x Parallax correction omitted. S1S3 R x Atmospheric refraction correction omitted. S1S4 P x , R x Parallax and refraction corrections omitted. S1S5 T ↓,x Total downward Atmospheric transmittance term omitted. S1S6 (cos i) x Cosine of the incidence angle omitted. S1S7 S x a Cast shadows using solar disk omitted. S1S8 S i p Cast shadows using point-source assumption. S1S9 E b Base for comparisons using Equation (2). It should be noted that the purpose of this study is to evaluate a limited number of key variables that represent topographic effects, and those that could be significant in terms of improving parameterization schemes of ARC methods that make use of a scaling factor to remove topographic effects in satellite imagery. In fact, our choice of variables is related to the assumptions that are frequently associated with utilizing empirical ARC methods and the choice of possible parameterization schemes that could be attempted in further research. In this way, our evaluation may provide insights into possible improvements for modifying empirical ARC parameterization schemes by including or excluding variables. The use of evaluating specific RTP parameters/sub-parameters based upon a Monte-Carlo simulation approach is far beyond the scope of our work, and this approach would dictate that a far greater number of variables would have to be evaluated, due to parameter function dependencies.

Direct Irradiance
A sample of cast shadow simulation results based on the solar point-source assumption for the upper region of the Diamir Basin are presented in Figure 2A. Figure 2B depicts a shadow distribution based upon using the solar disk to estimate the umbra and penumbra regions. The umbra region in the upper panel (A) is greater in length compared to the umbra region in the lower panel (B), although the spatial distributions are relatively similar. This is caused by θ 1 being less of an angle compared to θ s , as α/2.0 must be taken into consideration when accounting for the partial blockage of the sun. The penumbra region borders the extent of the umbra region, and various locations on the landscape within the general umbra region can also exhibit modulated irradiance due to relief conditions. The parameterization schemes dictates that changes in the Earth-Sun distance, solar geometry and relief can increase the length of the penumbra, thereby increasing the influence of this parameter on irradiance variability, and therefore spectral variability in imagery. Given the extreme relief at Nanga Parbat, we might expect this RTP to be significant in terms of governing the magnitude and variance structure of direct irradiance in various locations on the landscape. Direct irradiance simulation results are presented in Figure 3 and Table 4. As we would expect, local topographic conditions and cast shadows cause significant spatial variations in E b for the green region of the electromagnetic spectrum caused by the coupling of atmosphere-topography modulation. The magnitude and variance of E b for the red and near-infrared regions of the spectrum decrease given a systematic decrease in exoatmospheric irradiance and increase in atmospheric transmittance. These magnitude and spatial patterns are characterized in the E b statistics, where there is a systematic decrease in the mean and standard deviation for each region of the spectrum. These simulation results are consistent with results generated from Bishop et al. [11], although improved parameterization schemes are used in this study. Sensitivity analysis results of our direct irradiance simulations clearly demonstrate that all topographic effects that we have simulated govern the magnitude and spatial variability of irradiance in a significant way, irrespective of wavelength (Table 4). Specifically, we see relatively high RMSE values that indicate that specific topographic effects, if not accounted for, make a difference in the magnitude of irradiance and that specific parameterization schemes do not account for the collective topographic effects. There is clearly a wavelength dependence, such that RMSE values are much larger at shorter wavelengths, strongly suggesting that coupled atmospheric-topographic effects are the cause. This is also demonstrated by the decrease in the variability of irradiance with increasing wavelength (see σ values).
These results are also supported by the structural similarity index (SSI) that accounts for magnitude, global variability and correlation structure. We should expect an SSI value of one if the mean, variance and correction structure are identical. Clearly, SSI values indicate that there are magnitude, variance and correlation structure differences. The statistical significance of these differences can be ascertained by our inferential statistical results. We found that all inferential tests using the t-statistic and F-ratio were significant at the p = 0.00001 level for all regions of the spectrum, indicating that the magnitude and variance structure of direct irradiance is significantly different when RTP are not accounted for, or when different parameterization schemes are used. Collectively, these results clearly suggest that RTP that are generally thought not to be significant in terms of accounting for topographic effects may need to be accounted for, as one parameter can partially govern a multitude of other parameters considering the RTC. Table 4. Simulation and control scenario (SS) statistical results for direct irradiance. Control parameters (P c ) for a scenario were included ( i ) or excluded ( x ) for comparison to a more fully comprehensive parameterization scheme (E b ). Statistical parameter symbols include the minimum irradiance (Min, [W m −2 µm −1 ]), maximum irradiance (Max, [W m −2 µm −1 ]), mean irradiance (µ, [W m −2 µm −1 ], standard deviation (σ), root-mean-squared error (RMSE, [W m −2 µm −1 ]), 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 Appendix A (Tables 3 and A1). All inferential tests using the t-statistic and F-ratio were significant at the p = 0.0001 level for all regions of the spectrum, indicating that the magnitude and variance structure of direct irradiance is significantly different when parameters characterizing topographic effects are not accounted for, or when a different parameterization scheme is used. Our statistical results (Table 4) also provide insights into which RTP are most important. As we would expect, local topographic effects (i.e., cos i) are very significant, as solar and terrain geometry strongly governs irradiance and reflectance. This simulation and control scenario exhibited the highest and lowest RMSE and SSI values, respectively, irrespective of wavelength. Total atmospheric transmittance was also found to be important with relatively high RMSE values, while the third most important RTP appears to be cast shadows (solar disk parameterization scheme), while ARC methods make use of cos i, they generally do not account for the topographic effects that govern atmospheric transmittance or cast shadows (umbra and penumbra).

Diffuse-Skylight Irradiance
Simulation results demonstrate the anisotropic nature of E d over the sky dome and the anisotropic nature of topographic shielding (Figure 4). Panels A and B depict the entire sky dome and do not account for local or meso-scale topographic shielding.
The circumsolar irradiance zone governed by the solar geometry and the increase in scattering near the horizon are clearly depicted. The magnitude of E d is regulated by altitude, as lower altitude locations are associated with an increase in atmospheric mass that facilitates more atmospheric scattering. Panels A and B represent the lowest and highest altitudes over the Nanga Parbat Massif, with E d values of 446.895 and 301.155 W m −2 µm −1 , respectively. Local and meso-scale topographic effects would significantly decrease these magnitudes if they were taken into consideration, as topographic effects are highly location dependent.  Figure 4 demonstrate the anisotropic scale-dependent nature of topographic shielding. Local topographic effects were not accounted for in these simulations. We sampled locations in V-shaped and U-shape valleys to demonstrate meso-scale topographic effects. Panels C and D are located in the Diamir and Rakoit valleys, respectively, where rapid river incision and tectonic uplift produce well defined V-shaped valleys that have a strong terrain orientation. The orientation direction is a NW-SW trend for the Diamir Valley and a N-S trend for the Rakoit valley. The hemispherical relief structure in these basins clearly restricts the area of the skydome that contributes to E d . Panels E and F are located in the Rupal and Sachen valleys, respectively, where glaciation and glacierization have produced well-defined U-Shaped valleys. The relief structure also produces a terrain directional trend of NE-SW for the Rupal Valley and less of a directional trend for the Sachen Valley. Similarly, the hemispherical relief structure in these basins clearly restricts the area of the skydome that contributes to E d .

Panels C-F in
It is important to note that the anisotropic nature of the topographic shielding is relatively unique for each of our simulated sample locations, and this demonstrates that this topographic effect is highly variable over the landscape, such that each location has a relatively unique contributive skydome that is regulated by location (i.e., 3D position) in relation to the anisotropic relief structure of the topography. The relief structure is also highly variable across the landscape and strongly controlled by climate, surface processes and tectonics, which governs the orientation fabric of the topography. Consequently, E d is highly variable over the landscape, even though its magnitude may be less than E b .
Diffuse-skylight irradiance simulation results over Nanga Parbat are presented in Figure 5 and Table 5. As expected, the mean magnitude of E d decreases with increasing wavelength due to a decrease in atmospheric scattering. Simulated E d values exhibit considerable spatial variability caused by anisotropic diffuse irradiance, coupled with local and mesoscale topographic variations. Consequently, we note that the magnitude, variance and spatial distribution patterns of our E d simulations do not similarly correspond to E d estimates generated from other parameterization schemes that assume isotropic conditions, and those that do not formally account for local and mesoscale topographic effects. Our simulation results depict relatively high diffuse irradiance conditions on the southeastern slopes caused by relatively high circumsolar irradiance that is governed by the location and time of the simulation. Conversely, we see relatively low diffuse irradiance on the north side of Nanga Parbat, as diffuse irradiance is relatively low from the skydome geometry that faces the northwestern facing slopes. This is perhaps best depicted in the Indus valley, as the north-northwestern side of the valley (demarcated by the Indus River) exhibits relatively high irradiance, and the south-southeaster side of the valley exhibits relatively low irradiance. Such dramatic variation in irradiance over relatively short distances is also caused by local terrain geometry and topographic shielding where skydome irradiance is significantly modulated by these topographic effects. Furthermore, the significant relief within the region causes altitudinal variations in E d given an increase or decrease in atmospheric mass and therefore atmospheric scattering. These variations can be seen along the knife-edge Nanga Parbat ridge in the center of the scene.
Sensitivity analysis results of diffuse-skylight irradiance simulations clearly demonstrate that excluding a RTP and/or using a parameterization scheme that does not formally account for topographic effects, does not accurately characterize the magnitude and spatial variability of irradiance, irrespective of wavelength (Table 5). Specifically, we see that RMSE values are relatively high for control parameters, where we would expect a RMSE = 0.0 if anisotropy and multi-scale topographic effects were accounted for. We found that excluding the E a and cos I control parameters produced the highest RMSE values, irrespective of wavelength. This indicates the importance of the atmospheric aerosol scattering component and the coupled influence of anisotropic irradiance and local topographic effects. There is clearly a wavelength dependence, such that RMSE values are much larger at shorter wavelengths, strongly suggesting that significance of coupled atmosphere-topographic effects. This is also demonstrated by the decrease in the variability of E d with increasing wavelength (see σ values).
These results are also supported by the SSI values that accounts for magnitude, global variability and correlation structure. We should expect a SSI value of one if the mean, variance and correlation structure are identical. Clearly, all SSI values for control parameters are less than one, which indicates that all excluded RTP and included parameterization schemes exhibit differences in their magnitude, variance and correlation structure when compared to our base E d simulation, irrespective of wavelength. For excluded RTP, SSI values are lowest for the E a and cos I control parameters.
The significance of these differences can be ascertained by our inferential statistical results. We found that all inferential tests using the t-statistic and F-ratio were significant at the p = 0.0001 level for all wavelengths, indicating that the magnitude and variance structure of E d is significantly different when important RTP are not accounted for, or when different parameterization schemes are used. Collectively, these results clearly suggest that RTP that are generally thought not to be significant, or those topographic effects that are not generally considered (e.g., variation in solar geometry, secondary ground-atmosphere scattering, local topography, hemispherical shielding) may need to be accounted for, as one parameter can effect other parameters considering the RTC.
Our statistical results (Table 5) also provide insights into parameterization schemes that are frequently used by the Earth Science and remote sensing communities. As we might expect, the assumption of isotropic diffuse conditions, should not be used to characterize this irradiance component. The lowest SSI value was associated with E dh for all wavelengths. Coupled atmospheric-topographic effects are not effectively accounted for in this parameterization scheme. Consequently, researchers have also utilized the skyview factor to attempt to account for topographic shielding. Our results show, however, that this approach is not adequate, as SSI values for E dt are extremely low, irrespective of wavelength. Clearly, a topographic-shielding summary coefficient cannot represent the coupling of anisotropic irradiance and shielding conditions, as well as the coupling of anisotropic irradiance and local terrain conditions. Finally, we also note that a popular anisotropic diffuse irradiance scheme (E dp ) produces significantly different results because it cannot formally account for topographic shielding and uses a proxy parameter. Table 5. Simulation and control scenario (SS) statistical results for diffuse-skylight irradiance. Control parameters (P c ) for a scenario were included ( i ) or excluded ( x ) for comparison to a more fully comprehensive parameterization scheme (diffuse-skylight irradiance in bold). Statistical parameter symbols include the minimum irradiance (Min, [W m −2 µm −1 ]), maximum irradiance (Max, [W m −2 µm −1 ]), mean irradiance (µ, [W m −2 µm −1 ], standard deviation (σ), root-mean-squared error (RMSE, [W m −2 µm −1 ]), 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 Appendix A (Tables 3 and A1). All inferential tests using the t-statistic and F-ratio were significant at the p = 0.0001 level for all regions of the spectrum, indicating that the magnitude and variance structure of diffuse irradiance is significantly different when parameters characterizing topographic effects are not accounted for, or when a different parameterization scheme is used. It is widely known that topographic effects govern the direct irradiance that causes significant spectral variability in satellite imagery acquired over mountain environments. Our simulations demonstrate high spatial variability in E b that is governed by terrain variations in altitude, relief, slope and slope azimuth. In other mountain environments, the nature of terrain complexity will be different and irradiance magnitudes and spatial patterns will vary depending upon the magnitude and variability of these terrain properties. Researchers have been attempting to address this component for over 40 years, although they have not formally addressed topographic effects other than cos i.

SS
Our findings strongly suggest that other RTP that govern the RTC should be taken into consideration to account for spectral variation in imagery. Specifically, we found that a variety of topographic effects are significant and that they should be formalized in ARC parameterization schemes to better account for variations in E b . This is important, as researchers assume that the spatial variability in solar geometry, atmospheric refraction, parallax correction, atmospheric attenuation, and a generalized treatment of cast shadow do not independently govern significant variation in irradiance. Our sensitivity analysis demonstrates that they can, given the topographic conditions at Nanga Parbat. Furthermore, it is important to note that for cast shadows, the length of the penumbra region can be extended, as a function of time and relief, which can further increase E b variability compared to our static time simulation.
Specifically, the sub-parameters of T ↓ , and S a should be utilized. They are not accounted for in common ARC methods and are generally not utilized in modeling efforts because of the need for more complex modeling of atmospheric transmittance and the computational complexity associated with representing the umbra and penumbra regions. Clearly, given that the transmittance and cast-shadow parameters are coefficients, they could be integrated into a ARC parameterization scheme to account for the relative coupling of atmosphere-topographic effects, and the spatial and temporal dependencies associated with cast shadows. In this way, the relative variability in E b could be more comprehensively accounted for.

Diffuse-Skylight Irradiance
The application of commonly utilized ARC methods can result in overcorrection, which is thought to be the result of not accounting for the E d component [53,54]. This has prompted numerous investigators to develop new parameterization schemes in an attempt to address irradiance variability. Common approaches make use of a C parameter that was formally introduced as the C-Correction method [18]. As Bishop et al. [11] noted, however, a constant global parameter that is based upon empirical regression analysis cannot represent the irradiance variations caused by coupled atmospheric-topographic effects that includes local and meso-scale topographic variation. In essence, none of the commonly utilized ARC methods that are based upon scaling formally account for the E d component. This is problematic, as our results clearly demonstrate that E d can be highly variable over the landscape given significant variation governed by the nature of anisotropy associated with atmosphere-topography coupling.
Researchers have tried to compensate for this by utilizing other E d parameterization schemes in numerical modeling efforts, based upon assumptions of isotropy and proxy parameters to account for multi-scale topographic effects. Many researchers utilize the representation of E dh for a horizontal surface to represent this irradiance component, although our results demonstrate that the isotropic assumption cannot adequately characterize the anisotropic nature of this component coupled with multi-scale topographic effects. Attempts to compensate for this make use of a topographic-shielding coefficient (i.e., skyview factor), although we demonstrate that one coefficient cannot adequately account for the coupling of diffuse-irradiance anisotropy, local topography and shielding anisotropy. Consequently, all E dh related parameters could not account for the magnitude and spatial variability of our base E d parameterization scheme which accounts for this coupling. Even the popular parameterization scheme of Perez et al. [48] could not account for the coupling, as local topographic information is utilized to account for meso-scale topographic effects.
Our results also revealed the importance of local (cos I) and mesoscale (S t ) topographic effects that are required to adequately characterize atmosphere-topography coupling. The use of proxy topographic parameters to represent anisotropic topographic conditions should not be used. Unfortunately, full computational solutions are required to account for atmosphere-topography anisotropy.
This irradiance component can be integrated into modeling parameterization schemes designed to reduce topographic effects in imagery. The computational complexity involving integration suggest that it may be difficult to characterize this component in a relative way for incorporation into common ARC methods that make use of scaling coefficients, but it is clear that empirical regression analysis of imagery and the use of constant RTP cannot adequately account for atmosphere-topography coupling. Numerical modeling approaches to ARC are most likely required to account for anisotropy and scale dependencies associated with E d and E t .

Anisotropic Reflectance Correction
ARC represents a radiometric calibration approach that is required to remove topographic effects from multispectral imagery [11]. Researchers are still investigating common ARC methods that use scaling coefficients to accurately estimate surface radiance/reflectance, although the parameterization schemes do not account for a multitude of RTP that characterize surface irradiance and other topographic effects. Richter [8] and others have indicated that atmosphere-topography coupling needed to be accounted for, and our results demonstrate the need to account for other E b sub-parameters and the anisotropic nature of atmosphere-topography coupling associated with E d . It is also unclear the degree to which the E t component will have on image spectral variability, although we speculate that the variance in E t may be relatively high given the scale-dependent coupling of terrain and land cover conditions. Finally, we also speculate that the topographic effects governing the BRDF can cause significant variation in surface reflectance, as the anisotropic nature of E d , E t and land cover structure govern the BRDF. Common ARC methods do not fully account for surface irradiance and the BRDF. Furthermore, scaling methods result in under and/or over correction of surface radiance values that can increase or decrease (i.e., information compression) spectral variation in imagery.
The issue of information compression when using common ARC methods warrants additional research. Image information content is the result of the RTC, where the atmosphere and topographic effects govern irradiance and surface reflection, which is further modified by the atmosphere and sensor, such that a composite signal is recorded. Accurate scaling of the magnitude of surface reflectance in an upward or downward direction dictates that the collective topographic effects that govern lower or higher surface reflectance be accounted for. As previously indicated, ARC methods do not account for the collective topographic effects, given inappropriate parameterization schemes. Accounting for this from a spatial perspective, and assuming similar surface reflectance conditions, differential or anomalous scaling may result in an increase in spectral variation between a pixel and its local neighbors, as the magnitude of the scaling factors did not accurately account for the collective topographic differences between two locations. Conversely, assuming differences in surface reflection conditions, inaccurate scaling at different locations can cause a decrease in spectral variation, as the combination of over and under correction or various magnitudes of under or over correction can decrease spectral variation.
Ultimately, the scaling magnitude for a particular location needs to be based upon the collective topographic effects at that location and the differences between locations. ARC parameterization schemes cannot account for such differences, thereby resulting in problematic reflectance magnitudes and spatial patterns. Overcorrection is easy to identify when anomalous reflectance values exist, however, small/moderate magnitudes of spectral expansion and spectral compression may not be readily apparent from a visual interpretation perspective, and such conditions represent ARC errors. Investigators may incorrectly interpret spectral compression as removal of topographic effects, although it is not clear what topographic effects have been removed or remain (unseen), or if the compression represents a loss of surface reflectance information. More research involving the mapping of spectral expansion and compression locations is warranted, along with the development of new ARC parameterization schemes that account for more topographic effects, as revealed in this study.
Utilizing parameters that represent the relative influence of topographic effects is possible for the direct irradiance component related to atmospheric attenuation, local topographic effects and cast shadows. ARC parameterization schemes could be modified, but would need to be carefully evaluated. Using relative parameters to account for the complexities associated with representing E d and E t are more difficult, as they represent integrations of coupled conditions involving landcover, topography and atmosphere. Furthermore, accurate ARC will require accounting for the BRDF, and there is a strong interdependence between the BRDF and E t , as the BRDF must be known to accurately compute E t , as it is highly scale dependent. Similarly, the BRDF is governed by surface irradiance anisotropy. Consequently, we suggest that numerical modeling approaches may be more effective for addressing the issues associated with ARC, and that the common ARC methods involving scaling, have very limited potential to account for the collective topographic effects, as diagnostics and accuracy assessment for these methods can be highly subjective.
The primary advantage of numerical modeling for ARC is that a multitude of topographic effects can be reasonably simulated to generate irradiance components and anisotropic estimates of the irradiance related to the BRDF. This can assist in simulating more accurate BRDF's when combined with using empirical BRDF models that are calibrated with satellite imagery. Furthermore, numerical modeling also enables accounting for variations in atmospheric constituents which governs atmosphere-topography coupling. Common ARC scaling methods do not account for such complexities of the RTC.
Simulated imagery of surface reflectance can also be compared to topographicallycorrected imagery. In theory, if the collective topographic effects have been correctly accounted for, there should not be any significant correlation between irradiance/synthetic images and corrected imagery, and the magnitudes of surface radiance from the synthetic image should be nearly identical to the corrected imagery. In practice, this is notoriously difficult to accomplish for the follow reasons: (1) ARC scaling methods can result in information compression reducing surface radiance variation, resulting in no correlation between corrected imagery and topographic effects; (2) Different environments exhibit varying degrees of topographic complexity and variation which results in the number and significance of RT parameters that must be accounted for in models to assess topographic effects-ARC methods effectively only account for one RT sub-parameter (i.e., cos i), and most models do not account for anisotropic irradiance parameters and/or the BRDF; (3) Model parameterization schemes may not effectively account for topographic effects, and different parameterization schemes may produce a different variance structure, as simplifications and generalizations are often assumed to accurately characterize topographic effects for a specific RT parameter. Ultimately, comprehensive ARC models need to be developed that take into consideration sensor spatial, spectral and radiometric resolutions, that can govern the nature of topographic effects recorded by imaging sensor systems. More research on sensor modulation of environmental information is sorely needed.

Conclusions
Anisotropic reflectance correction of satellite imagery is required to remove the collective multi-scale topographic effects that govern the anisotropic nature of atmospheretopography coupling of the radiation-transfer cascade in mountain environments. Commonly utilized ARC methods are still being investigated to effectively remove topographic effects, although ARC parameterization schemes and those from numerical modeling efforts have not fully accounted for a variety of topographic effects related to the dominant irradiance components, given the computational complexity associated with spatial analysis and coupled integration of radiation-transfer parameters. Furthermore, it is not clear which topographic effects need to be accounted for, given the variability in locational topographic complexity, and the cumulative synthesis that one topographic effect can have on other topographic effects.
Consequently, we developed radiative-transfer parameterization schemes that account for a more comprehensive treatment of multi-scale topographic effects that govern the direct and diffuse-skylight irradiance components. We then conducted a parameter sensitivity analysis to determine if specific topographic effects or other parameterization schemes can account for the collective topographic effects, using our parameterization scheme as a basis for comparison. We specifically tested many RT parameters that have not been used in commonly used ARC methods, and those that have been used to addressed the limitations associated with ARC in numerical modeling efforts.
As expected, we found that the typical assumptions associated with ARC related to the direct irradiance component, such as the small-angle approximation, atmospheric attenuation and refraction, and solar point-source approximation for characterizing cast shadows, all make a difference in accounting for topographic effects in our simulations. When they were not more accurately accounted for, simulation results could not account for the collective topographic effects and generated statistically different results in terms of irradiance magnitude and variability. Although the cos i parameter is widely utilized in ARC, our simulations and analysis strongly suggest that coupled atmosphere-topographic effects should also be taken into consideration in high relief environments, as well as the variability in the irradiance in the penumbra region of cast shadows. Given that the length of a cast shadow is dependent upon the Earth-Sun distance and relief structure variability, the influence of this parameter can vary temporally.
We also found that the common assumptions of isotropic diffuse-skylight irradiance and the use of a topographic shielding coefficient to represent anisotropic topographic effects could not account for the collective topographic effects in our simulations. We also found this to be true with the Perez anisotropic irradiance parameterization scheme, as it utilizes local topographic conditions as a proxy to mesoscale anisotropic topographic effects. Local topographic conditions were also found to be important in characterizing the collective coupling of atmosphere-topography anisotropy. Although this irradiance component is not adequately accounted for in common ARC methods and numerical modeling efforts, our simulations and analysis strongly suggest that it should be accounted for, as its variability can be relatively high in mountain environments with moderate to high relief, and as it partially governs the anisotropic nature and topographic effects on the BRDF, which also governs spectral variability.
Our results reveal that a multitude of topographic effects govern surface irradiance variations in a synergistic way, and that assumptions, approximations and issues of ARC need to be formally addressed given atmosphere-topography coupling. Collectively, our simulations and analysis strongly suggest that empirical ARC methods that are based upon scaling cannot be used to address the problem of anisotropy and multi-scale topographic effects in mountain environments, given inadequate parameterization schemes. Characterizing and removing radiation-transfer parameter variance components from multispectral imagery will most likely require radiation-transfer modeling efforts. More research is warranted on developing and evaluating new parameterization schemes that more accurately characterize the anisotropic and scale-dependent nature of topographic effects on the surface irradiance and the BRDF. More research investigating the adjacent-terrain irradiance is also warranted in an attempt to understand the influence of this parameter on spectral variability, and in an attempt to better characterize the anisotropic nature of the BRDF. The next phase of our research will be to evaluate the scale-dependencies of the adjacent-terrain irradiance and the topographic effects of surface irradiance variability governing the anisotropic nature of the BRDF.