## ABSTRACT

Amante, C.J. and Eakins, B.W., 2016. Accuracy of interpolated bathymetry in digital elevation models. *In*: Brock, J.C.; Gesch, D.B.; Parrish, C.E.; Rogers, J.N., and Wright, C.W. (eds.), *Advances in Topobathymetric Mapping, Models, and Applications*. *Journal of Coastal Research*, Special Issue, No. 76, pp. 123–133. Coconut Creek (Florida), ISSN 0749-0208.

Digital elevation models (DEMs) are used to model numerous coastal processes, including tsunamis, contaminant dispersal, and erosion. In the bathymetric realm, the distance between measurements typically increases farther from shore (*i.e.*, deeper water), such that gridding interpolation to build a bathymetric DEM is often across large distances. This study examined the accuracy of interpolation in bathymetric DEMs using three common interpolation techniques: inverse distance weighting, spline, and triangulation. The goal was to examine the relationship between interpolation accuracy and cell sampling density, distance to the nearest depth measurement, and terrain characteristics. Kachemak Bay, Alaska, was chosen as the study area due to its heterogeneous terrain. A split-sample method was developed to randomly separate depth measurements to be used for interpolation from those used to quantify interpolation accuracy. Results show that the accuracy of the three evaluated interpolation techniques decreases (i) at smaller cell sampling densities, (ii) as the distance to the nearest measurement increases, and (iii) in areas of high slope and curvature. Spline was found to be the most accurate technique, though all techniques have approximately equivalent accuracy at large cell sampling densities and shorter interpolation distances. From these analyses, predictive equations were derived, for each interpolation technique, of the cell-level uncertainty introduced into bathymetric DEMs, as a function of the cell sampling density and interpolation distance. These equations permit the quantification of cell-level interpolation uncertainty in DEMs and, in turn, will aid in propagating that uncertainty into the modeling of coastal processes that utilize DEMs.

## INTRODUCTION

The oceans are truly Earth's last great unknown (Smith and Marks, 2014). Surprisingly, the surfaces of Mars, Venus, and the Earth's moon are mapped at a higher spatial resolution than the seafloor (Smith, 2004), as individual depth measurements of the Earth's ocean can have data gaps of hundreds of kilometers (Smith and Sandwell, 1997). Bathymetric digital elevation models (DEMs) are continuous representations of the seafloor that are derived from depth measurements and are commonly stored in a raster data format comprised of a matrix of square cells, with each cell representing the average depth of the area contained within that cell.

The National Oceanic and Atmospheric Administration's (NOAA) National Geophysical Data Center (NGDC) develops such bathymetric DEMs, as well as bathymetric models that are integrated with topography in the coastal zone, to support the modeling of tsunami propagation and coastal inundation (Eakins and Taylor, 2010). In order for the integrated bathymetric-topographic DEMs to retain high spatial resolution in the coastal zone to support detailed inundation mapping, interpolation across large distances is typically required between sparse bathymetric measurements. Such extreme interpolation fills data gaps in the models that might cause instabilities in tsunami modeling runs (*i.e.*, avoids null depth values).

The accuracy of modeling coastal processes, such as tsunami propagation and inundation, is partly dependent on the accuracy of the DEM, which in turn is dependent on the accuracy of the depth measurements. As a result, uncertainty in the depth measurements propagates into uncertainty in the DEM, and into the modeling of such processes. A detailed description of depth measurement uncertainty is provided in Hare, Eakins, and Amante (2011). However, the low spatial resolution of many bathymetric surveys, especially those predating the use of modern survey technologies (*e.g.*, multibeam sonar), requires interpolation to estimate areas of unknown depths to develop the continuous representation of the seafloor. Consequently, the accuracy of modeling such coastal processes is also partly dependent on the accuracy of the interpolation technique.

The accuracy of an interpolation technique itself refers to the agreement of interpolated values to known or accepted values. However, interpolation techniques are typically used in areas without measured depths; consequently, the interpolation accuracy cannot be rigorously quantified. The lack of knowledge about interpolation deviations from known or accepted values within such data gaps represents the uncertainty introduced into DEMs by interpolation techniques (Wechsler, 2007). This study quantified the accuracy of various interpolation techniques, using measured depths that are withheld from interpolation and subsequently used to quantify the interpolation accuracy. These depth measurements were considered the known or accepted values, despite any uncertainty of the measurements themselves, as previously mentioned. Predictive equations of cell-level interpolation uncertainty were then derived from the quantification of interpolation deviations from the measured depths.

### 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).

## METHODS

This study evaluated the accuracy of three common deterministic, local, exact interpolation techniques—inverse distance weighting (IDW), spline, and triangulation—in developing bathymetric DEMs of Kachemak Bay, Alaska. The overarching goal of the study was to examine the relationship between interpolation accuracy and the (i) cell sampling density, (ii) distance to the nearest depth measurement, and (iii) terrain characteristics, specifically the slope and curvature. These three deterministic techniques were analyzed in this study because unlike geostatistical techniques such as kriging, deterministic techniques do not provide an assessment of uncertainty along with the predicted values (Li and Heap, 2008). The accuracy of the three interpolation techniques was analyzed in this study using the Environmental Systems Research Institute (ESRI) ArcGIS 10.0 software (ESRI, 2010a).

### 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).

### 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.

### 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.

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.

## RESULTS

Results show that the accuracy of the three evaluated interpolation techniques decreases (i) at smaller cell sampling densities, (ii) as the distance to the nearest measurement increases, and (iii) in areas of high slope and curvature. Spline was found to be the most accurate interpolation technique, though all techniques have approximately equivalent accuracy at large cell sampling densities and shorter interpolation distances. From these analyses, predictive equations were derived, for each interpolation technique, of the cell-level uncertainty introduced into bathymetric DEMs, as a function of the cell sampling density and interpolation distance.

### 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).

### 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).

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).

### 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.

## DISCUSSION

These results support previous studies that investigated the accuracy of interpolation techniques in relationship to the sampling density (Aguilar *et al.*, 2005; Anderson, Thompson, and Austin, 2005; Chaplot *et al.*, 2006; Erdogan, 2009, 2010; Guo *et al.*, 2010; MacEachren and Davidson, 1987) and terrain characteristics (Aguilar *et al.*, 2005; Carlisle, 2005; Erdogan, 2009, 2010; Guo *et al.*, 2010). Previous research focused on topographic DEMs, and this study indicates that the accuracy of interpolating bathymetric DEMs produces similar results in spite of the different terrain characteristics. This study advances previous research by implementing the split-sample method 200 times and aggregating the interpolation deviations to show a clear relationship between interpolation deviations and terrain characteristics that was not biased by the distance to the nearest measurement. Furthermore, use of a data “buffer” around the study area, as shown in Figure 4, minimized edge effects in the interpolated DEMs. The buffer also minimized the influence of border measurements when deriving predictive equations of cell-level interpolation uncertainty.

### 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.

## CONCLUSIONS

This study evaluated the accuracy of three common interpolation techniques (IDW, spline, and triangulation) in developing bathymetric DEMs of Kachemak Bay, Alaska. Spline was found to be the most accurate interpolation technique, though all techniques have approximately equivalent accuracy at large cell sampling densities and shorter interpolation distances. Furthermore, the relationship between interpolation deviations from measured depths and the distance to the nearest depth measurement for each interpolation technique with various cell sampling densities were quantified. From this quantification, regression equations were derived at each cell sampling density that can be used to predict, on a cell by cell basis, the uncertainty introduced into DEMs by these interpolation techniques.

Kachemak Bay, Alaska, was chosen as the study area because it has heterogeneous terrain (*i.e.*, a wide range of spatially varying slope and curvature values). Similar interpolation accuracy analyses in other types of terrain would derive regression equations—based on the relationship between interpolation deviations and distance to the nearest depth measurement for various cell sampling densities—that are *roughly* consistent with these results for Kachemak Bay bathymetry. This research has indicated that interpolation deviations are positively correlated with both slope and curvature, and future research efforts should implement the robust and easily reproducible methodology described in this manuscript to derive more refined predictive equations of interpolation uncertainty for terrain that differs from Kachemak Bay. Future research should also investigate the impact of DEM cell size on derived predictive equations of interpolation uncertainty.

The modeling of numerous coastal processes, such as tsunami propagation and inundation, contaminant dispersal, and erosion, contain inherent uncertainty in their results that originates from many different sources. DEMs, which are used to model numerous coastal processes, are often one of those sources of uncertainty. The ability to quantify DEM uncertainty, which includes uncertainty from the source data, morphologic change, and from the uncertainty introduced by the interpolation technique described in this manuscript, will support efforts to quantify the uncertainty in coastal research results that utilize these models of the Earth's surface.

## ACKNOWLEDGMENTS

The research in this manuscript was part of Christopher Amante's Master's thesis at the University of Colorado Boulder (see Amante, 2012) and was supported by funding from the National Tsunami Hazard Mitigation Program (NTHMP). The authors would like to thank many at NOAA NGDC, especially Matthew Love, Lisa Taylor, Sue McLean, Chris Fox, Kelly Carignan, Pamela Grothe, Dorothy Friday, and Elliot Lim for all of their support. Christopher Amante would also like to thank Waleed Abdalati and Barbara Buttenfield for their guidance throughout his Master's thesis. The authors also thank Mike Sutherland at NOAA NGDC and the two anonymous reviewers for their valuable feedback that significantly improved the quality of this manuscript.