Collected over an 11-day mission in February 2000, the SRTM DEM—derived using C-band radar interferometry—has led to significant advances in (near) global topographic analysis (Farr et al., 2007). The ∼3 global passes of the C-band returns were used to generate DEMs at resolutions of 1 (∼30 m) and 3 (∼90 m) arcsec. Since the initial data collection in February 2000, this dataset has seen many releases and void-filling efforts (e.g., Jarvis et al., 2008), with the most recent being the reprocessed 1 arcsec SRTM-NASADEM (Crippen et al., 2016; Buckley et al., 2020). Absolute vertical uncertainties of SRTM from reference datasets on bare-earth range from ∼2–10 m depending on the chosen statistical metrics and topographic steepness (e.g., Rodríguez et al., 2006; Becek, 2008; Becek, 2014; Rexer and Hirt, 2014; Purinton and Bookhagen, 2017).

The ASTER optical sensor aboard the Terra satellite collects nadir and backwards pointing imagery at a 15 m resolution since 1999. These stereo pairs are used to generate 30 m DEMs (e.g., Kääb, 2002), which are stacked to form the ASTER-GDEM (ASTER, 2009; Tachikawa et al., 2011). The most recent ASTER-GDEMv3 was generated by stacking over 2.3 million individual DEMs (Abrams and Crippen, 2019; Abrams et al., 2020) with a varying number of stacked scenes for each row and column of the satellite orbit, with an average number of 23 (σ = 18) for the study area. Generally high absolute vertical uncertainties of the ASTER DEMs have been reported ranging from ∼6–20 m, again depending heavily on topography (e.g., Tachikawa et al., 2011; Becek, 2014; Rexer and Hirt, 2014; Baade and Schmullius, 2016; Purinton and Bookhagen, 2017).

The ALOS Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM), launched in 2006, provides another optical source of DEMs (Tadono et al., 2014). With along-track nadir, backward, and forward viewing cameras at 2.5 m resolution, the images are used in a similar manner to ASTER to produce tri-stereo DEMs at a resolution of 5 m, of which over 3 million were stacked to form the World 3D 5 m DEM available commercially (Takaku et al., 2014). A 1 arcsec resampled version of this dataset is available for public access, with the most recent edited and hole-filled version being the ALOS-W3Dv3.1 used in this study (EORC, 2021). Over the study area, the number of stacked PRISM DEMs used to generate the final product is on average 7 (σ = 5). Absolute vertical uncertainties for the ALOS World 3D DEMs have been reported to a more limited extent, but may be expected to range from ∼2–10 m (e.g., Purinton and Bookhagen, 2017) depending on terrain slope.

The TanDEM-X 0.4 arcsec (∼12 m) DEM is the next generation of radar-derived global topography following the SRTM. This DEM was generated by interferometric processing and stacking of >470,000 X-band radar TerraSAR-X/TanDEM-X satellite bistatic scenes collected between 2010 and 2015 (Krieger et al., 2013; Wessel, 2016; Rizzoli et al., 2017; Zink et al., 2021). These data were also resampled, without further processing via different multi-looking, to a 1 and 3 arcsec version. The bistatic scenes have also been used to create the 10 m commercial AIRBUS WorldDEM^TM (Zink et al., 2021) with careful manual editing (e.g., void filling, water-body flattening, smoothing). The 3 arcsec (∼90 m) TanDEM-X is now available for public access, but the 1 arcsec version used in the present study was received through a scientific proposal to the DLR (Zink et al., 2021). The number of individual stacked scenes in the final product is on average 7 (σ = 4) for our study area. The TanDEM-X absolute vertical uncertainty is in the range of ∼1–5 m (e.g., Purinton and Bookhagen, 2017; Wessel et al., 2018), representing a significant improvement over previous near-global DEMs.

The Copernicus DEM, publically released at 3 arcsec in 2019 and 1 arcsec in 2021, is a derived product from the TanDEM-X mission. This dataset is also available at 10 m resolution over Europe, and was generated from the commercial WorldDEM^TM by the European Space Agency (ESA) and AIRBUS, including additional editing and smoothing of the 1 and 3 arcsec products (Leister-Taylor et al., 2020). The recently released 1 arcsec Copernicus DEM used here should therefore represent an improvement over the similar 1 arcsec TanDEM-X (scientific product generated by the DLR from resampling the 0.4 arcsec version without editing); however, thus far, limited validation reporting exists (Fahrland et al., 2020; Guth and Geoffroy, 2021). We also emphasize that the processing and filtering steps of the commercial TanDEM-X WorldDEM^TM product leading to the Copernicus DEM are not open-source documentation. The Copernicus DEM has been assessed with ICESat-2 measurements, which indicate absolute vertical uncertainties of ∼1–3 m (Fahrland et al., 2020), and previous work to assess the WorldDEM^TM indicates similarly high accuracy of this dataset (Becek et al., 2016).

Within the outlined context, we explore three options for inter-pixel consistency metrics. The different metrics can be seen in Figure 3, and details of each are provided below. All metrics are generated in small analysis windows (3 × 3 or 5 × 5 pixels, or in length ∼ 90 × 90 or ∼ 150 × 150 m), which emphasizes the variability in adjacent pixels. We view these metrics as a normalization procedure to better compare pixels to their local neighborhoods. This also can be described as a detrending or local filtering step. This is necessary, because analysis directly on the elevation pixels may mask the inter-pixel consistency signal under large-scale topographic signatures.

Because of local topographic variations, these metrics also reflect true local topographic signatures. For example, they result in high values in areas of high curvature (e.g., ridges and valleys). But at the same time, they also capture the signal of inter-pixel consistency, which has a very different spatial pattern compared to the discrete, localized ridge and valley signals of real topography. We emphasize the bare-earth and nearly vegetation-free signal of this area (Figures 1C,D): all metrics rely primarily on the assumption of field and topographic knowledge of the area of interest. If vegetation were present, we may expect a rougher surface from the spaceborne DEMs, and this will differ depending on the optical or radar collection method. Our method thus relies on some a priori knowledge of the expected surface in the study area, which can be gained through field knowledge and/or additional remote sensing data (e.g., rainfall, vegetation). But despite this, metrics can be compared between different DEMs to provide a relative assessment of inter-pixel consistency and its geomorphic implication. In the comparison step, we use the power-frequency distributions of the consistency metric to exploit the frequency, as opposed to spatial, domain. A pixel-by-pixel comparison of the various DEMs in the spatial domain would require careful sub-pixel co-registration steps (e.g., Nuth and Kääb, 2011).

Having selected an inter-pixel consistency quality metric that accentuates the variability in adjacent pixels, we seek to quantify it. For this, we turn to frequency analysis using the two-dimensional discrete Fourier transform (2D DFT) on the HPHS grids. This converts gridded values from the spatial to the frequency domain, which provides information about the amplitude and periodicity of the grid. In other words, the 2D DFT quantifies the variance at discrete wavelengths (units of distance, where frequency = wavelength⁻¹). This technique has been used in prior assessments of DEM quality and artifact removal (Arrell et al., 2008; Purinton and Bookhagen, 2017; Yamazaki et al., 2017), and this and other wavelet-based frequency analysis are increasingly popular for characteristic landscape scaling and feature identification (e.g., Perron et al., 2008; Booth et al., 2009; Roering et al., 2010; Hooshyar et al., 2021; Struble et al., 2021; Wapenhans et al., 2021).

With the longer-wavelength peaks and higher-frequency variance quantified, we seek to compare their effects on geomorphic analysis. Calculations are always done on the UTM projected grids (30 m resolution, resampled via cubic spline).

In a first step, we calculate the slope distribution for each of the three test catchments (Figure 1). We focus only on slope as it is the first derivative of elevation, and inaccuracies caused by low inter-pixel consistency in this metric will manifest in other DEM derivatives, such as curvature and aspect. Gridded slope in degrees is calculated in a 3 × 3 pixel window following the standard Zevenbergen and Thorne (1987) algorithm. The slope distributions for each DEM in each catchment are then compared and connected to the inter-pixel consistency quantified with the Fourier analysis.

We also extract flow networks for each catchment via the D8 algorithm (O’Callaghan and Mark, 1984) using each DEM and compare the longitudinal river profile of the longest channel (trunk stream). We assess the difference in normalized channel steepness (k _sn ): k s n = S A − θ r e f (7) where S is channel gradient, A is drainage area, and θ _ref is a reference channel concavity. This is a standard analysis in tectonic geomorphology used to identify channel knickpoints and their tectonic and climatic forcings (e.g., Wobus et al., 2006). Drainage extraction and k _sn calculations (using a θ _ref = 0.45 and minimum drainage area threshold of 0.1 km²) were done with LSDTopoTools (Mudd et al., 2019), which implements the integral approach to channel steepness (Perron and Royden, 2013), with the addition of segmented fitting (Mudd et al., 2014).

In a final step, we explore the effect of inter-pixel consistency on local channel gradient calculations (channel node to channel node along the profile). We follow the method of Clubb et al. (2019), and apply linear regressions to local windows of connected channel nodes to extract slope. We vary the window size of gradient calculation over 3, 5, 7, 9, and 11 nodes (∼90–330 m channel segments) to explore its effect on reducing artifact-related variance in the gradient.

The peaks quantified in Figure 7 correspond only to the ≥2-arcsec (2-pixel) wavelengths in the HPHS grid on the un-projected 1 arcsec DEMs. The sinusoidal DEM artifacts at longer wavelengths (Figure 7) are notable, but the shorter-wavelength variance has the primary impact on inter-pixel consistency in local windows used for many topographic calculations. The variance (percent of total power) corresponding to the 42- to <60-m wavelengths (i.e., adjacent x, y, and diagonal pixels) for the 30-m UTM resampled DEMs is quantified in 96 square (∼20-km) tiles in Figure 8. Here, we are able to use a smaller tile size to extract a larger distribution (more tiles compared to 60-km tiling), since we are not interested in the longer-wavelength periodic features, which are more subtle (lower total percent of power) and only become prominent when considering larger areas in the frequency analysis.

Figure 8 only presents the results for the nearest neighbor, bilinear, and cubic spline resampling schemes, while Supplementary Figure S2 shows the average, cubic, and lanczos resampling methods. The percent of power at 42- to <60-m wavelengths can be interpreted as the relative inter-pixel consistency between the five DEMs. Lower values (median and interquartile range <7.5%) for the Copernicus and TanDEM-X correspond to relatively higher inter-pixel consistency (smoother appearance in Figure 2A), compared to the lower inter-pixel consistency implied by the higher percent of power values (median and interquartile range >7.5%, often >10%) for ALOS-W3Dv3.1, ASTER-GDEMv3, and SRTM-NASADEM. We note that the TanDEM-X has somewhat higher percent of power at < 60-m wavelengths compared to the Copernicus DEM (derived from the same source data as the TanDEM-X research product), which is a result of the careful editing and smoothing of the Copernicus DEM. This is difficult to discern from the plot, but, for example, the 25th percentile of the Copernicus boxplot for the cubic spline resampling is at ∼2.49%, and at ∼2.51% for the TanDEM-X. In most cases, variance is decreased (inter-pixel consistency is increased) going from nearest neighbor to bilinear to cubic spline resampling, with the SRTM-NASADEM being the exception. The variance is also decreased going from lanczos to cubic to average resampling in Supplementary Figure S2, but the cubic spline resampling in Figure 8 produces the greatest reduction in adjacent pixel variance.

As previously stated, the use of an inter-pixel consistency metric for comparing DEMs relies on some assumptions. First and foremost, the characteristics of the study area must be known via field and topographic knowledge or observation in satellite imagery or other remotely sensed datasets (e.g., vegetation, rainfall). The presence of vegetation or snow and ice would modify the performance of a given DEM. For instance, depending on the radar wavelength (e.g., C-band for SRTM and X-band for TanDEM-X), there will be different penetration depths of canopy (e.g., Carabajal and Harding, 2006; Hofton et al., 2006; Wessel et al., 2018) and snow (e.g., Rignot et al., 2001; Rossi et al., 2016), and for optical DEMs (e.g., ASTER and ALOS) only the top surface is recorded. Our study area presents an opportunity to examine inter-pixel consistency under ideal conditions: arid, nearly vegetation-free, and mixed steep (mountain ranges) and flat (salars) topography (Figures 1C,D). Field observations and field measurements of our study area show very low inter-pixel variations at 30 m spacing, and we can interpret deviations in the consistency metrics to be DEM artifacts and not topographic signal. We expect that the observed artifacts will be present in other regions, and these will be further impacted by land-surface characteristics, which may mask or amplify the artifacts in complex ways.

A second assumption was used for objective peak identification. Here, we assumed that the Copernicus DEM was free of longer-wavelength artifacts and used this to take a normalizing ratio to highlight anomalous peaks in spectral power in the other datasets. This step is justified by the high-quality source of the carefully edited AIRBUS WorldDEM^TM—based on the original TanDEM-X data. This said, the ratio step is only necessary for automatic peak identification. In many cases, the wavelength peaks could be graphically identified given their prominence against background values (high power, e.g., Figure 5), and the method is useful as a stand-alone (without reference data or another DEM) approach for the identification of ≥2-pixel wavelength artifacts. In a conservative sense, the automatic identification of repetitive signals on the DFT spectrum presented in this study is only relative to the Copernicus DEM. However, the second step to quantify the variability of high-frequency (<2 pixel) adjacent pixels in UTM 30-m reprojected DEMs does not rely on a reference surface, and instead acts as a comparative metric between DEMs, under the first assumption of field characteristics: are smoother or rougher local topographic surfaces expected? This also provides quantitative assessment of differences in resampling, which is a key first step in topographic analysis.

We note that the DFT can only be carried out on void-free tiles, and all DEMs used here are void-filled versions. This is accomplished by a combination of interpolation (e.g., Copernicus voids ≤16 pixels in size interpolated from surrounding terrain; Leister-Taylor et al., 2020) and/or replacement with other newer or older DEMs (e.g., SRTM-NASADEM voids filled by ASTER-GDEM and ALOS-W3D; Buckley et al., 2020). We do not expect the void-filling to alter the consistency metrics reported here, as these are discrete areas making up a small portion of the DEMs, particularly in our arid study area, where atmospheric and land-cover challenges to spaceborne DEM collection and processing are minimized.

Fourier analysis performed directly on elevation grids are taken to characterize landscape scales (Perron et al., 2008; Hooshyar et al., 2021), identify and quantify pit-and-mound (Roering et al., 2010) or landslide (Booth et al., 2009) topography, and even identify and remove striping artifacts in lidar (Arrell et al., 2008) and SRTM DEMs (Yamazaki et al., 2017; Purinton and Bookhagen, 2018). Akin to these studies, we previously calculated the 2D DFT directly on the gridded elevation values to highlight artifacts in 2–8 pixel wavelengths for an 8 × 14-km clip of the ALOS-W3D 5 m commercial DEM in Purinton and Bookhagen (2017). This may be appropriate for higher-resolution DEMs (<10 m) in small study areas, but the large-scale (60-km tiles) and diversity of topography in the 1 arcsec DEMs used here necessitates a neighborhood filtering approach.

In Figure 9, we present a normalized 2D power spectrum calculated directly on the elevation values for the same SRTM-NASADEM tile shown in Figure 4. The influence of longer (i.e., several dozen to hundred pixel steps) topographic wavelengths of valleys, ridges, mountain chains, and salt flats reduces and/or obscures the inter-pixel consistency, and the periodic signals of artifacts. By locally detrending the topography using a metric like HPHS, the wavelengths of nearby pixel variability are highlighted, preparing the DEM data for a more meaningful inter-pixel consistency analysis. Furthermore, this and other metrics (internal smoothing dR, plane-fit RMSE), provide an additional visual assessment of DEM quality and the spatial patterns of noise, closely tied to hillshade observations (Figure 2).

In another step, we tested the result of increasing the kernel size for high-pass filtering in the HPHS calculation (Supplementary Figure S1). Our 3×3-pixel filter particularly accentuates the variability of adjacent pixels—and notably also accentuates longer-wavelength periodic patterns—whereas increasing the high-pass kernel size to 5 × 5 pixels demonstrates that while the longer-wavelength patterns are still visible in the 2D DFT, the high-frequency pixel-to-pixel variance is reduced (Supplementary Figure S1).

A recent review by Polidori and El Hage (2020) highlights the need for different approaches of inter-pixel consistency reporting to improve understanding of DEM quality beyond vertical accuracy. This is key given increases in quantitative geomorphometry and dissemination of DEMs from many (sometimes poorly understood) sources (Sofia, 2020). The metrics developed here provide a more suitable assessment of DEM quality compared to point-based vertical accuracy, which does not account for the spatial variability of DEM vertical errors that impact derivatives of elevation. These metrics do not use reference surfaces (e.g., Kramm and Hoffmeister, 2019) to assess spatially continuous patterns of vertical error, which requires co-registration of the surfaces. Co-registration of the DEMs would be particularly important for the case of the ALOS-W3Dv3.1, which has the pixel edges aligned to integer coordinates in latitude and longitude (EORC, 2021), rather than the pixel centers as in all other DEMs, leading to a half pixel offset. Our Fourier steps transform the analysis from the spatial to the frequency domain, thus negating this co-registration requirement (Nuth and Kääb, 2011), which may introduce additional uncertainties from model fitting and/or resampling.

The Copernicus and TanDEM-X DEMs have the highest inter-pixel consistency (Figure 8), with a generally smooth representation of the arid landscape in our study area. The TanDEM-X DEM is a research-grade product and still contains some noise over water bodies and on steep slopes (particularly eastern and western facing slopes facing the TerraSAR-X/TanDEM-X satellite look direction). This manifests in slightly higher <60-m wavelength variance for the 30-m projected tiles, but the TanDEM-X does not have any longer-wavelength (>2 arcsec) anomalies in the 1 arcsec tiles.

The long-wavelength artifacts in the SRTM-NASADEM have previously been identified by many authors (e.g., Gallant and Read, 2009; Yamazaki et al., 2017; Purinton and Bookhagen, 2018; Grohmann, 2018). Here, we build on this previous work and develop another approach to quantify the wavelengths and orientations of these artifacts. The mast-oscillations during collection are responsible for the original artifacts (Farr et al., 2007), and this is clear from the north-northeast and north-northwest orientation of the waves along the ascending and descending passes of the shuttle mission. However, the multiple wavelengths (2, 12, and 23–24 arcsec; Figure 7) indicate possible interference-related harmonics from the mast oscillations, ∼3 stacked passes, bilinear resampling steps, and/or attempts to remove the ripples using ICESat measurements while producing the recent SRTM-NASADEM (Buckley et al., 2020). These artifacts could be removed via band-pass filtering on the selected wavelengths identified as spectral peaks (Yamazaki et al., 2017; Purinton and Bookhagen, 2018), but these peaks should be calculated for each 1° × 1° DEM tile separately prior to filtering as they may not be globally (or even regionally) consistent.

Aside from the SRTM-NASADEM (only ∼3 passes collected in 11 days), the other datasets represent multi-year collection efforts with hundreds-of-thousands (TanDEM-X and Copernicus) to millions (ASTER and ALOS) of individual stacked DEMs. All five DEMs are delivered with auxiliary rasters containing pixel-level information on their source. For instance, the coverage (COV) raster for TanDEM-X with the number of TerraSAR-X/TanDEM-X scenes (Wessel, 2016), and the mask (MSK) raster for ALOS with the number of PRISM scenes or filling source (EORC, 2021). Increasing the number of stacked scenes, especially in complex topography, can improve DEM quality (e.g., Purinton and Bookhagen, 2017), but inherent errors and biases will continue to limit quality. In any case, these void- and stack-masks, along with other auxiliary files, can be a useful check on DEM quality. Height error maps (HEM) delivered with the Copernicus, TanDEM-X, and more recently SRTM-NASADEM DEMs (Wessel, 2016; Buckley et al., 2020; Leister-Taylor et al., 2020) can be useful for gaining a first impression of accuracy, but these are based on interferometric coherence, which may experience decorrelation due to a number of surface conditions (e.g., vegetation cover, moisture). Our goal here was to consider only the elevation surface to extract internal error metrics without any other reference data or any background DEM processing data, which is how many end users receive and utilize the gridded DEMs.

Notably for the TanDEM-X, there has been recent efforts by González et al. (2020) to improve the quality via careful editing and smoothing (particularly over water bodies) of the 30 m version, which may bring this DEM closer to the Copernicus DEM. Additionally, ongoing bistatic scene collection to generate scientific research products, such as a new global DEM collected from 2017 to 2020, could lead to further improvements in the TanDEM-X DEM (Zink et al., 2021). We do note that observations in the study area show a smoother surface in the Copernicus DEM compared with TanDEM-X; however, in local areas this smoothing has noticeably flattened true topographic expression in the form of rough bedrock outcrops. This is an inevitable result of automatic DEM editing, which often has significant moving-window smoothing steps. Other recent efforts towards DEM editing and fusion (e.g., Yamazaki et al., 2017; Liu et al., 2021) may create new and improved DEM products, but these steps must be carried out carefully and the underlying DEM data (typically from SRTM, ASTER, ALOS, and/or TanDEM-X) can still propagate uncertainties into the final product. Any DEM created by multi-step editing of spaceborne data should be treated with care and analyzed for inter-pixel consistency, in addition to traditional vertical accuracy metrics.

High power in adjacent pixel steps measured on the HPHS grid (Figure 8) are a sign of low inter-pixel consistency in our geomorphically smooth, nearly vegetation-free study area. For the SRTM-NASADEM, the longer-wavelength errors are orbital- and, possibly, processing-related errors. The high-power in adjacent pixels can be accounted for by sensor errors (signal to noise ratio) and speckle associated with radar interferometric generation (Buckley et al., 2020). This radar speckle can be improved by multilooking, as in the case of TanDEM-X (Wessel, 2016; Rizzoli et al., 2017), but for the SRTM-NASADEM, the 30-m nominal resolution is likely already beyond the true sampling resolution of this sensor, reported as potentially 45–60 m (Sun et al., 2003; Farr et al., 2007; Tachikawa et al., 2011).The optical DEMs used here (ASTER-GDEMv3 and ALOS-W3Dv3.1) appear to suffer primarily from processing artifacts. In particular, the ALOS World3D 5 m DEM (underlying data for the 30-m product) and ASTER DEMs both experience short-wavelength artifacts (low inter-pixel consistency of adjacent pixels and ≤210-m wavelength artifacts without a consistent orientation for ASTER) likely related to the size of correlation kernels used in photogrammetric reconstruction and errors with tie point matching between the image pairs (ASTER) or triplets (ALOS). The ALOS-W3Dv3.1 was generated by average resampling of the original 5-m pixels (EORC, 2021), and this manifests in a distinct artifact with power peaks aligned orthogonal to the grid (due north and due east, Figure 7A) in 3-pixel steps. This again highlights the benefit of the HPHS calculation and Fourier analysis to identify resampling artifacts, which may vary for different resampling schemes.

Resampling during reprojection is a common pre-processing step prior to any topographic analysis and is a requirement by such software as TopoToolbox (Schwanghart and Scherler, 2014) and LSDTopoTools (Mudd et al., 2019). Thus, the differences in high-frequency (adjacent pixel) variance for different resampling schemes are particularly notable results of the analysis (Figure 8; Supplementary Figure S2). Often resampling is done using bilinear or nearest-neighbor approaches, as these are quick and thought to provide reasonable results. Resampling methods using larger or adjustable window sizes or higher-order polynomials, such as cubic spline resampling reduces the variance in adjacent pixels and increases inter-pixel consistency. However, cubic spline resampling may also be over-smoothing locations with high topographic variability at short distances. This is a natural result of DEM smoothing and points to a nuance of 30-m DEM usage: the topographic variability at short distances is convolved with the signal of non-topographic variance. Smoothing these DEMs to increase inter-pixel consistency while retaining topographic signatures will require adaptive resampling schemes. In any case, going forward with the geomorphic analysis we use the cubic spline resampling. Importantly, we note that different resampling schemes only change the magnitude and not pattern (relative difference between DEMs) of the results.

Moving from a tile-based approach to quantify the variance in each DEM at specific wavelengths (pixel steps), we turn to our selected catchments (Figure 1) to assess the impact of inter-pixel consistency differences for geomorphic research.

One interesting consideration from the analysis concerns the ALOS-W3Dv3.1. In Figure 8, the inter-pixel consistency is lower than for the Copernicus (higher percentage of power at 42- to <60-m wavelengths). However, the slope distributions are very similar to the Copernicus DEM (Figure 11B), compared to the ASTER-GDEMv3 and SRTM-NASADEM, which have much lower inter-pixel consistency (>∼10% in Figure 8). Notably, despite the lower inter-pixel consistency implied by the boxplot in Figure 8, the ALOS-W3Dv3.1 has the lowest variability in channel gradient of all DEMs (Figure 13A). This is a possible sign of a hydrologic preconditioning step of this dataset, though unreported in the technical documentation (EORC, 2021).

Although longitudinal river profile analysis at large (several km) scales in steep, high-relief catchments for tectonic and climatic forcings are less affected by DEM choice (e.g., k _sn ), other catchment analysis like slope distributions will be impacted by this choice. Furthermore, fine-scale analysis of channel gradients will be particularly impacted by the channel-node to channel-node variability, with implications for assessing reach-scale channel morphology and knickpoint detection for tectonic geomorphology (e.g., Neely et al., 2017; Clubb et al., 2019; Gailleton et al., 2019) and river ecology (e.g., Beechie and Sibley, 1997; Bisson et al., 2017), among others.

This study does not address DEMs of different resolution, as we focus only on the widely used and open-source (near) global 30 m DEMs. We do note that DEM resolution will impact geomorphic analysis and the calculation of derivatives of elevation (Purinton and Bookhagen, 2017; Grieve et al., 2016; Smith et al., 2019), but many studies in remote or extensive areas will continue to rely on the 30-m datasets. Although the spatial resolution of these DEMs is too coarse for channel-head detection (e.g., Passalacqua et al., 2010a,b; Clubb et al., 2014; Hooshyar et al., 2016), high inter-pixel consistency, especially in low slopes near drainage divides, will allow more accurate flow path derivation using divergent flow routing algorithms like D∞ (Tarboton, 2005). This is an important component in soil mantled and diffusional landscapes, irrespective of vegetation cover.

Previous work to quantify vertical error in DEMs (e.g., Becek, 2008; Becek, 2014; Becek et al., 2016; Yamazaki et al., 2017) highlights that error can be introduced by various sensor biases, quantization of elevation values (i.e., integer versus floating point), and land cover. In our study, a primary goal is to avoid reference data and develop an internal metric to compare DEM quality using only a priori knowledge of a smooth, bare-earth study area with mixed steep and flat terrain. While the methods developed by Becek (2008) demonstrate a novel use of airplane runways as continuous reference surfaces, this is limited to flat terrain of limited spatial extent and requires some assumptions, such as a uniform distribution of quantization and slope-induce errors. In any case, the results from the runway method on SRTM (Becek, 2008), ASTER GDEM (Becek, 2014), and WorldDEM^TM (Becek et al., 2016) support the relative quality between these datasets found in this study (where we use the Copernicus DEM based on the WorldDEM^TM).

Technical documentation available for spaceborne DEMs (Wessel, 2016; Abrams and Crippen, 2019; Buckley et al., 2020; Leister-Taylor et al., 2020; EORC, 2021; Zink et al., 2021) contain important information on processing and limitations, but this information is often neglected and end users of 30 m DEMs require fundamental metrics of DEM quality beyond vertical point-based accuracy. Our local pixel variability metric based directly on the hillshade image (which includes slope and aspect information), combined with quantification of inter-pixel consistency and demonstration of the impacts on geomorphic research, provides a quantitative explanation of spaceborne DEM quality.

From our analysis, it is clear that the newly released Copernicus DEM—based on the high-quality TanDEM-X derived WorldDEM^TM—has the highest inter-pixel consistency, and therefore most realistic representation of the topography in our study area. This is in agreement with the findings of other studies comparing TanDEM-X to previous DEMs (Pipaud et al., 2015; Purinton and Bookhagen, 2017; Boulton and Stokes, 2018; Grohmann, 2018), and recent reporting on Copernicus (Guth and Geoffroy, 2021). We do note that other authors have found better performance of the ALOS-W3D for stream profile analysis (Schwanghart and Scherler, 2017), possibly due to hydrologic preconditioning. Older DEMs like the SRTM-NASADEM and ASTER-GDEMv3 continue to be developed, and these can be useful as time-shots of the Earth surface during collection, but geomorphic analysis at the catchment and mountain belt scale should increasingly rely on the newer Copernicus and TanDEM-X global DEMs.