How does grid-resolution modulate the topographic expression of geomorphic processes?

Abstract

In many locations, our ability to study the processes which shape the Earth are greatly enhanced through the use of high-resolution digital topographic data. However, although the availability of such datasets has markedly increased in recent years, many locations of significant geomorphic interest still do not have high-resolution topographic data available. Here, we aim to constrain how well we can understand surface processes through topographic analysis performed on lower-resolution data. We generate digital elevation models from point clouds at a range of grid resolutions from 1 to 30 m, which covers the range of widely used data resolutions available globally, at three locations in the United States. Using these data, the relationship between curvature and grid resolution is explored, alongside the estimation of the hillslope sediment transport coefficient ($D$, in m${}^{\mathrm{2}}$ yr${}^{-\mathrm{1}}$) for each landscape. Curvature, and consequently $D$, values are shown to be generally insensitive to grid resolution, particularly in landscapes with broad hilltops and valleys. Curvature distributions, however, become increasingly condensed around the mean, and theoretical considerations suggest caution should be used when extracting curvature from landscapes with sharp ridges. The sensitivity of curvature and topographic gradient to grid resolution are also explored through analysis of one-dimensional approximations of curvature and gradient, providing a theoretical basis for the results generated using two-dimensional topographic data. Two methods of extracting channels from topographic data are tested. A geometric method of channel extraction that finds channels by detecting threshold values of planform curvature is shown to perform well at resolutions up to 30 m in all three landscapes. The landscape parameters of hillslope length and relief are both successfully extracted at the same range of resolutions. These parameters can be used to detect landscape transience and our results suggest that such work need not be confined to high-resolution topographic data. A synthesis of the results presented in this work indicates that although high-resolution (e.g., 1 m) topographic data do yield exciting possibilities for geomorphic research, many key parameters can be understood in lower-resolution data, given careful consideration of how analyses are performed.

Introduction

Geomorphologists have always made use of topographic data, from initial qualitative observations of surface morphology and its link to process e.g., to directly measuring landscape geometries from contour maps, constraining river dynamics and morphometric relationships (e.g., , , and ). Further quantitative analyses of the Earth's surface were facilitated through the advent of gridded topographic data. Work to generate digital elevation models (DEMs) from photogrammetry, contour maps, and active remote sensing platforms produced datasets at tens to thousands of meters' grid resolution, along with geomorphic analyses designed for such datasets . Algorithms have subsequently been developed which exploit the higher-resolution topographic data now available, predominantly from light detection and ranging (lidar), which not only refined existing techniques but also allowed the study of hitherto unresolvable features on landscapes .

Presently, lidar data coverage is predominantly focused around locations of particular scientific interest or infrastructural importance, as can be seen on many lidar data portals e.g.,. It is unlikely that global lidar coverage can be achieved in the near future, leaving the provision of commercially available 12 m TanDEM-X data and freely available 30 m Shuttle Radar Topography Mission (SRTM) data as the best available data options for many study sites.

As a consequence of this data availability it is crucial to understand the limitations of lower-resolution data when performing topographic analysis for geomorphic research. Extracting channels from topography is a common requirement of many analyses, and it is expected that the accuracy of extracted channel networks will be affected by increasing grid resolution . , , and used measurements of hillslope length and relief to identify signals of landscape transience. However, all such work was performed on high-resolution topography and the impact of grid resolution on these metrics is unknown. and demonstrated that the curvature of ridgelines measured from high-resolution topography can be used as a proxy for erosion rates in soil-mantled landscapes. This observation has been used in many studies in which cosmogenic radionuclide-derived erosion rates are unavailable . However, it can also be used in locations with an independent constraint on erosion rates in order to quantify a sediment transport coefficient that relates hillslope sediment flux to the topographic gradient, which is set by the material properties of soils . Therefore, understanding the effect of grid resolution on the extraction of curvature is crucial in order to evaluate the applicability of calculating the sediment transport coefficient from coarse-resolution data.

Here, we grid topographic data at a range of resolutions in order to test the sensitivity of these techniques to decreasing grid resolution, with the aim of placing constraints on the estimation of common geomorphic parameters when lidar topographic data are unavailable. Through an analysis of one-dimensional curvature and topographic gradient approximations, the changes in fidelity as grid resolution decreases for both curvature and topographic gradient are examined and placed within the context of the two-dimensional results of this study and the wider literature.

Previous work

It has long been recognized that the scale of topographic data used in an analysis or model will have an impact on the scale of the processes which can be measured . It is intuitive that in order to measure the properties of hillslope processes the resolution of the data must be high enough that variations in hillslope form can be captured adequately. The resolution of topographic data defines the Nyquist frequency, given as $(\mathrm{2}\text{Res}{)}^{-\mathrm{1}}$ where “Res” is the grid resolution of the dataset . The inverse of this frequency yields the minimum wavelength resolvable from a given dataset. In the example of a 1 m grid resolution, the smallest features that could be resolved would have a length scale of 2 m. Recognizing this, many authors have attempted to quantify this uncertainty, aiming to answer the following question: at what point does a dataset become unsuitable for a given analysis? e.g.,.

Many attempts to constrain the error content of topographic measurements have focused on comparisons between elevation values taken from differing resolution data products, often in conjunction with field survey data, with the aim of discriminating between DEM generation methods. performed a comparison of DEMs generated using cartometric and photogrammetric methods against field-surveyed elevation data. They demonstrated that at grid resolutions of 6.25, 12.5, and 25 m the cartometric DEM produced less error than the photogrammetric DEM when compared to the field-surveyed data, collected at 3.25 m intervals.

The advent of lidar-derived topographic data provided a new technique and increased the range of possible grid resolutions to evaluate. assessed the quality of high-resolution topographic data sourced from interferometry and lidar for a heavily vegetated catchment in North Carolina. This analysis demonstrated that, under such conditions, the lidar-derived DEM outperformed the interferometric data in addition to both classes of USGS DEM product. However, concerns were raised about the overall accuracy of the lidar data with a requirement for improved methodologies to be developed to process multistory vegetation. Further work was carried out in North Carolina to constrain the minimum number of lidar returns required to generate a DEM at a given grid resolution . This work indicated that a 5 m grid (the finest resolution used) required approximately 115 points ha${}^{-\mathrm{1}}$, whereas at 30 m grid resolution the requirement reduced to approximately 35 points ha${}^{-\mathrm{1}}$.

resampled a 1 m lidar-derived DEM to a range of grid resolutions up to 25 m and assessed the accuracy of elevation values for each of these resampled grids when compared to a 1 m resolution field survey. It was found that there was little variation in the distribution of elevation values between the resampled data sets. However, when the data was compared with 25 m DEMs generated from topographic maps and contour generalization, there were considerable errors, supporting earlier authors' conclusions that lidar-derived topographic data contain more useful geomorphic information than other methods of topographic data collection.

Topographic gradient (or slope) is one of the most fundamental topographic derivatives across the disparate disciplines which utilize topographic data. This measurement has been used in geomorphology e.g.,, ecology e.g.,, soil science e.g.,, and hydrology e.g.,. endeavored to constrain the accuracy with which this parameter can be calculated as grid resolution is increased from 100 to 1000 m and showed that as the grid resolution is decreased, there is a clear reduction in the slope values produced for a landscape. Similar wide-scale analysis has also been performed within the context of global hydrological analysis (e.g., , and ), indicating that from meter to kilometer scale the reduction in quality of slope measurements is an issue which must be considered when working with topographic data.

considered the accuracy of slope measurements at locations manually classified as valleys, peaks, and ridges. They found an initially small increase in the error of slope measurements at intermediate resolutions (10–20 m) and a much more rapid increase in error between 20 and 30 m resolution, suggesting a threshold minimum resolution for analysis of these landforms. More recent work has considered how high-resolution lidar data impact the quality of slope measurements. demonstrated a similar trend to previous authors working with lower-resolution data: as grid resolution is decreased from 1 to 25 m, there is a considerable reduction in the slope values generated for a landscape. evaluated the reliability of slope measurements by contrasting 10 methods of gradient calculation against field measurements of topographic gradient. The error between DEM and field-derived slope measurements was shown to increase with decreasing grid resolution (from 1 to 12 m), resulting in the recommendation to increase data resolution wherever possible to decrease errors in topographic analysis.

Numerous authors have considered the impact of grid resolution on hydrological applications, which often require slope calculation as a fundamental processing step. It has been demonstrated across many landscapes and scales that as grid resolution is decreased the upslope contributing area will increase and the local slope will decrease, which will have a significant impact on any hydrological analysis . Similarly, from the perspective of modeling global-scale sediment fluxes to the oceans, noted that measurements of slope dropped logarithmically with increasing grid resolution, and failing to account for this may lead to a substantial underestimate of the contribution of steep, montane regions.

performed analyses on the accuracy of hydrological networks generated through photogrammetry and radar interferometry at 5 and 30 m grid resolution, respectively. Their error analysis was extended to consider the vertical errors generated both through the downsampling of the topographic data, as well as from the techniques used to capture the topographic information. Predicted catchment runoff was up to 7 % larger in the lower-resolution datasets, considered to be driven by both the vertical errors and the reduction in spatial resolution increasing variables such as upslope drainage area.

Topographic wetness index (TWI), calculated as $\ln (A/S)$, where $A$ is the specific upslope area and $S$ is the slope, is used as a single variable to compare the hydrological setting of differing parts of the landscape, providing insight into factors including groundwater properties and overland flow rates. used lidar data to test the robustness of TWI calculations on spatial scales ranging from 5 to 50 m, concluding that the most sensitive part of the TWI calculation was the specific upslope area measurements. This sensitivity resulted in significant variation in the TWI values across the range of resolutions tested. Predicted slope stability, modeled in part as a function of TWI, was assessed by , who demonstrated that, for large-scale landsliding, a lidar-derived DEM downsampled to 10 m resolution was more suitable to identify landslide hazard than the highest-resolution data available. This highlights the requirement to consider the scale of the process being studied when selecting the appropriate grid resolution for a study and corresponds to the challenges of selecting the correct size of smoothing window to capture processes on a suitable scale (e.g., , , and ).

The accuracy of channel network extraction from topographic data was tested by , who tested a 1 m lidar DEM and a 10 m photogrammetrically generated DEM against a field-mapped channel network in a catchment in Alberta, Canada. The 1 m lidar-derived channel network was found to be the best representation of the field-mapped channel network, exceeding the quality of an additional channel network mapped by hand from aerial photographs. However, as no intermediate datasets were tested, it is not possible to understand at what resolution the degradation in channel network extraction quality occurs for this location.

As models of agricultural soil loss depend heavily on topographic variables such as slope, work has been carried out to understand the influence of grid resolution on calculated rates of soil loss. tested data resolutions from 1 to 81 m and demonstrated that in all cases, rates of predicted soil loss increased with grid resolution. However, the rates of soil loss were also influenced by the type of flow routing utilized, with the multiple flow direction algorithm (e.g., , and ) proving most sensitive to resolution decreases. Work by considering models of crop yields in Colorado, USA, demonstrated that on relatively flat surfaces, such as agricultural fields, the spatial resolution is less important than the vertical accuracy when predicting crop yields, with significant errors being produced due to centimeter-scale vertical displacements. Decreasing the grid resolution from 5 to 30 m had limited effect on the yield calculations.

Although considerable work has been carried out on the sensitivity of various factors to grid resolution, much of it has been focused on a specific application (e.g., , , , and ) with few studies considering the impact of DEM grid resolution within a geomorphic context. Here we aim to extend existing methodologies to constrain the utility of low-resolution data products across a suite of geomorphic analyses to understand the following: (1) how hillslope length, topographic curvature, and relief vary with grid resolution; (2) how best to extract channel networks in lower-resolution datasets in order to minimize errors; and (3) whether it is possible to estimate sediment transport coefficients from low-resolution topographic data, where an independent constraint on erosion rate is available.

Theory and methods

Generating topographic data

Previous studies that have explored the impact of changing grid resolution on topographic or geomorphic parameters have typically produced coarser-resolution topographic data by downsampling the highest-resolution data product available for their study sites (e.g., , , , and ). Work has been undertaken to understand the influence of various re-gridding schemes on topographic measurements , with focus placed upon understanding the use of downsampling high-resolution data in order to facilitate computationally expensive analysis on larger spatial areas with minimal loss in data fidelity. However, as computational power increases, cost decreases and more efficient algorithms are developed , the need to downsample data for computational convenience becomes reduced. Instead, it becomes more important to understand the limitations of available data products, to facilitate geomorphic analysis in locations in which high-resolution topographic data are not available. This is of particular importance in many studies of natural hazards (e.g., , and ) in which data quality is limited. It will also open geomorphic research up to communities which do not have the resources to acquire high-resolution topographic data.

Lidar point cloud metadata.

Location	Point density	Vertical	Horizontal
	(points m${}^{-\mathrm{2}}$)	accuracy	accuracy
		(m)	(m)
Santa Cruz Island	$8.27$	${0.067}^{*}$	${1.07}^{*}$
Gabilan Mesa	$5.56$	$0.20\pm 0.15$	$0.11$
Oregon Coast Range	$6.55$	$0.07\pm 0.03$	$0.06$

${}^{*}$ Accuracy is the 95 % confidence level of the root mean squared error of measurements compared to static GPS control points.

As a consequence of these constraints we have generated topographic data for our three study sites without downsampling or re-gridding high-resolution data products, as is commonly performed . Instead we have followed established techniques to grid the processed lidar point cloud data provided by OpenTopography (http://www.OpenTopography.org) at a range of data resolutions which span from 1 m, considered to be the limit of the Oregon Coast Range dataset by to 30 m, which is equal to the grid resolution of the global SRTM dataset and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) dataset and in excess of the TanDEM-X dataset and as such should span the vast majority of grid resolutions used in modern geomorphic research. The direct comparison between elevation products generated using differing methodologies is challenging e.g.,, and more work is required within the context of geomorphic research to understand limitations in topographic datasets, such as SRTM and TanDEM-X, which arise from data capture and processing rather than purely from resolution constraints. By generating the topographic data from the same source, we aim to isolate the signal of decreasing data resolution, without the introduction of new sources of error which may arise from data collected using a different instrument. The error estimates of the raw point clouds used in this re-gridding process are provided by OpenTopography and can be found in Table .

Example shaded reliefs of the same section of Santa Cruz Island at increasing grid resolutions. All coordinates are in UTM Zone 11${}^{\circ }$ N. Panels (a)–(f) represent resolutions of 1, 2, 5, 10, 20, and 30 m. Tick spacing is in meters. The red box outlines an extensively gullied first-order drainage, clearly visible in the highest-resolution data, but as the grid resolution is decreased, this feature and its internal structure become indistinguishable from the surrounding hillslopes.

[Figure omitted. See PDF]

The point clouds are gridded using Points2Grid, which employs a local binning algorithm, searching for points within a circular window of radius defined by as $$\text{Radius}=\lceil \sqrt{\mathrm{2}}\text{Res}\rceil .$$

An inverse distance-weighted averaging approach is then performed to assign an elevation value to each grid cell. This approach, which has been employed in previous studies , yields a reliable representation of the topographic surface, with few data gaps and a minimal amount of interpolation. The level of interpolation performed is controlled by the density of lidar ground returns within each search window, consequently more interpolation may be performed in areas of high vegetation density such as the Oregon Coast Range. This is an additional source of error which must be considered when processing lidar data, and this consideration informed the selection of 1 m as the maximum resolution used in this study as it is the highest resolution these datasets can have been gridded to in the past (e.g., , and ).

The topographic data used in this study have been gridded at 20 resolutions, and Fig.  provides representative hillshades of a section of Santa Cruz Island, highlighting the degradation of topographic information as grid resolution is decreased.

Measuring curvature from topography

Landscape curvature has long been recognized as a key geomorphic characteristic of landscapes, from Gilbert's () qualitative observations of hilltop convexity to more recent approaches to quantify landform curvature using digital topography (e.g., , and ). However, unlike other key landscape properties such as gradient , hydrology , or soil characteristics , the influence of grid resolution on curvature has not been fully explored, particularly within a geomorphic context.

This is particularly important with the proliferation of high-resolution topographic data from lidar, allowing the analysis of curvature on increasingly fine scales. Recent developments in channel extraction techniques typically require the identification of topographic convergence in high-resolution topography using a curvature threshold. and demonstrated that hilltop curvature scales with erosion rate and as such demonstrated the importance of accurately constraining the impact of grid resolution on this landscape parameter. Its importance is highlighted by an increasing number of studies using this relationship as a proxy for erosion rate . Hilltop curvature can also be used to constrain the sediment transport coefficient of a landscape where an independent constraint on erosion rate is available .

The measured curvature of a topographic surface depends on the orientation of the measurement. Here, we consider two common types of curvature, with the following definitions: (1) total curvature ($C_{\text{Total}}$) – the curvature of a surface calculated in two dimensions – and (2) tangential curvature ($C_{\text{Tan}}$) – the curvature calculated normal to the slope gradient . These two measures are employed to extract hilltop curvature and channel networks, respectively. However, these definitions vary between studies and software packages; see for a full review of the varying nomenclature and definitions of curvature measurements used in the literature.

Work by utilized 10 m resolution DEMs to evaluate the most accurate method for calculating curvature from digital topographic data. It was concluded that curvature could be most accurately calculated when a nine-term polynomial was fitted to the elevation surface, with the caveat that this will only be effective where the data quality is high enough. In cases in which the data are of lower accuracy, recommended using quadratics to fit the elevation data. This work was extended by to consider whether these patterns held for high-resolution topographic data, and it was found that fitting a six-term quadratic or nine-term polynomial yielded similar results. Therefore, chose to use the six-term quadratic to compute curvature. For this study we also chose to use the six-term quadratic in order to reduce computation time and, more importantly, to provide more robust curvature values as the data quality is degraded to resolutions below 10 m .

We calculate curvature using a circular window passed across the landscape, with a radius defined by identifying scaling breaks in the standard deviation and interquartile range of curvature calculated at increasing window sizes, consistent with the length scales of individual hillslopes . Consequently, curvature measurements on the hillslope scale can only be considered at data resolutions high enough to resolve individual hillslope features, considered here to be no more than 10 m, based on the window sizes identified for each landscape. A quadratic function of the form $$\mathit{\zeta }=a{x}^{\mathrm{2}}+b{y}^{\mathrm{2}}+cxy+dx+ey+f$$ is then fitted to the elevation values within the window by least squares regression , where $\mathit{\zeta }$ is the elevation, $x$ and $y$ are horizontal coordinates, and $a$ through $f$ are fitting coefficients. The fitted coefficients of this polynomial can be used to calculate differing types of curvature: $$C_{\text{Total}}=\mathrm{2}a+\mathrm{2}b$$ and $$C_{\text{Tan}}={\displaystyle \frac{\mathrm{2}a{e}^{\mathrm{2}}-\mathrm{2}cde+\mathrm{2}b{d}^{\mathrm{2}}}{(d^{\mathrm{2}}+e^{\mathrm{2}})\sqrt{(\mathrm{1}+d^{\mathrm{2}}+e^{\mathrm{2}})}}}.$$

From the measure of $C_{\text{Total}}$ for every cell in a DEM, we can also extract a subset of curvature values from the hilltops. The value of curvature at a hilltop ($C_{\text{HT}}$) can be readily evaluated if the positions of the hilltops are known. To extract hilltops we follow in defining a hilltop as the boundary between two drainage basins of the same stream order. These points in the landscape can be algorithmically extracted once a channel network is defined through the identification of points in the landscape where two channels of the same Strahler order meet and the identification of that point's upslope contributing area. Each of these areas defines a basin of a given order, and by repeating this process across the range of Strahler orders found in the landscape, a network of hilltops can be defined. This network is then used to sample the curvature values at these locations to provide the $C_{\text{HT}}$ values across the landscape. To ensure consistency between $C_{\text{HT}}$ measurements at changing grid resolutions, the same channel network, generated using the geometric method described in Sect.  from 1 m resolution data, is used as the basis of the hilltop extraction algorithm.

For our data on hilltop curvature, $C_{\text{HT}}$, hilltops with a gradient exceeding 0.4 are excluded as demonstrated that this gradient is the point at which $>\mathrm{\hspace{0.17em}15}\mathit{\%}$ of sediment transport is nonlinear. Under nonlinear sediment flux hilltop curvature scales nonlinearly with erosion rate and consequently cannot be used as a proxy for erosion rates. As hilltops have a convex form, their curvature should be negative, so as a final step any points identified as hilltops which have a positive curvature are excluded from further analysis.

Channel extraction

Extracting channel networks from digital topographic data remains a fundamental challenge for many areas of topographic analysis. Without the ability to discriminate between fluvial and hillslope domains, it is not possible extract many topographic metrics such as hillslope length , mean basin slope , or hilltop curvature , and the accuracy of each of these metrics will be influenced by the accuracy of the channel network extracted. At a more fundamental level, the ability to identify where channels initiate will facilitate better understanding of the processes acting at the transition between diffusive (hillslope) and advective (fluvial) sediment transport .

Many authors have made use of field-mapped channel heads both as a basis for geomorphic analysis and as a method for evaluating channel extraction methods . Prior to the availability of high-resolution topographic data, contributing area and slope-area scaling thresholds were commonly used to define the location of channel heads directly from DEMs . The influence of decreasing grid resolution on such channel extraction methods was evaluated by , who demonstrated a strong sensitivity in predicted channel head location to grid resolution, suggesting that coarser-resolution data may not be suitable for channel extraction through an area threshold. We apply the method described by to quantify the accuracy of an extracted channel network, detailed in Sect. .

Several methods have been proposed to identify channel heads from high-resolution topography. Typically these methods exploit the high-resolution nature of topographic data to resolve morphometric or process-based signatures of channel initiation or the transition between the hillslope and fluvial domain . Here we evaluate how two techniques – one geometric method built upon work by and and one process-based method, the DrEICH algorithm, developed by – are influenced by decreasing grid resolution.

The DrEICH method was selected for evaluation as the technique on which it is based has been shown to operate successfully in lower-resolution data . The DrEICH method makes use of $\mathit{\chi }$ analysis, performed by integrating drainage area along a river profile to facilitate comparisons between river profiles of differing drainage area, with fewer uncertainties than traditional slope-area analysis . When plotting the $\mathit{\chi }$ value against elevation for a river profile, river channels will plot as linear segments, whereas hillslopes will display nonlinear segments. The DrEICH algorithm identifies the transition between these linear and nonlinear segments as the best-fit location of the channel head.

The geometric method, used by , removes noise from the raw topographic data using a Wiener filter , as recommended by . This smoothed topography is then processed to identify channelized portions of the landscape using a tangential curvature threshold e.g.,, selected using the deviation of the probability density function of curvature from a normal distribution on a quantile–quantile plot (e.g., , and ). The identified areas of channelization are then combined into a contiguous channel network by employing a connected-components algorithm and thinned into a final channel network skeleton using the algorithm of .

Channels were extracted from the 5, 10, 20, and 30 m DEMs generated in Sect.  using both of the channel extraction methodologies. Parameters required in the operation of each algorithm were selected based on values used in previous studies , and these values can be found in Appendix .

Comparing channel networks

To assess the accuracy of the channel networks extracted using both methods, we employ two measures of quality described by . These measures operate on classifications of the predicted location of channel heads placing each channel head into one of three categories: true positives (TPs), false positives (FPs), and false negatives (FNs). A TP is where a predicted channel head from low-resolution data occupies the same spatial location as the channel head derived from 1 m resolution topography. An FP is where a predicted channel head is placed in a location where there is no channel head in the high-resolution data. An FN is when a channel head from high-resolution topography does not have a predicted channel head from low-resolution topography in the same spatial location.

We follow in employing a 30 m search radius around the 1 m derived channel heads and consider a low-resolution channel head falling within this radius to be spatially coincident. This has the effect of normalizing the size of each channel head point, to ensure that we can perform comparisons between predictions made at different spatial resolutions.

The reliability, $r$, of a channel extraction method is the ability of a method to not predict channel heads in areas where none are located and is calculated as $$r={\displaystyle \frac{\sum \text{TP}}{\sum \text{TP}+\text{FP}}},$$ where $\sum$ TP is the total number of true positives and $\sum$ FP is the total number of false positives. The sensitivity, $s$, of a channel extraction methodology is given by $$s={\displaystyle \frac{\sum \text{TP}}{\sum \text{TP}+\sum \text{FN}}},$$ where $\sum$ FN is the total number of false negatives. The sensitivity is the ability of a method to predict all of the channel heads expected. Using these two indexes it is possible to quantify the quality of channel heads predicted using low-resolution data, as well as understand why a particular method fails, by distinguishing between methods which fail due to either over- or underpredicting the number of channel heads in a landscape or by simply placing channel heads in the wrong spatial location.

Estimating the hillslope sediment transport coefficient from hilltop curvature

The sediment transport coefficient, $D$ [${\mathrm{L}}^{\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}$] (dimensions of mass [$\mathrm{M}$], length [$\mathrm{L}$], and time [$\mathrm{T}$] denoted in square brackets), of a landscape is a measure of its sediment transport efficiency and was demonstrated by to be controlled by the material properties of soil such as grain size, cohesion, and thickness. The value of $D$ within a landscape will exert a control on the morphology of hillslopes e.g.,. Diffusion-like hillslope evolution can be modeled using a 1-D mass conservation equation, assuming that the contribution to surface lowering from chemical processes is negligible when contrasted with the signal of surface lowering from physical processes (e.g., , and ): $${\mathit{\rho }}_{\text{s}}{\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}=-{\mathit{\rho }}_{\text{s}}{\displaystyle \frac{\partial q_{\text{s}}}{\partial x}}+{\mathit{\rho }}_{\text{r}}U,$$ where $\mathit{\zeta }$ $\left[\mathrm{L}\right]$ is the elevation of the land surface, ${\mathit{\rho }}_{\text{s}}$ and ${\mathit{\rho }}_{\text{r}}$ $\left[\mathrm{M}{\mathrm{L}}^{-\mathrm{3}}\right]$ are densities of dry soil and rock, respectively, and $U$ $\left[\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}\right]$ is the uplift rate. In steady-state landscapes, where $U=E$ and $\partial \mathit{\zeta }/\partial t=\mathrm{0}$, Eq. () simplifies to $${\displaystyle \frac{{\mathit{\rho }}_{\text{r}}}{{\mathit{\rho }}_{\text{s}}}}E={\displaystyle \frac{\partial q_{\text{s}}}{\partial x}},$$ with $E$ $\left[\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}\right]$ denoting the erosion rate. To solve this equation, a statement of the volumetric sediment flux per unit contour length, $q_{\text{s}}$ $\left[{\mathrm{L}}^{\mathrm{2}}{\mathrm{T}}^{-\mathrm{1}}\right]$, must be derived. A nonlinear relationship between sediment flux and topographic gradient has been proposed by a number of authors . Support for such models has been found from both tests of the resulting topographic predictions as well as through independent measurements of sediment flux across hillslopes .

The nonlinear model proposed by and is of the form $$q_{\text{s}}=\text{DS}{\left[\mathrm{1},-,{\left({\displaystyle \frac{\left|S\right|}{S_{\text{c}}}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}},$$ where $S_{\text{c}}$ is a critical gradient, and as the hillslope gradient approaches this threshold, $q_{\text{s}}$ asymptotes towards infinity.

At low hillslope gradients (e.g., on hilltops), the term within brackets in Eq. () approximates to unity . Equation () can therefore be substituted into Eq. () and can be solved for $D$ on low-gradient hilltops, assuming that an independent constraint on $E$ is available, $$D=-{\displaystyle \frac{E{\mathit{\rho }}_{\text{r}}}{C_{\text{HT}}{\mathit{\rho }}_{\text{s}}}}.$$

Hillslope length and relief

Hillslope length ($L_{\text{H}}$) is a crucial landscape parameter to constrain as it controls the rate of mass flux by overland flow within catchments , influences rates of soil erosion , and presents a first-order control on the maximum source area of landslides . Furthermore, it may be used to demonstrate nonlinearity in hillslope sediment flux .

Many studies have attempted to calculate hillslope length through the inversion of drainage density , analysis of plots of local slope against drainage area , direct measurements from topographic maps , and by measuring the length of overland flow from ridgeline to channel . demonstrated that the most geomorphologically suitable technique to use, particularly in the context of hillslope sediment transport, was that of measuring the length of overland flow. An additional measure which can be derived from this technique is the topographic relief, which is the difference in elevation between a hilltop and channel connected by a hillslope flow path. Topographic relief has been estimated in a number of ways and is frequently used in studies of tectonic geomorphology (e.g., , , , and ). Furthermore, topographic relief may be used to generate dimensionless erosion and relief plots , which can be used to identify landscape transience . Consequently, we intend to test the robustness of measuring hillslope length and relief as grid resolution decreases, with the aim of facilitating increased confidence in geomorphic analyses performed in locations where high-resolution topography is unavailable.

(a) Map showing the location of each of the study sites within the USA. (b–d) Shaded reliefs of representative sections of each study site, generated from 1 m resolution data. Tick spacing is in meters. All coordinates are in UTM. (b) Gabilan Mesa, California, UTM Zone 10${}^{\circ }$ N. (c) Santa Cruz Island, California, UTM Zone 11${}^{\circ }$ N. (d) Oregon Coast Range, Oregon, UTM Zone 10${}^{\circ }$ N.

[Figure omitted. See PDF]

Using the 20 topographic datasets generated in Sect.  for each of the three landscapes, hillslope length measurements were generated following the methods outlined in . We measured hillslope length on each dataset using two different channel networks. Firstly, channel heads were extracted from the highest-resolution data set, in each case 1 m, using the geometric method outlined in Sect. . These high-resolution channel heads were mapped onto the coarser-resolution topographic data, to ensure that changing channel extraction results will not have an influence on the measures of hillslope length. This allows improved isolation of the factors driving variations in hillslope length as grid resolution is decreased. Secondly, the analysis was performed using coarser-resolution channel networks extracted using the geometric method of channel extraction. We use the geometric method as opposed to the DrEICH method because, as we will show below, the geometric method is less sensitive to grid resolution. These two channel networks effectively provide upper and lower bounds for the accuracy of hillslope length and relief measurements.

Study sites

Three study sites from the United States have been selected for this study: Santa Cruz Island, California; Gabilan Mesa, California; and the Oregon Coast Range, Oregon. The first two sites have regularly spaced valleys at a range of length scales, particularly Gabilan Mesa, which has been the focus of previous work in this context . Santa Cruz Island, while less studied in the context of topographic analysis than Gabilan Mesa, has a wider range of hilltop curvatures (Fig. ). The Oregon Coast Range has been considered to be very regular, with uniform first-order drainage areas . However, more recent work has demonstrated the spatial variability of many topographic measurements in this landscape , and as such it provides a more challenging test case for our analyses. Furthermore, these sites were selected as they each have high-resolution lidar data covering a large spatial area and have been the subject of many previous studies , which should provide a good basis for the evaluation of the results of this study in a wider geomorphic context.

Gabilan Mesa

Gabilan Mesa, a section of the Central Coast Ranges in California, USA (Fig. b), is a highly regular landscape with very gentle transitions between hillslopes and channels, which correspond to topographic predictions of diffusion-like sediment transport . The area's semiarid climate supports a range of vegetation from oak savanna to chaparral shrubland . The nature of this lower-density vegetation allows a larger proportion of lidar pulses to reach the ground, requiring less processing and interpolation to generate a final bare-earth DEM for analysis .

A series of large, linear canyons running northeast to southwest are fed by parallel tributaries which flow perpendicular to the main trunk channel. These regularly spaced valleys present two distinct length scales in the landscape which have been observed both qualitatively and quantitatively through measurements of hillslope length distributions . Relationships between dimensionless erosion rate and relief, the uniformity of hilltop curvatures, and the regularity of valley spacing have all been used to assert that much of this landscape is in steady state , although localized observations of a relict plateau surface add complexity to this steady-state observation.

Maps showing the spatial variation in total curvature measurements as grid resolution is decreased for the same section of Santa Cruz Island as displayed in Fig. . All coordinates are in UTM Zone 11${}^{\circ }$ N. Panels (a)–(f) represent resolutions of 1, 2, 5, 10, 20, and 30 m. Tick spacing is in meters. The black boxes outline the same features as highlighted in Fig. , showing the reduction in the curvature signal with grid resolution for such a feature.

[Figure omitted. See PDF]

Santa Cruz Island

Santa Cruz Island (Fig. c), the largest of the eight California Channel Islands located to the west of California, USA, is divided by a large east–west trending valley, which follows the Santa Cruz fault . Parallel to this valley are two large ridges – one to the north and one to the south – which exhibit regularly spaced parallel channels draining north to south ; this regular pattern is particularly evident in the northwest section of the study area. The Santa Cruz Fault has been demonstrated to have left-lateral strike slip motion, which deflects channels away from the perpendicular to the main valley in the center of the island . Studies of marine terraces in the region suggest that the Channel Islands have been steadily uplifted through the late Quaternary .

The island has a Mediterranean climate similar to that of Gabilan Mesa , supporting extensive grassland with occasional patches of pine forest and chaparral vegetation . Human activities led to overgrazing across the island at the turn of the 19th century, causing a period of gullying and rapid erosion, particularly evident in the southwest of the island . The lidar data collected for this location have been extensively tested and ground truthed, ensuring that they are suitable for use in a geomorphic context and for performing topographic analysis at high spatial resolutions.

Plots of the distribution of $C_{\text{Total}}$ (a, c, e) and $C_{\text{Tan}}$ (b, d, f) measurements as resolution is decreased for each of the study landscapes. Whiskers are the 2nd and 98th percentiles; the box covers the 25th and 75th percentiles; the blue bar is the mean and the red bar is the median. The gray outline is the probability density function of each dataset.

[Figure omitted. See PDF]

Oregon Coast Range

The Oregon Coast Range in Oregon (Fig. d), USA, is a densely vegetated upland landscape, dominated by coniferous and hardwood forests , with a humid climate . Qualitative observations of the landscape suggest that the valleys are regularly spaced, with a particular uniformity found in the dimensions of first-order drainage basins . Such observations have been supported by measurements of hillslope length across the landscape . However, comparisons of the dimensionless relief and erosion rate performed by highlight the small-scale topographic variability inherent in this otherwise regular landscape. The Oregon Coast Range is considered to be in steady state due to the correlation between uplift rates from marine terrace data and erosion rates from cosmogenic radionuclides . The hillslopes are steeper and the ridgelines sharper than in Gabilan Mesa, consistent with observations of debris flows and shallow landsliding across the range , which have the potential to create a distinct topographic signature .

Results

Curvature

Figure  illustrates the variations in total curvature with grid resolution for a section of Santa Cruz Island. As the grid resolution is decreased, the range of $C_{\text{Total}}$ measurements are reduced, with much of the landscape becoming apparently planar. Within the black box, which covers the same spatial area as the boxes in Fig. , the impact of degrading resolution on small topographic features is observed, with the curvature signal of this first-order feature being lost as the grid resolution approaches 30 m.

Plots of the reduction in range between the 2nd and 98th percentiles (blue triangles) and the interquartile range (red circles) of $C_{\text{Total}}$ (a, c, e) and $C_{\text{Tan}}$ (b, d, f) measurements as resolution is decreased for each of the study landscapes.

[Figure omitted. See PDF]

Representative sections of each landscape's channel network displaying the extent of each network as grid resolution is decreased. Panels (a), (b), and (c) are generated using the DrEICH method of channel extraction. Panels (d), (e), and (f) are generated using the geometric method. All coordinates are in UTM. Tick spacing is in meters. The left column is from Santa Cruz Island, UTM Zone 11${}^{\circ }$ N, the central column is from Gabilan Mesa, UTM Zone 10${}^{\circ }$ N, and the right column is from the Oregon Coast Range, UTM Zone 10${}^{\circ }$ N.

[Figure omitted. See PDF]

Figure  displays the variations in the distribution of total and tangential curvature measurements with grid resolution for each of the study landscapes. Santa Cruz Island shows little variation in mean and median curvature with resolution, with the majority of the changes in each distribution with resolution occurring at the extremes of the curvature distribution for each dataset, as the representation of ridgelines and channel bottoms becomes increasingly diffuse. As resolution is decreased, the range between 2nd and 98th percentiles and the 1st and 3rd quartiles decreases, with a more rapid reduction in the more extreme values than in the quartiles (Fig. ). While this effect is most marked at the extremes, the distributions are condensed across all percentile intervals as grid resolution is increased beyond 3–4 m. This behavior is observed for both $C_{\text{Total}}$ and $C_{\text{Tan}}$ as grid resolution is decreased.

In the Oregon Coast Range for both measurements of curvature, there is little variation between the 1, 2, and 3 m datasets, with a broad range of measurements shown in the probability distributions. Beyond this point the mean and median do not significantly change, but as in Santa Cruz Island, the overall distribution of measurements compresses towards the average value for the landscape. The Gabilan Mesa data show similar trends to those of Santa Cruz Island but exhibit less variability at lower resolutions. The probability distributions of each measurement also exhibit less change with resolution than the other two datasets, indicating a reduced sensitivity to grid resolution at this location.

The variations in reliability (Eq. ) and sensitivity (Eq. ) of each channel network with decreasing grid resolution. Panels (a), (c), and (e) are generated using the geometric method of channel extraction. Panels (b), (d), and (f) are generated using the DrEICH method. The top row is from Gabilan Mesa, the middle row is from Santa Cruz Island, and the bottom row is from the Oregon Coast Range. The full results from this analysis can be found in Tables  and .

[Figure omitted. See PDF]

Channel networks

Figure  provides a qualitative overview of the changes of channel network extent with decreasing grid resolution for both methods, across the three test landscapes. In each case the general patterns are that as the grid resolution is decreased, the lowest-order channels are lost, as they exist on a spatial scale below that of the data resolution. In contrast, large parts of the predicted networks appear to occupy similar spatial locations in larger, higher-order channels where the topographic signal of a channel is more pronounced. The geometric method shows less reduction in drainage density than the DrEICH method, as data resolution is decreased.

Figure  provides a quantitative assessment of channel extraction quality by presenting the indexes of reliability and sensitivity for both the geometric channel extraction and extraction based on DrEICH, as the grid resolution is decreased. In Gabilan Mesa the channels extracted by the geometric method exhibit a high reliability which does not decrease considerably with decreasing grid resolution, suggesting that for each resolution step a large proportion of the predicted channel heads are spatially coincident with the channel heads generated from the 1 m data. The sensitivity values for this method and location are lower and decline more steadily with decreasing grid resolution, suggesting an increasing number of channel heads being missed by the algorithm as grid resolution is decreased. The DrEICH method does not perform as well in Gabilan Mesa, with lower index values for the 5 m data than the geometric method, and a rapid decline towards index values of 0, suggesting that the predicted channel heads bear little relation to the channel heads from the 1 m data.

In Santa Cruz Island the geometric method's reliability index is similar to Gabilan Mesa; however, the sensitivity index is not as high, which indicates that a large number of channel heads are being missed, but where a prediction is made, it is typically accurate. The DrEICH method exhibits a similarly large reliability initially but again shows more rapid degradation in the index value as grid resolution is decreased. The sensitivity values again decline more rapidly and reach a 0 value at 20 m grid resolution.

Changes in the estimated sediment transport coefficient, $D$, calculated using Eq. () and parameters in Table  for each of the three study landscapes, with decreasing data resolution. The error bars on each data point represent the uncertainties reported for each landscape's erosion rate data.

[Figure omitted. See PDF]

Published parameters used to calculate diffusivity.

Location	Soil density	Rock density	Erosion rate	Reference
	(kg m${}^{-\mathrm{3}}$)${}^{*}$	(kg m${}^{-\mathrm{3}}$)${}^{*}$	(mm yr${}^{-\mathrm{1}}$)
Santa Cruz Island	$1.4$	$2.4$	$0.069\pm 0.007$
Gabilan Mesa	$1.4$	$2.4$	$0.36\begin{array}{c}+0.38\\ -0.22\end{array}$
Oregon Coast Range	$1.4$	$2.4$	$0.1\pm 0.05$

${}^{*}$ Soil and rock densities are representative of typical measurements of the field sites and are taken from .

Reliability and sensitivity metrics for the DrEICH method of channel extraction.

Location	Resolution (m)	$\sum$ TP	$\sum$ FP	$\sum$ FN	$r$	$s$
Gabilan Mesa	$\mathrm{5}$	$555$	$982$	$1489$	$0.36$	$0.27$
	$10$	$210$	$879$	$1875$	$0.19$	$0.1$
	$20$	$42$	$734$	$2088$	$0.05$	$0.02$
	$30$	$13$	$609$	$2122$	$0.02$	$0.01$
Santa Cruz Island	$\mathrm{5}$	$3295$	$1971$	$4799$	$0.63$	$0.41$
	$10$	$2454$	$793$	$6865$	$0.76$	$0.26$
	$20$	$69$	$838$	$8235$	$0.08$	$0.01$
	$30$	$27$	$915$	$8284$	$0.03$	$0.0$
Oregon Coast Range	$\mathrm{5}$	$507$	$1718$	$1131$	$0.23$	$0.31$
	$10$	$144$	$445$	$1462$	$0.24$	$0.09$
	$20$	$16$	$105$	$1623$	$0.13$	$0.01$
	$30$	$\mathrm{2}$	$442$	$1639$	$0.0$	$0.0$

Reliability and sensitivity metrics for the geometric method of channel extraction.

Location	Resolution (m)	$\sum$ TP	$\sum$ FP	$\sum$ FN	$r$	$s$
Gabilan Mesa	$\mathrm{5}$	$1019$	$519$	$987$	$0.66$	$0.51$
	$10$	$712$	$380$	$1301$	$0.65$	$0.35$
	$20$	$448$	$332$	$1592$	$0.57$	$0.22$
	$30$	$292$	$333$	$1775$	$0.48$	$0.14$
Santa Cruz Island	$\mathrm{5}$	$4280$	$991$	$3109$	$0.81$	$0.57$
	$10$	$2473$	$777$	$4998$	$0.76$	$0.33$
	$20$	$334$	$505$	$7861$	$0.4$	$0.04$
	$30$	$475$	$470$	$7659$	$0.5$	$0.06$
Oregon Coast Range	$\mathrm{5}$	$792$	$1438$	$788$	$0.36$	$0.5$
	$10$	$562$	$602$	$938$	$0.48$	$0.37$
	$20$	$276$	$374$	$1275$	$0.42$	$0.18$
	$30$	$475$	$277$	$1418$	$0.38$	$0.11$

The data for the Oregon Coast Range show similar patterns for both methods, although the geometric method exhibits systematically larger index values. In each case the reliability increases slightly from 5 to 10 m resolution and then declines gradually towards 30 m resolution. The sensitivity indexes for both methods begin at a larger value than the reliability indexes and steadily decline towards 0. A sensitivity value exceeding the reliability value suggests that in this landscape there are fewer missed channel heads in the 5 m data but at the expense of too many predicted channel heads in locations where there are none predicted in the 1 m data.

Plots of the distribution of hillslope length (a, c) and relief (b, d) measurements as resolution is decreased for Santa Cruz Island. Whiskers are the 2nd and 98th percentiles; the box covers the 25th and 75th percentiles; the blue bar is the mean and the red bar is the median. The gray outline is the probability density function of each dataset. The top row presents the best-case scenario, where an independent constraint on the channel network is available for the lower-resolution data, and the bottom row uses the channel networks extracted using the geometric method outlined in Sect.  for each resolution step.

[Figure omitted. See PDF]

Sediment transport coefficient

Using the values for hilltop curvature generated in Sect. , published parameters for erosion rate and material properties outlined in Table  and Eq. (), the average sediment transport coefficient ($D$) of each landscape can be calculated as a function of grid resolution. Figure  displays the relationship between diffusivity and grid resolution for each of the three study sites. The data for Santa Cruz Island and Oregon Coast Range both show a gradual increase in diffusivity with decreasing grid resolution, the rate of which reduces with decreasing grid resolution. The Gabilan Mesa data do not exhibit the same trend, with little variability in calculated $D$ values as resolution is decreased. Although the Oregon Coast Range and Santa Cruz Island datasets exhibit an increase in estimated $D$, all of the values for each location fall within the range of values for $D$ compiled by .

Hillslope length and relief

The hillslope length measurements for Santa Cruz Island calculated using 1 m channel heads (Fig. a) show little variation in the distribution of the data up to 10 m resolution, with the main difference being the decrease with grid resolution in the 2nd percentile measurements, which is a trend observed within each of the datasets. The mean and median values also gradually decrease towards the 10 m resolution dataset, before gradually increasing towards the 30 m resolution step. However, these variations are very small, with the overall distributions of hillslope length and relief not varying considerably between resolution steps. When the same hillslope length algorithm is applied using channel networks extracted using the geometric method for each resolution step (Fig. c), there is little change in the distribution or average values of $L_{\text{H}}$ until beyond the 10 m resolution step. Beyond this point the measurements of hillslope length are clearly affected by the reduction in accuracy of the channel network. The relief measurements for both channel head methods (Fig. b, d) in Santa Cruz Island exhibit little resolution dependence up to 10 m grid resolution, beyond which point the values increase steadily. In the case of the 1 m channel heads, the distribution becomes compressed around the average values at lower resolutions, whereas with the variable channel head dataset the distribution of values increases with decreasing resolution.

In Gabilan Mesa the hillslope length measurements calculated using 1 m channel heads (Fig. a) show a gradual reduction in mean and median values between the highest-resolution data and the 8 m resolution data before a small plateau and then a small increase until the 30 m dataset. The average relief values calculated for the same dataset increase steadily by approximately 20 m between the highest- and lowest-resolution datasets (Fig. b). The distribution of relief measurements are broadly consistent between 1 and 5 m resolutions before reducing about the median as grid resolution is decreased. The same trends are apparent in the hillslope length and relief data calculated using the variable channel heads (Fig. c, d) with little change between the two pairs of datasets.

Plots of the distribution of hillslope length (a, c) and relief (b, d) measurements as resolution is decreased for Gabilan Mesa. Whiskers are the 2nd and 98th percentiles; the box covers the 25th and 75th percentiles; the blue bar is the mean and the red bar is the median. The gray outline is the probability density function of each dataset. The top row presents the best-case scenario, where an independent constraint on the channel network is available for the lower-resolution data, and the bottom row uses the channel networks extracted using the geometric method outlined in Sect.  for each resolution step.

[Figure omitted. See PDF]

Plots of the distribution of hillslope length (a, c) and relief (b, d) measurements as resolution is decreased for the Oregon Coast Range. Whiskers are the 2nd and 98th percentiles; the box covers the 25th and 75th percentiles; the blue bar is the mean and the red bar is the median. The gray outline is the probability density function of each dataset. The top row presents the best case scenario, where an independent constraint on the channel network is available for the lower-resolution data, and the bottom row uses the channel networks extracted using the geometric method outlined in Sect.  for each resolution step. At higher-resolution steps the 98th percentile data is not shown in the plot, to better highlight the distribution of measurements between the 25th and 75th percentiles, which make up the majority of the data points.

[Figure omitted. See PDF]

The hillslope length measurements for the Oregon Coast Range with channel heads from the 1 m data (Fig. a) again show a gradual reduction in the median values with a gradual increase in the mean values until 20 m grid resolution. Beyond this point the data become considerably more variable, with a large increase in both the mean and median results. The relief data shown in Fig. b are the most consistent of the three landscapes, with very little variation in the values until they begin increasing with grid resolution at approximately 20 m resolution. The data presented in Fig. c and d show the most sensitivity to grid resolution of the three landscapes. Average hillslope length values reduce towards 10 m before stabilizing and then rapidly increasing in the same manner as the fixed channel head data. The relief measurements show a gradual decline in mean relief across the range of resolutions from 1 to 10 m, where the fixed data show much less variation.

Discussion

Curvature and the problem of resolution-dependent filtering

Across the three landscapes the variance of the distributions of both total and tangential curvature values are systematically reduced as resolution is decreased, an effect that is particularly notable after the grid resolution exceeds 3–4 m (Fig. ). In each of the three datasets, the interquartile ranges remain relatively constant, whereas beyond 4 m resolution in each case the range between the 2nd and 98th percentiles reduces rapidly (Fig. ), demonstrating that the majority of the loss of curvature information occurs at the extremes of the distribution.

In producing a DEM, we are sampling a complex two-dimensional elevation signal, in which spatial variations in geomorphic processes drive variations in topographic amplitude at different wavelengths . Decreasing the grid resolution of DEMs acts as a low-pass filter on this topographic signal, which preferentially degrades features in the topography that have significant amplitude at small wavelengths, such as sharp ridgelines, narrow valley bottoms, and local topographic roughness generated by, for example, landslides, tree throw, and rock exposure (Figs.  and ). While the position of ridges and valleys is preserved in coarser-resolution data, the magnitude of their associated curvature values is reduced as resolution decreases; this effect is particularly marked for hillslopes in which curvature is focused at the ridge crest and valley bottoms, a common characteristic of more rapidly eroding landscapes . For first-order landscape features, such as gullies, landslide scars, and first-order channels, decreasing grid resolution eventually results in the complete loss of topographic information, as highlighted in Figs.  and .

Topographic filtering and its implications for curvature and slope measurements

We can explain some of the observed behavior in Figs.  and  through spectral analysis. Spectral analysis assumes that data can be approximated as the sum of sine waves of varying frequency. One can apply a spectral filter to any dataset: this simply means that one transforms input data into output data using linear functions (that is, we can multiply the input data by a series of weights). Any filter will have a gain, which is the ratio between the filtered amplitude and the original amplitude. A filter will also have a fidelity, which is the ratio between the continuous gain and the discrete gain. We are using discrete data, so the fidelity measures how well our discrete filter is able to reproduce a theoretical signal that is continuous. We can never have continuous data since lidar is not continuous: our filters will always represent an imperfect version of nature and fidelity quantifies just how imperfect it is. Hopefully our readers will not be put off by this foray into jargon, and we can move on to practical application of spectral filters for use in topographic applications.

We will examine the spectral behavior of a simplified one-dimensional system. We acknowledge that a 1-D approach cannot fully describe complex two-dimensional topography of real landscapes, but a one-dimensional system is amenable to mathematical treatment that can at least give us qualitative insight into trends observed in our data. In addition, some of the features of interest, for example ridgelines and channels, can be roughly approximated as one-dimensional structures within a two-dimensional landscape.

Curvature in one dimension, $C_x$ [${\mathrm{L}}^{-\mathrm{1}}$], is often approximated with the differencing equation: $$C_x={\displaystyle \frac{{\mathit{\zeta }}_{(x-\mathrm{\Delta }x)}-\mathrm{2}{\mathit{\zeta }}_x+{\mathit{\zeta }}_{(x+\mathrm{\Delta }x)}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}}},$$ where $\mathit{\zeta }$ [$\mathrm{L}$] is the elevation of the land surface, $x$ [$\mathrm{L}$] is a location in space, $C_x$ is the curvature at location $x$, and $\mathrm{\Delta }x$ [$\mathrm{L}$] is the grid interval. The subscripts denote the discrete locations where elevation is evaluated. Equation () is in fact a spectral filter. The original data is $\mathit{\zeta }$, which is distributed in space, and the weights in the filter are $(\mathrm{\Delta }x{)}^{-\mathrm{2}}$, $-\mathrm{2}(\mathrm{\Delta }x{)}^{-\mathrm{2}}$, and $(\mathrm{\Delta }x{)}^{-\mathrm{2}}$ for data points at $(x-\mathrm{\Delta }x)$, $x$, and $(x+\mathrm{\Delta }x)$, respectively. From this filter, we can calculate the wave number response function. A full description of the theory and significance of a wave number response function can be found in . For our purposes, it is sufficient to know that this function must be calculated if we are to calculate the gain and fidelity of the filter (which here is a measure of curvature of our elevation data). The wave number response function ($H(\mathit{\omega };\mathrm{\Delta }x)$) from this filter, given by in their Eq. (7.3.7), is $$H(\mathit{\omega };\mathrm{\Delta }x)={\displaystyle \frac{\mathrm{2}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}}}\left[\cos \right(\mathit{\omega }\mathrm{\Delta }x)-\mathrm{1}],$$ where $\mathit{\omega }=\mathrm{2}\mathit{\pi }/L[{\mathrm{L}}^{-\mathrm{1}}$] is the wave number with wavelength $L$ [$\mathrm{L}$]. Higher wave numbers correspond to shorter wavelengths. Using this function, we can calculate the gain, $G(\mathit{\omega };\mathrm{\Delta }x)$. Again, the gain measures the ratio of the amplitude of the filtered signal (in this case curvature) to the amplitude of the original signal (in this case elevation) at the wave number $\mathit{\omega }$. The theoretical gain for continuous waveforms of curvature (i.e., not discrete filters like Eq. ) is ${\mathit{\omega }}^{\mathrm{2}}$. The gain of a discrete filter is the modulus of the wave number response function see p. 296 in, so in the case of Eq. () the resultant gain, $G(\mathit{\omega };\mathrm{\Delta }x)$ is $$G(\mathit{\omega };\mathrm{\Delta }x)={\displaystyle \frac{\mathrm{2}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}}}[\mathrm{1}-\mathrm{2}\cos (\mathit{\omega }\mathrm{\Delta }x)+{\cos }^{\mathrm{2}}(\mathit{\omega }\mathrm{\Delta }x){]}^{\mathrm{1}/\mathrm{2}}.$$

In the case of our curvature filter (Eq. ), the gain function reveals how high-frequency waveforms (e.g., ridge crests, tree throw mounds, local roughness) in the elevation data involve relatively large values of curvature, whereas low-frequency elevation waveforms (e.g., ridge–valley features or geologic folds) with the same amplitude involve relatively small curvatures. Crucially, however, the discrete filter does not retain all of the high-frequency information. Some of this information is lost in the discretization process (i.e., it is lost because we are sampling the data at fixed intervals rather than having continuous information about the surface). We can calculate what information is lost by calculating the fidelity, which is the ratio between discrete gain (Eq. ) and the theoretical gain (${\mathit{\omega }}^{\mathrm{2}}$): $$\begin{array}{cc}& & F(\mathit{\omega };\mathrm{\Delta }x)=& {\displaystyle \frac{\mathrm{2}}{(\mathrm{\Delta }x{)}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}}[\mathrm{1}-\mathrm{2}\cos (\mathit{\omega }\mathrm{\Delta }x)+{\cos }^{\mathrm{2}}(\mathit{\omega }\mathrm{\Delta }x){]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Plot of fidelity ($F$) of two one-dimensional differencing operations: curvature (Eq. ) and topographic gradient (Eq. ) as a function dimensionless wave number $\mathrm{\Delta }x/L$ to the Nyquist wave number, $\mathrm{\Delta }x/L=0.5$.

[Figure omitted. See PDF]

Again, fidelity is a measure of how closely our discrete filter (here curvature measured at discrete points in the landscape) reflects the true curvature (that is, the curvature measured if we had a perfectly continuous dataset). Fidelity is a function of the ratio between the grid interval and the wavelength (Fig. ). When the fidelity is unity, the discrete filter exactly reproduces the underlying continuous function. Again, the landscape (and its derivative metrics like curvature and gradient) has features at different wavelengths, such as long-wavelength ridges and valleys and short-wavelength tree throw mounds.

As the frequency approaches the Nyquist wave number, defined as $\mathrm{\Delta }x/L=\mathrm{1}/\mathrm{2}$, fidelity decreases (Fig. ); a fidelity of only approximately 0.4 is achieved at the Nyquist wave number itself. To achieve a fidelity, $F$, of 0.9 requires that $L/\mathrm{\Delta }x$ is equal to approximately six grid points per wavelength. A fidelity $F=0.95$ requires eight points per wavelength, and $F=0.99$ requires 18. Therefore, while the grid resolution imposes a minimum wavelength that can be resolved (defined by the Nyquist wave number), the behavior of the fidelity function (Fig. ), clearly illustrates that curvature information will be lost when calculated for features with wavelengths greater than but still close to the minimum resolvable at the Nyquist wave number.

What does this mean in practical terms? In our simple, one-dimensional example, if we use 1 m resolution data we can only capture the curvature of a one-dimensional ridgeline that had a wavelength of 3–4 m (one does not need the entire wave to capture the peak of the waveform) but with a loss of fidelity on the magnitude of the curvature. Or, in other words, we would underestimate the magnitude of the curvature.

Another landscape metric that is widely measured is topographic gradient. In our study we have not computed how topographic gradient varies as a function of grid resolution because this has been examined by many previous authors (e.g., , , and ). However, our treatment of the properties of a one-dimensional filter can give some insight into previous results. Consider a simple central-difference approximation of the topographic gradient ($S_x$, dimensionless): $$S_x={\displaystyle \frac{{\mathit{\zeta }}_{(x+\mathrm{\Delta }x)}-{\mathit{\zeta }}_{(x-\mathrm{\Delta }x)}}{\mathrm{2}\mathrm{\Delta }x}}.$$

Equation () is yet another spectral filter, with weights of $\mathrm{2}(\mathrm{\Delta }x{)}^{-\mathrm{1}}$ at $x+\mathrm{\Delta }x$ and $-\mathrm{2}(\mathrm{\Delta }x{)}^{-\mathrm{1}}$ at $x-\mathrm{\Delta }x$. We can follow the same series of operations that we performed on Eq. () to arrive at the fidelity of Eq. (), denoted as $F_{\text{S}}$, taking into account that the theoretical gain is $\mathit{\omega }$ see Eq. 7.3.8 in: $$F_{\text{S}}(\mathit{\omega };\mathrm{\Delta }x)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x\mathit{\omega }}}\left[\sin \right(\mathit{\omega }\mathrm{\Delta }x\left)\right].$$

Equation () formally illustrates why estimates of slope tend to systematically decrease with increasing grid interval $\mathrm{\Delta }x$ (Fig. ). Namely, an increasing $\mathrm{\Delta }x$ is able to resolve less local (high wave number) elevation structure while picking out the slope of more regional structure. The fidelity increases as the ratio of the grid interval to the wavelength, $\mathrm{\Delta }x/L$, decreases (Fig. ). To achieve a fidelity $F_{\text{S}}=0.9$, for example, requires $L/\mathrm{\Delta }x$ or approximately eight grid points per wavelength. A fidelity $F_{\text{S}}=0.95$ requires 11 points per wavelength, and $F_{\text{S}}=0.99$ requires 18. The fidelity of the one-dimensional gradient operator goes to 0 when approaching the Nyquist wave number ($\mathrm{\Delta }x/L=\mathrm{1}/\mathrm{2}$). These results explain the pronounced loss of gradient information in coarse-resolution data observed by many authors (e.g., , , and ).

Total and tangential curvature

Having explored simplified one-dimensional filters, we now return to our two-dimensional results. Although real landscapes are two-dimensional and we use polynomial fitting rather than simple differencing as in Eq. (), we can still use Eq. () as a qualitative indicator of the grid resolution required for appropriate curvature estimates. In the Gabilan Mesa, where ridgelines are broad, lower-resolution data can still capture the curvature with relatively high fidelity. However, in locations with sharper ridgelines, such as Santa Cruz Island, the narrowest ridgelines are no longer adequately resolved as the grid resolution is decreased, as can be seen in Fig. .

The loss of fidelity predicted by the simple one-dimensional system (Eq. ) qualitatively predicts the pattern observed in Figs.  and , namely that the curvature values are smeared over a greater length scale leading to apparently broader ridges with resolution and a systematic underestimation of their peak elevations. This highlights that in conjunction with data quality, landscape morphology also exerts a control on the optimal resolution to use for a given study, where landscapes with more gradual hillslope to valley transition morphologies can be analyzed using coarser-resolution topographic data with more confidence. Although the identification of landscape morphology is often achieved through observations of high-resolution topography, it can be achieved through field observations and the use of ancillary datasets, which allow the qualitative checking of results obtained from a low-resolution dataset.

Santa Cruz Island and the Oregon Coast Range have the highest tangential curvature at 1 m resolution. High tangential curvature at Santa Cruz Island corresponds to observations of extensive gullying and hillslope erosion . In the Oregon Coast Range, features such as pit and mound topography produced by tree throw and other biotic activity are resolved in the lidar dataset , which manifests itself as an increase in values of curvature. However, this could also be indicative of non-topographic noise in the DEM surface produced during the processing of the point clouds, which is particularly required in heavily forested locations such as the Oregon Coast Range. This suggests an unfortunate collinearity between the two causes of small-wavelength topographic noise and warrants further testing in future to disentangle synthetic and natural noise from high-resolution topographic measurements. However, high curvature is not solely a manifestation of stochastic disturbance in local topographic roughness but is also generated at narrow valley bottoms and at ridgelines where erosion rates are rapid relative to the hillslope sediment transport coefficient . Gabilan Mesa exhibits much lower curvature values than the other two locations, which is a consequence of high landscape diffusivity, indicating that sediment transport at Gabilan Mesa is dominated by diffusion-like processes , smoothing the landscape and reducing the tangential curvature of the hillslope surface.

Channel extraction

It is intuitive to consider that when extracting channel networks at any data resolution, regardless of method, the higher-order, larger channels will be more accurately constrained than lower-order channels. This pattern is observed in each of the study landscapes, with the majority of the variations in channel locations occurring in first- and second-order channels. Such loss of low-order channels from datasets has implications for studies focusing on upland areas, in particular where detailed measurements which depend on channel network position are performed.

The contrast between the extent of channel networks and their indexes of quality for the two methods outline that a geometric method of channel extraction outperforms the process-based DrEICH algorithm. Due to the relative simplicity of the geometric method of channel extraction, errors inherent in the DEM are not compounded on the same scale as the DrEICH algorithm, which performs more operations on topographic data. As the geometric method identifies channels based on their tangential curvature, although channel head features may be smoothed out of the DEM as resolution is decreased, the channel will still express some positive curvature in lower-resolution data. The initiation point may be located downslope of the true channel head but even in this worst case most of the channel network will be extracted correctly. This is observed in Fig.  which shows a gradual reduction in drainage density as the grid resolution is decreased.

The indexes of quality defined by provide a clear framework to understand the quality of channel head predictions using these two methods as data resolution is decreased. In each case, the geometric method outperforms the DrEICH method, both in the accuracy of the channel heads which are predicted, and in the ability of the method to not predict channel heads in locations where no channel exists. These indexes are influenced by the size of the search radius around each channel head, and reducing this radius would decrease the index values. However, the use of a 30 m search radius allows comparisons to be drawn between predictions made at different data resolutions, and also between this study and that of .

This assessment of high-resolution methods with degraded-quality data demonstrates the ongoing challenges that channel extraction poses to the geomorphology community. performed extensive testing on channel extraction using threshold channel extraction methods and demonstrated similar limitations when channels were extracted using lower-resolution data. Our results suggest that a geometric method of channel extraction will provide an optimal channel network as data quality is reduced, particularly in uniform landscapes such as Gabilan Mesa. However, the only way to ensure the highest-quality results is to employ high-resolution data in conjunction with field mapping of channel network extents.

Sediment transport coefficient

The predicted values of the sediment transport coefficient ($D$) for the 1 m data fall within the range of values compiled by and estimated for the Oregon Coast Range and Gabilan Mesa by and . This suggests that this method can produce useful estimates of $D$ when employing high-resolution topography.

The sediment transport coefficients calculated at the Oregon Coast Range and Santa Cruz Island locations both increase with grid resolution, reflecting the sensitivity of $C_{\text{HT}}$ to grid resolution in each of these locations. Despite the Oregon Coast Range eroding 45 % more rapidly (Table ) than Santa Cruz Island, the rate of increase in $D$ measurements remains similar between the two landscapes. Gabilan Mesa data are generally insensitive to a decrease in grid resolution, as the scale of hilltop widths measured in Gabilan Mesa is on the order of tens of meters. This allows datasets with grid resolutions approaching half the width of a hilltop to provide an accurate estimate of hilltop curvature and, thus, the sediment transport coefficient.

These data suggest that estimating $D$ from low-resolution topographic data is possible in many landscapes, particularly those which have average ridgelines broader than the grid resolution of the topographic data. In the case of landscapes with sharper ridgelines such as Santa Cruz Island and the Oregon Coast Range, it is more challenging to constrain $D$ effectively as the grid resolution is decreased. The magnitude of the overestimation of $D$ between the highest- and lowest-resolution diffusivity estimates, 0.0023 m${}^{\mathrm{2}}$ a${}^{-\mathrm{1}}$ in the case of the Oregon Coast Range, will be a product of the uncertainty within the calculation of the erosion rate and material densities in addition to the local variations of $D$ within each landscape.

Hillslope length and relief

Measurements of hillslope length and relief have been used to test sediment flux laws and to identify landscape transience . Such analyses have previously been restricted to high-resolution topographic data. When considering hillslope length, we must select a grid resolution that is at least half the median hillslope length in order to resolve any useful information. However, in reality more than two pixels are required if any meaningful information is to be extracted from topographic data. As the median hillslope length for many landscapes has been shown to be in excess of 100 m , this requirement for several pixels per hillslope falls well within the range of many lower-resolution data products. Therefore, our results show that meaningful hillslope length measurements can be made from lower-resolution topographic data, with data products approaching 30 m resolution proving suitable in some cases.

The relief measurements for each landscape, however, show more sensitivity to grid resolution, with a systematic increase in the median values in each location beyond 10 m grid resolution. As decreasing grid resolution acts as a low-pass filter on the landscape, the elevation of ridges are expected to be reduced, whilst the elevation of channel beds are raised, producing a net reduction in topographic relief. However, the increased relief observed with decreasing grid resolution is produced by the decrease in drainage density with decreasing resolution observed in Fig. ; this produces fewer channels reaching up towards ridgelines and leading to hillslope flow paths traveling further downslope before reaching a channel.

By contrasting the $L_{\text{H}}$ and $R$ results computed using fixed and variable channel heads, it is clear that the optimal method for measuring hillslope length and relief is to employ as accurate a channel network as possible. However, the variable channel head data show that the signal of average hillslope length and relief is broadly insensitive to data resolution up to grid resolutions of at least 10 m. This would facilitate the analysis of landscape transience using these measurements on a global scale, using high-resolution satellite-derived DEMs, such as TanDEM-X . This relationship is again strongest in Gabilan Mesa, the landscape with the least topographic complexity which demonstrates the least sensitivity to curvature measurements and the estimation of diffusivity. However, even in the more noisy landscape of the Oregon Coast Range, meaningful hillslope length and relief measurements can still be made through the use of a geometric channel extraction algorithm and lower-resolution topographic data.

Conclusions

Through the generation of topographic data spanning the range of grid resolutions currently used in much of geomorphic research, a number of key metrics have been evaluated for their sensitivity to grid resolution. We have demonstrated the reduction in the range of total and tangential curvature values as grid resolution is decreased, across three test landscapes. These curvature measurements are important in the estimation of the hillslope sediment transport coefficient ($D$), in their use as a proxy for erosion rate, and in the extraction of channel networks from topographic data. We demonstrate that the estimation of $D$ from low-resolution topographic data is possible, particularly in landscapes such as Gabilan Mesa where hilltops are broad. Higher resolutions are required to extract meaningful curvature information in steep landscapes with sharp ridges and narrow gullies.

The extraction of channel networks from digital topographic data is a significant challenge on all spatial scales, as the definition of a channel network is integral in the execution of many analyses (e.g., , , and ). We demonstrate that the use of a geometric channel extraction algorithm produces channel networks for all three of our landscapes which correspond well to networks extracted from high-resolution topography. This correspondence is tested through the computation of quality indexes for each predicted network, which outline the suitability of this algorithm over a process-based method at coarse DEM resolutions.

Average values of hillslope length and relief for each landscape are shown to be broadly insensitive to grid resolution up to grid resolutions which correspond to the highest-resolution topographic data globally available. This indicates that these measurements can be used to identify landscape transience in locations where lidar data are unavailable. The accuracy of these measurements is dependent on the accuracy of the channel network used, however, as using a geometric method of channel extraction from the 1 m DEM still provides robust measurements of hillslope length and relief.

The relationships between decreasing grid resolution and the geomorphic parameters explored here demonstrate the influence of the spatial scale of the topographic expression of process on the quality of results which can be extracted from lower-resolution topography. From these analyses it is challenging to identify a clear threshold below which data become unsuitable for use in geomorphic analysis. Rather, it is important to highlight the influence of landscape morphology and the dominant processes acting upon it in the selection of an appropriate data resolution for a study. Using this work as a framework, it is now possible to place constraints on the accuracy of results derived from coarse-resolution topographic data, particularly where non-topographic or field data can be used to provide insight into general landscape morphology.

Code availability

All of the code used in this analysis is open source and the topographic analysis routines are available at http://github.com/LSDtopotools/LSD_Resolution; the code to generate the figures in this paper, alongside the raw plot data, can be downloaded from http://github.com/sgrieve/Resolution_Paper_Figs.

Data availability

The topographic data used in this study are freely available from http://www.OpenTopography.org, and the specific point clouds used can be downloaded from http://hdl.handle.net/10283/2071.

Channel extraction parameters

This table provides the parameters used to generate channel networks both using the geometric method and the DrEICH method. The drainage area value is used to thin the initial extracted network by removing channels which have a drainage area below the threshold value. The connected-components value defines the point at which a group of contiguous channel pixels are considered to be connected. The $\frac{m}{n}$ ratio is determined using software provided by , and its use within this context is discussed in detail in .

Parameters used by the geometric and process-based techniques in the extraction of channel networks.

Location	Window	Drainage	Connected	$\frac{m}{n}$	Reference
	radius ($m$)	area ($m^{\mathrm{2}}$)	components (Pixels)	ratio
Santa Cruz Island	4	4	5	0.50	This study
Gabilan Mesa	5	4	5	0.45
Oregon Coast Range	4	4	5	0.45

Stuart W. D. Grieve, Simon M. Mudd, David T. Milodowski, and Fiona J. Clubb wrote the software. Stuart W. D. Grieve performed the analysis. David J. Furbish and Simon M. Mudd resurrected the spectral filtering analysis from an unpublished 2002 manuscript because they are lovers of the long game. Stuart W. D. Grieve wrote the paper with contributions from the other authors.

Acknowledgements

Simon M. Mudd and Stuart W. D. Grieve are funded by NERC grant NE/J009970/1 and Simon M. Mudd is funded by US Army Research Office contract number W911NF-13-1-0478. Fiona J. Clubb is funded by the Carnegie Foundation for the Universities of Scotland. David T. Milodowski was funded by a NERC Doctoral Training Grant NE/152830X/1 and NE/J500021/1. David J. Furbish was funded by US National Science Foundation grant EAR-1420831. We thank Marie-Alice Harel, Kristin Sweeney, Wolfgang Schwanghart, and two anonymous reviewers for comments on earlier versions of this manuscript. Edited by: R. Gloaguen Reviewed by: K. Sweeney, W. Schwanghart, and two anonymous referees