Background

This research builds on numerous previous studies that investigated the accuracy of interpolation techniques used to develop DEMs (e.g., Aguilar et al., 2005; Akkala, Devabhaktuni, and Kumar, 2010; Ali, 2004; Anderson, Thompson, and Austin, 2005; Burrough and McDonnell, 1998; Carlisle, 2005; Chaplot et al., 2006; Declercq, 1996; Desmet, 1997; Erdogan, 2009, 2010; Fisher and Tate, 2006; Guo et al., 2010; Li and Heap, 2008; MacEachren and Davidson, 1987; Smith, Holland, and Longley, 2005). Studies have found that the accuracy of all interpolation techniques is related to the sampling density and distribution of measurements (Aguilar et al., 2005; Anderson, Thompson, and Austin, 2005; Chaplot et al., 2006; Erdogan, 2009, 2010; Guo et al., 2010; MacEachren and Davidson, 1987). In these previous studies, sampling density referred to either a percentage of original measurements (Aguilar et al., 2005; Anderson, Thompson, and Austin, 2005; Guo et al., 2010; MacEachren and Davidson, 1987) or a count of measurements per area (Chaplot et al., 2006; Erdogan, 2009, 2010). In this study, cell sampling density is defined as the percentage of the total number of cells in a DEM constrained by depth measurements at a given cell size. Previous studies have also found that the accuracy of all interpolation techniques is related to terrain characteristics (Aguilar et al., 2005; Carlisle, 2005; Erdogan, 2009, 2010; Guo et al., 2010).

While almost all previous studies focus on the accuracy of interpolation techniques in developing DEMs from dense topographic data, this research is unique in that it addresses the accuracy of interpolation techniques used to develop bathymetric DEMs from sparse depth measurements. The morphology of the seafloor differs from the morphology of land above sea level because of differences in natural processes including sedimentation and erosion. The differences between topographic and bathymetric DEMs regarding the cell sampling density, subsequent interpolation distances, and in terrain characteristics, make it worthwhile to specifically assess the accuracy of interpolation techniques in developing bathymetric DEMs. Furthermore, this research addresses the current lack of literature on predictive equations of cell-level interpolation uncertainty.

Interpolation

Interpolation is a mathematical process of predicting unknown values, such as ocean depths, using surrounding measured values (Burrough and McDonnell, 1998). Several different interpolation techniques exist for creating DEMs, and all techniques are based on the same assumption: that bathymetry is spatially autocorrelated. The notion of spatial autocorrelation is largely attributed to Tobler's 1^st Law of Geography, “Everything is related to everything else, but near things are more related than distant things” (Tobler, 1970). That is, the depth at one location is more similar to depths nearby than the depths far away. However, the different mathematical algorithms used in each interpolation technique produce divergent DEMs, even when developed from the same source data (Aguilar et al., 2005; Ali, 2004; Declercq, 1996; Erdogan, 2009).

Interpolation techniques can be classified into general groups based on the assumptions and features used to estimate the depths of unknown areas using surrounding known measurements. Interpolation techniques can be classified as any combination of geostatistical or deterministic, local or global, and exact or inexact (Li and Heap, 2008). Geostatistical techniques, such as kriging, use both mathematical and statistical functions to estimate depths, while deterministic techniques, such as spline, IDW, and triangulation, use the measurements directly and only mathematical functions to predict unknown values (Childs, 2004). Local techniques use a subset of measurements surrounding the location to be predicted and global techniques use all available measurements (Burrough and McDonnell, 1998). Last, exact interpolators always honor the measurements by creating a surface that has the same values at the locations of measurements, in contrast to inexact interpolators, which are not constrained by the measurement values at those locations (Burrough and McDonnell, 1998). The interpolation technique should be chosen based on the measurement uncertainty, cell sampling density, cell sampling distribution, terrain characteristics, and computational resources. More specifically, each interpolation technique has particular mathematical constraints for predicting unknown values.

Study Area

The bathymetry of Kachemak Bay, Alaska (Figure 1) was used to quantify the accuracy of the interpolation techniques. Kachemak Bay was chosen because of its varying terrain that ranges from areas of constant depths to areas of rapidly changing depths, and because of the availability of high-resolution bathymetric data. An area of the bay with heterogeneous terrain was chosen in order to evaluate the effect of terrain on interpolation accuracy, with the initial assumption that interpolation accuracy would decrease in areas of high terrain slope and curvature. Kachemak Bay is fed by glacial runoff, and the area surveyed has a variety of glacial sediments covered by a layer of sand that produces several types of bedforms and changes in depths (Bouma et al., 1980).

Figure 1.

Kachemak Bay, AK (derived from DEM of Friday et al., 2010). The red box in the bay outlines the study area (see Figure 2).

Interpolation Techniques

IDW predicts unknown values using a linearly weighted combination of sample points (Liu et al., 2007). The weight is a function of inverse distance and is raised to a user-defined mathematical power that controls the significance of known measurements based on their distance from the location to be predicted. A larger power results in distant measurements having less influence on the interpolated values, which produces a less smoothed surface. The number of measurements used for prediction can also be defined, with more depths resulting in smoother surfaces and more similar to global interpolation.

There are numerous types of spline interpolation, including thin plate spline, regularized spline, and spline with tension. Regularized spline was evaluated in this study. The regularized technique creates a smooth, gradually changing surface and allows for interpolated values outside the range of the measurement values (ESRI, 2010b). A user can further fine-tune this spline interpolation technique by the number of points and a weight parameter. The greater the number of points specified, the greater the influence of distant points, and the smoother the interpolated surface. For regularized spline, a greater weight results in a smoother surface (ESRI, 2010b).

The third technique evaluated is a vector representation of terrain: triangulation. Triangulation creates triangular irregular networks (TINs) by forming triangles between depth measurements. In this study, the Delaunay triangulation technique was evaluated. This technique ensures that no measurements lie within the interior of any of the circumcircles of the triangles in the network, while also maximizing the smallest interior angle of all triangles to avoid long, thin triangles (Soucy and Laurendeau, 1996).

Terrain Characteristics

The terrain characteristics—specifically the slope and curvature—were used to assess the ability of each of the three previously described interpolation techniques to accurately represent heterogeneous terrain. The slope (also referred to as gradient) represents the maximum rate of change of depths, commonly in a three by three cell moving window (Burrough and McDonnell, 1998). A greater slope (units in degrees) corresponds to steeper surfaces. Curvature is the rate of change of the slope. There are three types of curvature: profile, plan, and total curvature. Profile and plan curvatures are the curvatures on a line formed by the intersection of a plane and the terrain surface. The curvature of a line is the reciprocal of the radius of curvature, so a gradually changing curve has a small curvature, while a tight curve has a large curvature (Galant and Wilson, 2000). The total curvature represents the curvature of the surface itself, not the curvature of a line across the surface in a given direction (Galant and Wilson, 2000). A positive total curvature value represents an upwardly convex surface, while a negative total curvature value represents an upwardly concave surface. A curvature value of zero represents an area of constant slope. Both the slope and the total curvature were used to characterize the relationship between interpolation accuracy and terrain characteristics for each interpolation technique, and the sign and strength of these relationships were quantified by calculating the Spearman's rank correlation coefficient (Spearman, 1904).

Data

Kachemak Bay, Alaska, was surveyed in 2008 with high-resolution multibeam swath sonar as part of the NOAA Hydropalooza project. One of the surveys, H11934, spans the entire study area and was downloaded from NGDC's National Ocean Service (NOS) Hydrographic Data Base in Bathymetric Attributed Grid (BAG) format (Brennan et al., 2005). This BAG had a 4 m cell size and was referenced to UTM Zone 5 North and mean lower low water. It was converted to XYZ (eastings, northings, depth) format, and then averaged at a 10-m by 10-m cell size raster. This averaging ensured that every cell in the measured bathymetry raster (Figure 2A) had at least one contributing survey depth measurement, i.e., no values from interpolation. This 10-m cell size raster will be referred to as the measured depth raster. Slope and curvature rasters were derived from this measured depth raster (Figure 2B and C).

Figure 2.

The measured depths (A), slope (B), and curvature (C) used in this study. The measured depth raster, at a 10-m cell size, was the source of all split-sample depths used for interpolation and subsequent accuracy evaluation. The slope and curvature rasters were derived from the measured depth raster in order to analyze the effect of heterogeneous terrain on interpolation accuracy. The pixilation evident in the curvature raster (C) reflects amplification of fine-scale noise in the measured depth raster and original survey.

Accuracy Assessment

There are numerous methods to assess the accuracy of interpolation techniques. One approach is to quantify the differences between the interpolated surface and a higher resolution independent dataset (Chang and Tsai, 1991; Erdogan, 2009; Fisher and Tate, 2006; Robinson and Metternicht, 2005). This is often not possible due to the time and monetary resources needed to collect the independent dataset. Furthermore, it makes little sense to develop a DEM using interpolation and sparse measurements if a higher resolution DEM is available.

More often, cross-validation and split-sample methodologies are used to determine the accuracy of interpolation techniques. The most common form of cross-validation is the “leave one technique” (Erdogan, 2009). This method consists of omitting one depth measurement prior to using an interpolation technique; the difference between the interpolated depth and the omitted measured depth is then calculated. This process is repeated so that all depth measurements are omitted once. The differences between interpolated depths and omitted depths are then aggregated and an indication of the accuracy of the interpolation technique is provided by a global, non-spatial statistic, such as mean absolute error (MAE) or root mean square error (RMSE).

A split-sample methodology uses the same method as cross-validation, except that a split-sample method randomly omits a percentage of depth measurements and calculates the differences between the omitted depths and depths interpolated using the retained depths. The split-sample method is often used to assess changes in the accuracy of an interpolation technique when using various sampling densities (Declerq, 1996; Erdogan, 2009; Smith, Holland, and Longley, 2005). A split-sample method was used in this research to simulate sparse bathymetric datasets due to an interest in interpolating across many cells. The method splits depths randomly at cell sampling densities of 1%, 5%, 10%, 25%, and 50% (Figure 3). The Kachemak Bay measured depth raster (Figure 2A) was the source for all depths used for interpolation and evaluation.

Figure 3.

Flowchart depicting the split-sample method for quantifying interpolation deviations. (A) The original data, averaged to have exactly one depth value per grid cell, are randomly split by a fixed percentage (e.g., 50%) into training (“retained”) and test (“omitted”) data. (B) An interpolation technique (e.g., spline) is applied to the training data to build an interpolated DEM (interpolated cells represented by hollow squares). (C) The interpolated DEM is compared to the test data to quantify the interpolation deviations (represented by filled squares). Steps (A) to (C) are repeated for each interpolation technique at a consistent split-percentage to account for spatial bathymetric variability (i.e., randomness producing different training and test data). The interpolation technique is run iteratively using different split-percentages to evaluate the change in magnitude of interpolation deviations of the chosen interpolation technique at various cell sampling densities.

Implementation of Split-Sample Methodology

The split-sample methodology was implemented 200 times, for each technique, at each of the five cell sampling densities. During each split-sample routine, the retained depths were gridded using the specified interpolation technique, and the resulting interpolated raster was compared, on a cell-by-cell basis, to the measured depth raster to quantify the interpolation deviation. For each cell, the Euclidean distance to the nearest measurement was calculated—measured in map units (meters) but converted to raster cell units by dividing by the cell size (10 m). The interpolation deviations were then statistically compared, on a cell-to-cell basis, to the distances to the nearest measurement, as well as the slope and curvature rasters. The cell comparisons were aggregated, at each of the five cell sampling densities, and used to assess the relationships between these variables. They were also used to develop predictive equations of interpolation uncertainty, using ordinary least squares (OLS) linear regression, which may be applicable in other areas with sparse depth measurements.

Numerous steps were taken to limit bias in the analysis. First, reasonable parameters for each interpolation technique were identified, as inappropriate parameters can cause significant artifacts in the interpolated surface, especially with IDW (Akkala, Devabhaktuni, and Kumar, 2010; Burrough and McDonnell, 1998) and spline interpolation (Cebecauer, Hofierka, and Suri, 2002; Mitasova and Mitas, 1993). Numerous combinations of the interpolation parameters were evaluated using a brute-force methodology. The parameters for each interpolation technique that resulted in the lowest median percent deviation for the entire study area were identified and then used in all subsequent analyses. The parameters for IDW were “Power” of 2 and a “Variable Search Radius” of 8 points. The parameters for spline were regularized “Spline Type,” a “Weight” of 1, and 100 “Number of Points.” The triangulation parameter used in analysis was a linear interpolation “Raster Conversion Method.” Each TIN vector representation of the Kachemak Bay seafloor was converted to a 10-m cell size raster for direct comparison to the other two interpolation techniques. There can be loss of accuracy when converting a TIN to a raster format, but a comprehensive evaluation of such loss was outside of the scope of this research and was assumed to be negligible. In general, the raster better represents the TIN surface as the cell size decreases (ESRI, 2010c).

Depth measurements along the border of the study area were used to guide bathymetric interpolation and minimize edge effects during analysis. The area of analysis was also restricted to a subset of the original study area by 12 cells (120 m) on each side. This ensured that the edge of the area of analysis had a range of distances to the nearest measurement without biasing toward depths along the outermost border of the study area (Figure 4). This subset area used for analysis also ensured slope and curvature values were accurate by having a full three by three window for calculation.

Figure 4.

An example of the distance to the nearest depth measurement raster. The data “buffer” of depth measurements along the border of the study area guided interpolation across the region of analysis outlined by the inner box (map shown at 30-m cell size for visualization purposes).

The deviations of the interpolated depths were quantified as percent deviations by dividing the deviation by the original measured depth. Tsunami modeling requires more accurate, higher resolution DEMs in shallow waters near-shore than in deeper waters off-shore (Titov et al., 2003). Normalizing the deviation by the measured depth balances the importance of interpolation accuracy for both shallow and deep waters. Furthermore, this normalization allows for the predictive equations of interpolation uncertainty developed from the percent deviations to be used irrespective of water depth.

To isolate the effect of terrain on interpolation deviations, the distance to the nearest measurement needed to be similar for every raster cell in the study area. The split-sample method was a random process; therefore, the median distance to the nearest measurement for every cell wasn't exactly the same, but instead showed a Gaussian distribution. For example, after aggregating 200 iterations of the split-sample method at the 1% cell sampling density, the median distance to the nearest measurement for each cell ranged from 4.3 to 5.7 cells, with a median distance of 5 cells. With a similar median distance to the nearest measurement, variability in the interpolation deviations could be attributed to the terrain characteristics. Interpolation median percent deviation maps were created from the aggregation of 200 split-sample routines at each cell sampling density to address the inadequacy of global, non-spatial measures of accuracy by visualizing clusters of high and low accuracy resulting from varying terrain characteristics.

Influence of Cell Sampling Density on Interpolation Deviations

The largest differences in absolute median percent deviation between the three interpolation techniques occurred when using a 1% cell sampling density. The median value from the aggregation of the absolute median percent deviations after 200 split-sample routines when using a 1% cell sampling density was 0.27% for spline, 0.78% for IDW, and 0.39% for triangulation. As the cell sampling density increased, the median distance to the nearest measurement decreased, and the accuracy of all interpolation techniques increased (Figure 5).

Figure 5.

Absolute median percent deviation as a function of cell sampling density for the entire study area. The largest differences in absolute median percent deviation between the three interpolation techniques occurred when using a 1% cell sampling density. As the cell sampling density increased, the median distance to the nearest measurement decreased, and the accuracy of all interpolation techniques increased.

Influence of Terrain on Interpolation Deviations

The interpolation median percent deviation maps created from the aggregation of 200 split-sample routines illustrate the non-stationarity of interpolation deviations from the measured depths, which is the result of heterogeneous terrain with varying slope and curvature (Figure 6). In areas of flat terrain (low slope and low curvature), the deviations for interpolation techniques were small, whereas the interpolation deviations were larger in areas of high slope and curvature. A negative median percent deviation indicated that the interpolated depths were shallower than the measured depths, and a positive median percent deviation indicated that the interpolated depths were deeper than the measured depths. The median percent deviation maps also visually illustrate that spline interpolation is the most accurate of the three evaluated techniques (Figure 6).

Figure 6.

The median percent deviation after 200 IDW (A), spline (B), and triangulation (C) interpolation routines using a random 1% cell sampling density. A negative median percent deviation indicated that the interpolated depths were shallower than the measured depths, and a positive median percent deviation indicated that the interpolated depths were deeper than the measured depths. IDW and triangulation resulted in larger deviations from measured depths, especially in areas of high slope and curvature. Spline resulted in smaller deviations, with its largest deviations in areas of high curvature. See Figure 7 for a more detailed figure on the influence of terrain on interpolation deviations.

All interpolation techniques show a similar pattern of deviations that is strongly associated with the terrain. When visually comparing the median percent deviation maps in Figure 6 to the slope and curvature rasters of the study area shown in Figure 2B and C, it is evident that the interpolation deviations for all techniques are positively correlated with both slope and curvature. The Spearman's rank correlation coefficient between the absolute median percent deviation and the slope using the 1% cell sampling density was 0.46 for IDW, 0.26 for spline, and 0.42 for triangulation, and all correlations were statistically significant with P-values less than 0.001. The Spearman's rank correlation coefficient between the median percent deviation and the curvature at the 1% cell sampling density was 0.44 for IDW, 0.56 for spline, and 0.49 for triangulation, and all correlations were statistically significant with P-values less than 0.001. The influence of terrain slope and curvature on the differences between the interpolated depths and the measured depth raster can be illustrated by extracting a profile across the south-central portion of the original study area (Figure 7D).

Figure 7.

The influence of terrain on the differences between the median interpolated depths and the original measured depth raster after 200 iterations using a random 1% cell sampling density. A profile across the south-central portion of the original study area is represented by the black line overlaying the color images of the measured depths (A), the slope (B), and curvature (C) derived from the measured depths. (D) shows the profile and the differences between the interpolation techniques in representing areas of high slope (dashed arrow) and high curvature (solid arrow). IDW and triangulation produced local maxima that are “pulled” down compared to the measured depths and local minima that are “pulled” up compared to the measured depths, and spline can produce “overshoots” near areas of high curvature.

Predictive Equations of Cell-level Interpolation Uncertainty

Predictive equations of cell-level interpolation uncertainty were derived from the quantification of interpolation deviations from the measured depths. Figure 8 shows the relationship between the absolute median percent deviation and the distance to the nearest measurement at the various cell sampling densities (1%, 5%, 10%, 25%, and 50%) with best-fit OLS regression lines (Table 1). The relationship between these variables diverges from the stable linear trend when the statistical sample size becomes small (less than ~300 counts). At smaller sample sizes, the absolute median percent deviation may be biased by the local terrain effects rather than being representative of the median slope and curvature across the entire study area. Those statistically insignificant results due to low sample size were eliminated prior to deriving the best-fit linear regression equations (Table 1). The y-intercepts of the regression equations are zero as all three techniques are exact interpolators.

Figure 8.

Relationship between the absolute median percent deviation and the distance to the nearest measurement at the 1%, 5%, 10%, 25%, and 50% cell sampling densities shown by the best-fit OLS regression lines for IDW (A), spline (B), and triangulation (C). For each interpolation technique, the smallest cell sampling density (1%) resulted in the highest slope and highest absolute median percent deviation at all interpolation distances (regression line on top of each panel). The slope and absolute median percent deviation decreases at larger cell sampling densities. Also, the length of the regression lines decreases at larger cell sampling densities as the distances to the nearest measurement decreases due to smaller interpolation distances.

Table 1

Regression equations^* for predicting cell-level interpolation percent uncertainty in other areas of sparse depth measurements. The equations predict the cell-level interpolation percent uncertainty (y) as a function of distance to the nearest measurement (x), in cell units, at various cell sampling densities.

Influence of Terrain on Interpolation Deviations

Guo et al. (2010) established that interpolation deviations were correlated with the elevation coefficient of variation of the terrain, but indicated that more research was needed to identify the specific terrain variables, such as slope or curvature, that will have the greatest effect on interpolation accuracy. Research indicates that the deviations for each interpolation technique are positively correlated with both slope and curvature and that there are also differences with the strength of the correlations for each technique. The Spearman's rank correlation coefficients indicate that deviations for all three interpolation techniques are most strongly correlated with curvature. Deviations with spline interpolation are especially correlated with curvature as the spline-interpolated surface exhibits “overshoots” near areas of high curvature due to its minimum curvature algorithm (Figure 7D). The Spearman's rank correlation coefficients also indicate that deviations for all three interpolation techniques are correlated with slope, although spline deviations exhibit a lower correlation coefficient than both IDW and triangulation. Both IDW and triangulation are linear-weighted algorithms and are less accurate near areas of high slope, as local bathymetric maxima and minima cannot be represented unless they have been directly measured. Accordingly, local maxima are “pulled” down and local minima are “pulled” up by surrounding measurements with both IDW and triangulation interpolation techniques (Figure 7D).

Erdogan (2010) modeled interpolation deviations with slope and curvature as the independent variables at various sampling densities with both OLS and a local regression technique—geographically weighted regression (GWR). Erdogan (2010) used these regression techniques as exploratory tools to assess the relationship between interpolation deviations and terrain variability, and to assess the differences between the regression techniques in modeling the non-stationary magnitude of interpolation deviations that result from heterogeneous terrain. These comparative regression analyses were not performed in this study, though this may be worthy of further investigation.

Predictive Equations of Cell-level Interpolation Uncertainty

The predictive equations of cell-level interpolation uncertainty based on the cell sampling density and distance to the nearest measurement (Table 1) all show linear trends in Figure 8, but the slope of the regression line decreases with larger cell sampling densities. The differences in the slope of the line indicate that it is not only the distance to the nearest measurement to an unconstrained cell that determines the accuracy of the interpolated cell, but that the accuracy is also related to the cell sampling density. Even if the distance to the nearest measurement is the same when using different cell sampling densities, the distance to the next nearest measurements that are also used for interpolation will likely be closer at larger cell sampling densities. These closer depths that are also used for interpolation are more similar to the unknown depth to be interpolated according to Tobler's 1^st Law of Geography (Tobler, 1970), resulting in smaller interpolation deviations at larger cell sampling densities.

More specifically, the effect of cell sampling densities is evident when comparing the median percent deviations for a given interpolation technique at a specific distance to the nearest measurement for the various cell sampling densities. For example, the median percent deviation at a distance of one cell to the nearest measurement decreases with increasing cell sampling density. For IDW, the median percent deviation at a one cell distance away from a depth measurement is 0.28%, 0.22%, 0.19%, 0.14%, and 0.1% at 1%, 5%, 10%, 25%, and 50% cell sampling density, respectively. Spline and triangulation exhibit the same decreasing median percent deviation at a distance of one cell to the nearest measurement when increasing the cell sampling density, and all three interpolation techniques also exhibit this decreasing median percent deviation for other distances to the nearest measurement when increasing the cell sampling density.

Accordingly, both the cell sampling density and distance to the nearest measurement were used to develop predictive equations of the interpolation uncertainty on a cell-to-cell basis. The cell sampling density can be calculated for an area of sparse depth measurements as the percentage of cells in the DEM constrained by a depth measurement at a given cell size. A DEM developer could therefore change the cell sampling density for an area with depth measurements by adjusting the cell size. For example, decreasing the cell size would result in a smaller percentage of cells in the DEM constrained by depth measurements, and therefore a smaller cell sampling density. Also, if the cell sampling density is between one of the percentages evaluated in this study, a weighted-average of the two closest linear regression equations provided in Table 1 could be calculated to predict the median percent interpolation uncertainty.

The relationship between the cell sampling density and the absolute median percent deviation shown in Figure 5 provides a global, non-spatial assessment of interpolation accuracy. Such assessments assume a similar distance to the nearest measurement and terrain variability for every cell, which typically is not the case with sparse depth measurements. The inadequacy of these accuracy assessments regarding the typically irregular distribution of measurements, and subsequent varying interpolation distances, is addressed by deriving equations to estimate interpolation uncertainty using distance to the nearest measurement and the cell sampling density. Methods to address the inadequacy of these global, non-spatial assessments of accuracy regarding terrain variability are proposed toward the end of the following subsection.

Limitations

The equations in Table 1 were developed from quantifying the interpolation deviations for only one DEM cell size (10 m). Further research using other cell sizes may reveal the impact of DEM cell size on derived predictive equations of interpolation uncertainty and determine if the equations in Table 1 are more widely applicable across a range of DEM cell sizes. Furthermore, these regression equations are best suited for Kachemak Bay and presumably for areas of similar terrain, as the slope and the curvature were found to be positively correlated with the interpolation deviations. The regression equations provided in Table 1 can be applied to other areas with the caveat that these equations may underestimate interpolation uncertainty in areas of higher slope and curvature than Kachemak Bay, and overestimate interpolation uncertainty in areas of lower slope and curvature (e.g., plains). It would be useful to incorporate terrain characteristics, such as slope and curvature, directly into predictive equations of interpolation uncertainty. However, such efforts fall outside the focus of this study, and terrain characteristics would be presumably unknown in regions with sparse measurements requiring large interpolation distances anyway.

Slope and curvature could indirectly be incorporated into the predictive equations of interpolation uncertainty by implementing the methodology presented in this research on a number of different terrains (e.g., low slope and low curvature; high slope and high curvature; high slope and low curvature) to create a range of uncertainty equations. If there was a densely mapped area nearby with presumably similar terrain to an area requiring interpolation, as expected when considering Tobler's 1^st Law of Geography (Tobler, 1970), then this information could be used to further constrain the predictive uncertainty equation. While potentially useful, such analyses were beyond the scope of the research presented in this manuscript.

The results of this research do suggest that slope and curvature would aid in predicting interpolation uncertainty when interpolating across small distances (e.g., one or two cells) in topographic DEMs. In these cases, the slope and curvature could be reasonably estimated from the interpolated DEM. Future research could advance the results presented in this manuscript, as well as the work of Aguilar et al. (2005), Carlisle (2005), Erdogan (2009), Erdogan (2010), and Guo et al. (2010), by improving the incorporation of terrain characteristics in predicting interpolation uncertainty in both bathymetric and topographic DEMs.