Error and Uncertainty Degrade Topographic Corrections of Remotely Sensed Data

Chemical and biological composition of surface materials and physical structure and arrangement of those materials determine the intrinsic reflectance of Earth's land surface. The apparent reflectance—as measured by a spaceborne or airborne sensor that has been corrected for atmospheric attenuation—depends also on topography, surface roughness, and the atmosphere. Especially in Earth's mountains, estimating properties of scientific interest from remotely sensed data requires compensation for topography. Doing so requires information from digital elevation models (DEMs). Available DEMs with global coverage are derived from spaceborne interferometric radar and stereo‐photogrammetry at ∼30 m spatial resolution. Locally or regionally, lidar altimetry, interferometric radar, or stereo‐photogrammetry produces DEMs with finer resolutions. Characterization of their quality typically expresses the root‐mean‐square (RMS) error of the elevation, but the accuracy of remotely sensed retrievals is sensitive to uncertainties in topographic properties that affect incoming and reflected radiation and that are inadequately represented by the RMS error of the elevation. The most essential variables are the cosine of the local solar illumination angle on a slope, the shadows cast by neighboring terrain, and the view factor, the fraction of the overlying hemisphere open to the sky. Comparison of global DEMs with locally available fine‐scale DEMs shows that calculations with the global products consistently underestimate the cosine of the solar angle and underrepresent shadows. Analyzing imagery of Earth's mountains from current and future spaceborne missions requires addressing the uncertainty introduced by errors in DEMs on algorithms that analyze remotely sensed data to produce information about Earth's surface.

2 of 16 all gases is lower (Shugart et al., 2001). Combinations of drought and fire affect mountain forests and sources of water (Moody & Martin, 2001). The critical role that mountains serve as water towers and vegetation hotspots may change under climate change, contributing to hazards to people living in or relying on mountain resources (Kirschbaum et al., 2020).
The recent National Academies' Decadal Survey for Earth science and applications, Thriving on our Changing Planet, reflects these multiple concerns, with recommendations calling for observations "at scales driven by topographic variability" to reflect the heterogeneity of ecological, hydrological, and geological dynamics in Earth's mountains (National Academies of Sciences, Engineering, & Medicine, 2018). Investigating these processes via remote sensing requires spatial resolutions fine enough to characterize the variability, recognizing that the topography affects the reflected signals, thereby affecting the retrieval algorithms that interpret the state variables and fluxes of energy and mass.
Analysis of the topographic effect requires information in digital elevation models (DEMs) of the bare surface, usually but not universally meaning DEMs, as distinct from digital surface models (DSMs) that include vegetation, buildings, or other features. We consider two globally available DEM datasets: the NASADEM (Buckley, 2020) and the Copernicus DEM (European Space Agency, 2021), both distributed at a resolution of 1 arcsecond (∼30 m at the Equator). Locally or regionally, finer-resolution DEMs are available, so we consider three of those, which were derived by lidar, interferometric synthetic aperture radar, and structure-from-motion stereo photogrammetry from fine-resolution images. Our analysis considers the fine-resolution DEMs, in three different terrains, to provide the best assessment of the topographic effects on solar illumination geometry, and we compare those assessments to those derived from the two globally available datasets.
Characterization of the quality of DEMs typically assesses the vertical accuracy of the elevation. Uuemaa et al. (2020), through comparison of globally available products with fine-resolution lidar elevations, estimated root-mean-square (RMS) errors of 8-10 m for the NASADEM and TanDEM-X datasets (TanDEM-X is the primary source of data for the Copernicus DEM). Guth and Geoffroy (2021) compared several datasets with airborne lidar and ICESat-2 data and preferred the Copernicus DEM based on its ability to penetrate vegetation canopies and retrieve bare-Earth elevations.
However, the focus on elevation errors misses the effect of the topography on remotely sensed information in the wavelengths of the solar spectrum, which lies with the solar illumination geometry. The cosine of the local solar angle and the shadows cast by neighboring terrain are the most important variables for remote sensing of Earth's surface in the reflective domain. On clear days, most of the irradiance is direct, but the diffuse component is significant (∼30%) at the blue end of the solar spectrum. The surrounding landscape causes multiple reflections, which can be represented by the sky view factor-the fraction of the overlying hemisphere open to the sky. The topographic variables that affect the solar angles also affect the viewing angles from the sensor.
We therefore assess the DEMs based on their ability to provide insight into the ways that topography affects our ability to retrieve properties of the surface important to the study of Earth science. Fundamentally, retrievals that are sensitive to the magnitude of the spectral reflectance will be most affected. Examples include snow albedo (Bair et al., 2021;Painter et al., 2013) and ecosystem composition (Bogan et al., 2019). Retrievals that utilize the shape of the reflectance spectrum, characterized for example, by the spectral angle (Kruse et al., 1993), will be less affected but not entirely immune because the fraction of the incident irradiance that is diffuse versus direct is sensitive to the topography. Finally, retrievals that depend on the wavelength of absorption features do not depend on the magnitude of the reflectance. The primary example is mineral identification in soils and vegetation, which requires enough illumination to identify spectral features (Clark et al., 2003;Mulder et al., 2013).

Elevation Data
We consider two spatial resolutions of DEMs: fine and coarse. The coarse-resolution datasets are available globally, whereas the fine-resolution data are available selectively in specific locations. Table 1 summarizes the information sources for three fine-resolution and two global coarse-resolution datasets.
For the fine-resolution imagery, data are derived from three different methods: lidar altimetry, interferometric synthetic aperture radar, and structure-from-motion using fine-resolution commercial satellite imagery.
1. Airborne Snow Observatories Inc. (Painter et al., 2016) maps snow depth with lidar altimetry over drainage basins in the Western U.S., Switzerland, and Norway. The operation acquires elevation data during the snow-free summer and then periodically measures the snow-on elevation during the winter and derives snow depth by subtraction.  For the coarse resolution imagery, we used two global data sources at 1 arcsecond resolution distributed in geographic (latitude-longitude) format. In cropping to the boundaries of each fine-resolution area, we added 5 km to each edge to minimize edge effects in calculating topographic parameters.
1. We spliced 1° × 1° tiles from the NASADEM (Buckley, 2020) together because our areas of interest crossed latitude or longitude tile boundaries. The NASADEM combines information from the Shuttle Radar Topography Mission (Farr et al., 2007) and stereo-photogrammetry from ASTER imagery (NASA & METI, 2019). 2. We downloaded Copernicus DEMs (European Space Agency, 2021) that were spliced and distributed by Open Topography. The Copernicus DEM is derived from TanDEM-X imagery. Figure 1 shows the Copernicus DEM on the left and the ASO DEM on the right for the Carson River Watershed. The red rectangle in the left-hand image shows the area that the ASO DEM illustrates the detail of the topographic data at 3 m spatial resolution.

Notation
We selected or calculated the following variables for each grid point in each elevation dataset. 0, 0, , , and vary with date; the other variables are independent of date and thus the solar illumination. Deep snow can smooth the topography, but our comparisons of snow-off with snow-on elevations find only a few grid cells with significantly different slope and azimuth. At the 3 m spatial resolution of the Carson River DEM, the RMS difference in slopes between the snow-on and snow-off elevations is 1.8°, and the 99th percentile absolute difference is 6.8°. Resampled to 10 m spatial resolution, the RMS difference is 1.3°, and the 99th percentile absolute difference is 5.3°. 0, 0 Solar zenith and azimuth angles, 0 = cos 0 Cosine of solar illumination angle on a slope, set to zero for slopes that are in shadow, either by adjacent terrain or when is negative Spectral directional-hemispherical or bihemispherical reflectance, depending on subscripts (Schaepman-Strub et al., 2006) Fraction of incoming spectral irradiance that is diffuse Angle to the topographic horizon, upward from horizontal, in azimuth direction Spectral irradiance, incoming or reflected depending on subscript

RMS
Root-mean-square value ( ) = √ 1 ∑ =1 | | 2 Slope angle, upward from horizontal, and slope azimuth, south at 0°, eastward positive and westward negative, consistent with a right-hand coordinate system We assume the fine DEM is more accurate, particularly because variables derived over multiple points are compared to those derived from an individual location in the coarse DEM; therefore, the RMS of the difference between the coarse and fine estimates of a variable is considered the RMS error in the coarse-resolution data.
We calculated for seven dates between the winter and summer solstices, spaced so that the intervals between the solar declinations were equal ( Figure 2). For every date, we chose 10:45 in the local time zone to match typical mid-morning acquisition times of satellites: Pacific Standard (UTC-8:00) for the Carson River, Alaska Standard (UTC-9:00) for the Wrangell Mountains, and India Standard (UTC+5:30) for the Himachal Pradesh. Figure 3 shows cosine of solar illumination angles for the Himachal Pradesh on the seven dates in Figure 2. The pixels are more illuminated as the solar illumination angle gets closer to zenith, and the fraction of diffuse illumination decreases, thereby affecting the relationship between intrinsic and apparent reflectance.  Figure 4b shows a calculation of , the cosine of the solar illumination angle at the same date and time as the Landsat image, using elevation data from the NASADEM (Buckley, 2020). The cosines are calculated from the slope and aspect of the terrain and the solar zenith and azimuth angles on a flat surface. Where shaded by local horizons or by the slope itself, the cosines are set to zero. Superficially, the two images appear to match, allowing that some illuminated areas are dark because that surface material is dark, whereas shadows are dark even if the surface material is bright. The bright areas in the Landsat image correspond to highly illuminated pixels. However, the scatter density plot (Figure 4c) indicates some problematic values. The high reflectance values in the upper left corner of Figure 4c correspond to pixels either in shadow or with highly oblique solar illumination, indicating that the solar angle calculated from the DEM is wrong. The low reflectance values in the lower right corner of the scatter plot tell a similar but more ambiguous story. These dark pixels are well illuminated; they could represent a dark surface, or they might not truly be well illuminated. Figure 4d shows probability density functions (pdf) of the reflectance values in areas of low ( < 0.2) and high ( > 0.87) illumination (each threshold represents 14% of the image). Each pdf has a long tail. Those in the tail of the low illumination category indicate that is not correctly estimated and is too small. With a correct DEM, we would not see such high reflectance values at highly oblique solar illumination, because of the low values of incident irradiance at those locations.
Algorithms that retrieve land surface properties analyze the spectral "reflectance," broadly defined to cover the several possible angular reflectance configurations that Schaepman-Strub et al. (2006) articulate. In their study of the effects of surface roughness on snow albedo, Bair et al. (2022) defined intrinsic reflectance of a substance independent of effects of roughness or topography. The corresponding apparent value, as one might measure at a plot, incorporates artifacts caused by roughness or topography. For example, the intrinsic reflectance of clean snow in the visible wavelengths cannot drop to 0.2, but the apparent reflectance of shadowed snow can reach such low values. Corrected for atmospheric effects, measurements of spectral top-of-atmosphere radiance by a satellite sensor can be converted to apparent values of, for example, bihemispherical reflectance or bidirectional reflectance (Schaepman-Strub et al., 2006). Retrieval of a parameter of interest at Earth's surface using these data requires an estimate of the intrinsic spectral reflectance, interpreted by a combination of topographic information and the apparent value derived from the satellite sensor data (Brodrick et al., 2021). For some pixels, however, incorrect or imprecise topographic information could cause those retrievals of the surface properties to produce incorrect interpretations, an issue addressed in Section 5. This study characterizes the errors in the solar angles in globally available DEMs and recommends steps to mitigate these uncertainties in retrieval of Earth's properties in mountainous terrain.

Results
Tables 2 and 3 summarize results for all fine-and coarse-resolution datasets analyzed. Figures 5 and 6 illustrate examples of the results, comparing pairs of variables derived from a fine-and a coarse-resolution image. We include examples from each of the three study sites: Carson River Watershed, Himachal Pradesh in the Indian Himalaya, and Wrangell Mountains in Southeast Alaska.

Topographic Variables Independent of Solar Illumination
Variations in elevation across topography create sloping terrain, so we characterize each pixel by its slope S upward from the horizontal and its aspect A as the direction the slope faces. Slope and aspect combine with the solar angles to create variability in local illumination. The varying terrain also creates the view factors Ω , the fraction of the sky hemisphere open above a point. The view factor controls the re-reflection of solar radiation that strikes the surface and the fraction of the diffuse irradiance that reaches the surface. The view factor is also important in modeling the thermal infrared radiation in the mountains (Robledano et al., 2022). The terrain geometry affects the incoming irradiance and the reflected radiation, so the errors in elevation itself are less important than errors in slope, aspect, and view factor. Based on the differences between the fine-resolution and coarse-resolution DEMs, Table 2 shows the RMS error for elevation, differences in elevation between neighbors, slope, aspect, and view factor. Because the differences between the datum sources (Table 1) for elevation exceed 25 m and because we are mostly interested in the internal differences within an elevation grid, we subtract the mean elevation of each grid from that grid's values before calculating the RMS errors for elevation.
Errors in elevation are small fractions of the elevation values themselves, but the errors in slope and aspect indicate significant errors in the differences between elevations of neighboring points. Calculation of the slope S and aspect A of a topographic pixel considers the spatial derivative of elevation Z in two or more directions x and y (Dozier & Frew, 1990), which could be projection coordinates or longitude and latitude distances computed from the coordinates and the dimensions of the ellipsoid: Including the signs for the numerator and denominator separately enables calculation of aspect A over the full circle. From topographic data, the derivatives are calculated numerically. In a matrix Z with grid spacing Δℎ representing topography and the rows running west-east and columns north-south, the derivatives at point [ ] , as calculated by a central difference method, are: The grid spacing Δℎ is usually known accurately, so assessment of errors that affect topographic radiation depends on the error distribution of the differences between neighboring elevations. We estimated the RMS error of the differences by calculating the numerators of Equation 2 in each direction and then the hypotenuse of the x-and y-direction differences in each pixel.  elevations are smaller than the RMS errors in the elevations themselves, thereby indicating some spatial coherence in the elevation errors. Otherwise, if the RMS errors of the elevations were indeed independent, then the variance of the differences would be the sum of the variances in the elevations themselves (Weisstein, 2021) and the RMS error of the differences would be √2 × the RMS error of the elevations. However, the RMS errors of the differences are much smaller.
Results for the NASADEM and the Copernicus DEM are similar, but both show outliers that translate into outliers in calculating illumination angles.
The variability in the data indicates variation within the topographic grid. Figure 5 shows the scatter diagrams for the row in Table 2 that summarizes the statistics for the Copernicus DEM for the Carson River Watershed in the Sierra Nevada. The x-axes represent values from the fine-resolution ASO DEM, and the y-axes the values from the globally available Copernicus DEM. For elevation, the spread around the regression in Figure 5a is small. For the other variables, however, the spread is much larger. The outliers in the scatter plots for slope and aspect imply that outliers are present in the local solar angles. Figure 5c shows that some slopes less than 20° in the ASO 3 m DEM correspond to slopes greater than 40° in the Copernicus 1 arcsecond DEM, and conversely some slopes greater than 50° in the finer-resolution DEM correspond to slopes less than 20° in the Copernicus DEM. Similar differences occur in the aspects and view factors.
The regression lines in Figures 5b-5d for differences in elevation between neighbors, slope, and aspect are constrained to go through the origin. For elevation ( Figure 5a) the datasets do not use the same datum (  so the regression includes an intercept. For the view factor (Figure 5c) all values are above about 0.6, so that regression is constrained to go through (1,1) instead of (0,0). In all cases except elevation, the slopes of the regression lines that characterize the relationship between the coarse-and fine-resolution variables are less than 1.0, indicating generally that the Copernicus DEM and NASADEM slightly underestimate the magnitudes. The section on Bias in Table 2 therefore indicates a negative bias in the coarser-resolution datasets (i.e., a regression slope of 0.90 corresponds to a bias of −10%).
Aspect values and their RMS errors must be treated with caution, because aspect has negligible effect on solar radiation when the slope is small but a huge effect when the slope is steep. In our formulation, we follow the right-hand convention that 0° aspect represents south, from which eastward aspects are positive and westward aspects are negative.

Effect of Topography on Illumination and Reflection
The two crucial topographic variables in order of importance are , the cosine of the local illumination angle measured from normal to the slope, and Ω , the fraction of the hemisphere over a point that is open to the sky. The equation for Ω uses the horizon angles for all directions (Dozier, 2022b): For slopes facing in the direction toward the Sun, that is, cos( − ) ≥ 0 , the limits of integration [ 1,2] being constrained to those azimuths: (3) Figure 5. Detailed illustration supporting one row in Table 2 for the Copernicus digital elevation model (DEM) in the Carson River Watershed in the Sierra Nevada. The x-axes show data for the Airborne Snow Observatories (ASO) 3 m DEM; the y-axes show the same information derived from the Copernicus DEM, with both DEMs reprojected to a common size and projection. Aspect angles represent south as 0°, eastward positive, westward negative, and therefore consistent with a right-hand coordinate system. Regression lines in the figure and statistics in Table 2 are based on the whole topographic grid, but just 100,000 points are randomly selected for the illustrative scatter plots.
For the slopes where cos( − ) < 0 , the slope itself might obscure the horizon, so in integrating across those values with the limits of integration corresponding to those azimuths, for each azimuth is set to Over a flat unobstructed surface, Ω = 1 .
The local illumination angle is related to the topography and the solar illumination geometry as: The max function accounts for slopes facing away from the sun by setting = 0 in situations where the equation would yield < 0 . To account for points where neighboring horizons block the Sun, we also set = 0 where sin 0 ≥ 0 .

The variables
and Ω affect the relationship between the apparent reflectance of the surface and its intrinsic reflectance that would be measured independent of any topographic effects (Bair et al., 2022). The apparent reflectance of a topographic surface involves multiple reflections, especially for bright surfaces such as snow. Let indicate spectral reflectance, omitting a wavelength identifier, and as the fraction of the spectral irradiance that is diffuse. Set the initial irradiance on a horizontal surface to I. The spectral radiation that initially escapes into the overlying hemisphere without being re-reflected is:  Table 3 for the Copernicus digital elevation model (DEM). All axes show values of , the cosine of local illumination, varying with the dates that Figure 2 shows. Points along either the x-or y-axis identify locations that are shadowed in one DEM and illuminated in the other. Regression lines in the figure and statistics in Table 3 are based on all pixels in the data, but just 100,000 points are randomly selected for the illustrative scatter density plots. Note that the yellow (bright) values in the scatter density plots migrate to higher values of as the solar declination moves northward.
The equation assumes that is approximately isotropic averaged over the field of view, it neglects atmospheric attenuation within the topography, and it ignores variation in albedo and irradiance within the neighborhood. The superscripts designate the reflectance to direct versus diffuse irradiance. The right-most term inside the brackets accounts for reflected radiation within a point's field of view impinging on the point. The direct and diffuse spectral albedos might differ slightly, for example, for snow.
Not all the initially reflected radiation escapes into the overlying hemisphere. Instead, some of it re-reflects and eventually escapes or is trapped by the topography, in which case it is subject to internal reflection. At the first iteration, its value is: To account for multiple reflections, at each reflection the value of the incident radiation is multiplied by the fraction (1 − Ω) that accounts for the reflection remaining within the topography, the fraction Ω that escapes, and the intrinsic spectral reflectance. An orders-of-scattering approach to the multiple reflections lets some reflected radiation escape at each iteration n and some remains available for re-reflection: This series converges in a half dozen iterations because ( ) declines in proportion to (1 − Ω) . The apparent reflectance for the pixel is = ∑ ∕ .

Errors in Estimating μ S , the Cosine of Local Illumination
RMS errors and outliers in the topographic variables combine with the solar illumination geometry to propagate into the calculation of each pixel's illumination. The most important variable whose accuracy affects the interpretation of the remotely sensed signal is the cosine of the local illumination angle. The ratio ∕ 0 appears in Equation 5, but 0 is usually known accurately. The view factor Ω affects the diffuse irradiance from the sky and the internal reflections within the topography.
Therefore, the accuracy of the cosine of illumination from the DEM affects our ability to calculate or correct for the topographic effects. For example, attempting to invert Equation 5 to solve for would involve the ratio 0∕ ; uncertainty in the denominator of a fraction often has significant consequences, especially if the denominator is small (Richter & Schläpfer, 2021, chapter 7). Table 3 shows the RMS errors for the cosine of illumination, along with the fraction of the terrain that is shadowed, for the dates in Figure 2 that extend from the winter to the summer solstice in equal changes of the solar declination. The RMS error for varies inversely with the value of 0 ; the errors in slope S and aspect A have a greater effect when 0 is smaller. The full extent of errors in the results indicates issues with outliers that the RMS errors do not reveal. Figure 6 shows scatter diagrams of calculated from the Copernicus DEM versus calculated from the Alaska IFSAR DEM. On all dates but particularly early in the year, some pixels that are illuminated ( ≫ 0) in the Copernicus DEM are in the dark ( < 0.1) in the Alaska IFSAR DEM. Similarly, some pixels that the Alaska IFSAR DEM shows to be illuminated are dark in the Copernicus DEM. A popular text on surveying published six decades ago (Davis et al., 1966) calls these kinds of mistakes "blunders" rather than "errors," because they cannot be characterized by an error distribution.

Discussion
Although errors or blunders in the NASADEM and Copernicus DEM are minor compared to the elevation values, their impact on remote sensing can be large. Thus, the small dispersion around the 1:1 line in the scatter diagram for elevation in Figure 5a translates to much greater dispersion in the slope, aspect, and view factor (Figures 5c-5e), which in turn translates to large dispersion in the illumination angles that Figure 6 shows. Therefore, small errors in slope or aspect can then significantly affect estimated reflectance, especially wherever is small.
Algorithms to retrieve surface properties differ in their sensitivities to topographic uncertainty. The effect is mostly a shift in spectral reflectance magnitude, so algorithms that rely on relative spectral shapes may escape significant harm. These include detection of materials based on diagnostic spectral absorptions, as in mineral identification (Clark et al., 2003). On the other hand, studies that rely on absolute radiometry, such as surface energy balance investigations (Wang et al., 2015) or retrieval of snow properties (Bohn et al., 2021), could be more severely affected. Moreover, errors in change the estimated balance between diffuse and direct illumination onto the surface. Therefore, they can distort the estimated reflectance spectrum in visible wavelengths, harming snow or vegetation studies that rely on features in this spectral range. Solar illumination geometry in mountains affects current satellite imagery from Landsat 8/9 and Sentinel-2A/B, it affects data from imaging spectrometers EnMAP (Chabrillat et al., 2020) and EMIT (Connelly et al., 2021), and it will affect data from future missions SBG (Cawse-Nicholson et al., 2021;Stavros et al., 2022) and CHIME (Rast et al., 2021). Locally, fine-resolution DEMs will be available from lidar, InSAR, or structure-from-motion deployed from drones or aircraft, and slightly coarser DEMs will be available using structure-from-motion from spaceborne data. However, the prospect is unlikely for globally available data to accurately estimate the solar illumination geometry for these imaging satellites. A chapter in Thriving on our Changing Planet (National Academies of Science, Engineering, & Medicine, 2018, p. 513) identifies applications that "would benefit from multibeam, space-based lidar to obtain global coverage of bare-earth topography and of the biomass/canopy at ≪5 m spatial and 0.1 m vertical resolutions." However, no such recommendation carried through to that report's Executive Summary, and no future NASA mission is in the planning stages.
Therefore, we face a future where the globally available DEMs at ∼30 m resolution are what we have now, at least through the launches and initial few years of the spectrometers SBG and CHIME and future versions of Landsat and Sentinel. If we could trust the variables calculated from DEMs and consider only the RMS errors, we could implement topographic correction algorithms that estimate from measurements of atmospherically corrected and thereby recover the geophysical and biological properties of the surface that govern spectral reflectance, with known uncertainty. However, we face the problem of outliers in the calculations of and less crucially Ω , so applying any correction algorithm globally on entire images would produce some incorrect, thus misleading, retrievals.
Strategies to mitigate the impact of topographic errors in processing and distributing image data and products must be considered. The list is deliberately terse; any bullet point could be expanded to a whole journal article: • In the basis documents for algorithms for geophysical and biological products, assess their sensitivity to uncertainty in illumination geometry and distinguish between topographic effects that change the spectral shape of the signal versus those that change the magnitude only (Lamare et al., 2020). • Gain a better understanding of the use of shade endmembers (Adams et al., 1986) in spectral mixture analysis, which implicitly acknowledge the limitations of available DEMs by solving for an illumination adjustment on modeled values of a pixel's reflectance. • Understand the relative magnitudes of topographic effects on angular properties of the reflectance versus the effects of illumination and viewing geometry on the intrinsic reflectance (Roupioz et al., 2014;Schaepman-Strub et al., 2006). • Develop and validate image processing methods that identify pixels where errors in the underlying DEM would lead to incorrect calculations of the illumination geometry, for example, detection of shadowed terrain (Hagolle et al., 2017;Hollstein et al., 2016;Shahtahmassebi et al., 2013). • Avoid exclusively prescribing global topographic correction solutions. Preserve the flexibility, within the mission science data system, for investigators to apply new regional DEMs of higher accuracy as these become available, or to ignore topography.
In the longer term, future research may reduce DEM-induced reflectance errors through strategies such as the following: • Implement topographic corrections in superpixels, thereby smoothing out the errors in individual pixels (Gilmore et al., 2011). • Continue efforts to improve DEMs globally, especially in mountainous areas, for example, the USGS 3D elevation program in the U.S. (Stoker & Miller, 2022).
• Examine and validate novel methods to estimate illumination geometry directly from images, for example, by simultaneously solving for unknown atmospheric and topographic properties in retrieval of surface reflectance from top-of-atmosphere radiances.

Conclusions
Our analyses show that calculations in the globally available DEMs miss shadows and consistently underestimate cosines of solar illumination angles, RMS error increasing with solar zenith angle. Analyzing imagery of Earth's mountains from current and future missions requires addressing the uncertainty introduced by errors and outliers in the DEMs on algorithms that retrieve surface properties from measurements of the apparent spectral reflectance. Intriguing potential improvements lie in assessing the uncertainties in retrievals of geophysical and biological properties and in novel methods to gain information about topography from the imagery itself.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
We have assembled all elevation data used in this research in Dryad  under the Creative Commons Zero Waiver. Those elevation files include splicing and cropping to match areas of fine and coarse resolution, so within each region (Carson River, Himachal Pradesh, Wrangell Mountains) the DEMs at different resolutions cover the same area, thereby enabling the comparison of the same topographic variables calculated from different data sources. Public sources of the data are: • NASADEM tiles are available from the U.S. Geological Survey Land Processes DAAC Data Pool (NASA JPL, 2020). Registration is required but is free. • Copernicus DEMs customized to specific latitude-longitude quadrilaterals are available from Open Topography (European Space Agency, 2021). • Airborne Snow Observatories Inc. provided the snow-off elevation data at 3 m spatial resolution for the Carson River Watershed. The data are available in Dozier et al. (2022). • The Alaska elevation data, acquired by airborne interferometric synthetic aperture radar, are available from the U.S. Geological Survey (USGS EROS Archive, 2018). • Tiles for the High Mountain Asia 8 m DEM are available at the National Snow and Ice Data Center (Shean, 2017). • Global grids of the EGM96 and EGM2008 Geoids are available from Agisoft (2008).
Computer codes for calculating solar illumination geometry (Dozier, 2020) and topographic horizons and other terrain parameters (Dozier, 2022a) are available from the MATLAB Central file exchange. Code for reprojecting raster data is on GitHub (Dozier, 2021). All codes are published and copyrighted under a free re-use license, even for commercial purposes.