EN中文

Journal of Geographical Sciences >

2015 , Vol. 25 >Issue 2: 196 - 210

DOI: https://doi.org/10.1007/s11442-015-1162-2

Orginal Article

Reconstructing pre-erosion topography using spatial interpolation techniques: A validation-based approach

  • Rafaello BERGONSE ,
  • Eusébio REIS
Expand
  • Center for Geographical Studies, University of Lisbon, Portugal

Author: Rafaello Bergonse, E-mail:

Received date: 2013-06-05

  Accepted date: 2014-04-20

  Online published: 2015-02-15

Copyright

Journal of Geographical Sciences, All Rights Reserved
Fold

Abstract

Understanding the topographic context preceding the development of erosive landforms is of major relevance in geomorphic research, as topography is an important factor on both water and mass movement-related erosion, and knowledge of the original surface is a condition for quantifying the volume of eroded material. Although any reconstruction implies assuming that the resulting surface reflects the original topography, past works have been dominated by linear interpolation methods, incapable of generating curved surfaces in areas with no data or values outside the range of variation of inputs. In spite of these limitations, impossibility of validation has led to the assumption of surface representativity never being challenged. In this paper, a validation-based method is applied in order to define the optimal interpolation technique for reconstructing pre-erosion topography in a given study area. In spite of the absence of the original surface, different techniques can be nonetheless evaluated by quantifying their capacity to reproduce known topography in unincised locations within the same geomorphic contexts of existing erosive landforms. A linear method (Triangulated Irregular Network, TIN) and 23 parameterizations of three distinct Spline interpolation techniques were compared using 50 test areas in a context of research on large gully dynamics in the South of Portugal. Results show that almost all Spline methods produced smaller errors than the TIN, and that the latter produced a mean absolute error 61.4% higher than the best Spline method, clearly establishing both the better adjustment of Splines to the geomorphic context considered and the limitations of linear approaches. The proposed method can easily be applied to different interpolation techniques and topographic contexts, enabling better calculations of eroded volumes and denudation rates as well as the investigation of controls by antecedent topographic form over erosive processes.

Key words： pre-erosion topography; surface reconstruction; spatial interpolation; spline interpolation; triangulated irregular networks; erosive landforms; gully erosion

Cite this article

Rafaello BERGONSE , Eusébio REIS . Reconstructing pre-erosion topography using spatial interpolation techniques: A validation-based approach[J]. Journal of Geographical Sciences, 2015 , 25(2) : 196 -210 . DOI: 10.1007/s11442-015-1162-2

1 Introduction

Topographic form is both a cause and a consequence of erosive processes. Concerning water erosion, slope (i.e. the first derivative of altitude) controls the velocity of flow, and hence the shear stress it exerts. Relations between slope and shear stress are explicit in the shear stress equation (e.g. Poesen et al., 2003; Yao et al., 2008):
where ρ is water density (g·cm-3), g is gravitational acceleration (cm·s-2), R is the hydraulic radius①(①This formulation assumes channeled flow. For unconcentrated runoff, mean flow depth can be used.), and S is the sine of the soil surface gradient. τ is expressed in Pa.
Topography also exerts control over sub-surface flow, promoting concentration in plane concavities (Dunne et al., 1975; Anderson and Burt, 1978). The regolith in these areas will have a relatively high water content, be saturated for more time, and produce relatively more runoff. The increased water content will at the same time reduce its shear strength (Wood et al., 2001).
The control exerted by topographic form over water erosion has justified the numerous slope and drainage area combinations associated in the literature to the location of incipient gullies (e.g. Zucca et al., 2006; Kakembo et al., 2009), the magnitude of surface flow (e.g. stream power index, Daba et al., 2003) or the tendency for elevated soil water content (e.g. wetness index, Bou Kheir et al., 2007). Simple properties such as profile and plan curvature have been used, for example, to predict lateral expansion in large gullies (Martínez- Casasnovas et al., 2009) or gully location and dimension (Bou Kheir et al., 2007). More complex indicators combining different properties have also been applied, such as the morphometric slope index proposed by Buccolini et al. (2012).
In parallel to water erosion, topographic properties are strongly related to the occurrence of mass movements. Slope is a major factor on the shear stress exerted over hillslope materials, being inversely proportional to their shear strength (by reducing stress normal to the shear plane, Summerfield, 1991).
At the same time, the loss of shear strength associated to a higher water content (e.g. Wood et al., 2001) will lower the resistance of the regolith against forces promoting collapse, with these forces being increased by the mass of water itself.
Apart from its significance as a control on erosion, characterizing the topography preceding the development of erosive landforms is a condition for estimating the volume of the eroded mass. If ages can be associated to existing surfaces, mean rates of denudation can be derived from these volumes (e.g. Alexander et al., 2008).
Due to their enormous variability in terms of dimension, from rills to continuous multi-hectare gully complexes, erosive landforms frequently have expression on the topographic information available (e.g. contour maps, DEMs, LIDAR points). In these cases, characterizing surface form in eroded areas from this information would imply characterizing the erosive features rather than their antecedent topography. This makes clear that any articulation between pre-erosion surfaces and large erosive landforms (i.e. either for volume estimation or investigating causal relations) must therefore begin by a reconstruction of those same surfaces using some method of spatial interpolation, i.e. using knowledge of the currently unincised topography to estimate the surface corresponding to incised areas.
Published research includes several reconstructions of pre-erosion surfaces with different purposes, of which a synthesis is presented in Table 1.
Table 1 Contexts and methodologies of some published pre-erosion surface reconstructions
Author Purpose of reconstruction Interpolation method
Wells and Gutiérrez (1982) Estimation of eroded volumes and combination of results with current mean erosion rates in order to date badland initiation Undefined
Daba et al. (2003) Estimation of eroded volume in a large gully system for two different dates; comparison of results in order to quantify temporal evolution Undefined1
Alexander et al. (2008) Understanding the geomorphic evolution of a badland site from a set of remnant surfaces, and estimating rates of denudation Linear interpolation (Triangulated Irregular Networks)
Perroy et al. (2010) Estimating volumetric soil loss from a set of gully channels Linear interpolation
(grid based)
Buccolini et al. (2012) Estimating the volume eroded by a set of gully systems (calanchi), and relating pre-erosion topography to gully system properties Linear interpolation (manual2)

1 Authors used the SCOP (Stuttgard Contour Program) software to interpolate surfaces, but the specific method is not identified.2 Straight contour lines were drawn connecting points with the same height on both sides of the watershed of each calanco.

Any reconstruction implies assuming that the resulting surface represents the original one to the point that differences between both will not bear significantly on research results. Although Table 1 shows that most methodologies found have been based on linear interpolation, this method seems to be particularly inadequate, not only because if one regards real hillslopes it is obvious that some degree of curvature (profile, planform or both) is extremely frequent, but also because curvature is strongly associated to erosive features, as the previous considerations have made clear. As linear interpolation methods, regardless of specific algorithm, always assume linear variation between a set of known points, it is obvious that they cannot generate curved surfaces in areas without data (i.e. incised areas), a limitation that will be most clear when erosive features have an areal rather that a linear configuration, as gully systems frequently do (e.g. Wells et al., 1991; DeRose et al., 1998; Betts et al., 2003; Parkner et al., 2006; Buccolini and Coco, 2010).
Another shortcoming of linear methods is that resulting values are constrained by the range of variation of known ones, implying an incapacity to reproduce valley bottoms that are lower and hilltops that are higher than, respectively, the minimum and maximum known altitudes (e.g. contours). As such, in a context of volume calculation, a linear approach will result in underestimation if the original surface was convex and overestimation if concave, with error increasing with the degree of original curvature and the size of the erosive feature under analysis.
Although it seems clear that linear methods should be inadequate for reconstructing hillslope surfaces, their use has not yet been contested, a fact certainly stemming from the lack of an original surface against which to validate results.
Against this theoretical background, we propose that in spite of the absence of original topography, interpolation methods can be nonetheless evaluated by quantifying their capacity to reproduce known surfaces in unincised locations within the same topographic context of existing erosive landforms. This implies the sole assumption that the measured error, obtained by comparing interpolation results to reality, will be equivalent to that which would be obtained for the eroded areas had the original surface still existed.
Adopting this approach, this paper has three objectives:
(1) To apply a new methodology with the purpose of defining the optimal spatial interpolation method for reconstructing pre-erosion topography in a given geomorphic context;
(2) To evaluate the suitability of the linear interpolation approach (as the method most frequently adopted in the literature) for reconstructing pre-erosion topography;
(3) To reconstruct the topographic surface preceding a set of 90 large gullies and gully complexes presently evolving in two tributary basins of the Lower Tagus River, South Portugal.
The gully dataset used was obtained through interpretation of digital ortophotos produced in 2004 (Portuguese Geographic Institute/General Direction of Forest Resources; 0.5 m resolution) and subsequent field validation (75.6% of all features were directly verified). Gully outlines were vectorized into a polygon shapefile using ArcGIS 9.1.
The area of each studied feature varies from a minimum of 196 m2 to a maximum of 33152 m2, with an average of 2692 m2 and square deviation of 4973 m2.

2 Study area

The study area is constituted by two small basins on the left margin of the Lower Tagus River (Figure 1). The Ulme river basin extends for 138.4 km2 and the Vale do vasal velho river basin occupies just about 12.9 km2. Climate is temperate with a warm and dry summer (Csa in the Köppen system), and annual rainfall varying between 600 and 800 mm. Lithologically, both basins are composed of variably consolidated Tertiary clastic formations (clays, silts, sandstones) topped by extensive conglomeratic deposits, as well as Quaternary alluvium, terraces and deposits. The conglomerates define a very regular interfluve that has been strongly incised by the drainage network during the Quaternary, leading to relatively steep hillslopes in both studied basins. For the Ulme basin, a set of four profiles defined over 1:25000 contour lines at equal intervals along the river’s longitudinal profile allowed to define an average hillslope inclination of 11.1º, with similar topography characterizing the neighbouring Vale do casal velho basin.
Figure 1 The two study basins in the context of the lower Tagus. The set of 90 gullies and gully complexes subjected to pre-erosion surface reconstruction is identified
In contrast to the unincised flat valley bottoms, resulting from Holocene aggradation, hillslopes on the two basins are affected by numerous large gullies and gully complexes, the latter being herein defined as ‘systems of gully channels separated by interfluves that are themselves degraded in relation to the topographic surface exterior to the system’ (Bergonse and Reis, 2011). These systems are sometimes more than 20 m deep and extend through several hectares. Although their bottom sectors show no signs of activity, being mostly colonized by vegetation, walls and headcuts are very steep and unvegetated, currently retreating through mass movements. An example is shown in Figure 2.
Figure 2 A large gully complex in the Ulme river basin. Walls show signs of active retreat (note deposits at the base and Quercus suber with root system completely exposed), contrasting with bottom colonized by vegetation

3 Methodology

The adopted methodology followed the five generic stages outlined in Figure 3.
It begun (1) with the selection of a set of interpolation methods and respective parameters, to be applied and critically compared. These methods were selected either because their characteristics made them suitable for topographic modelling, or because of their frequent usage in the literature. In order to assess and compare their capacity to model topography no longer in place, a set of test areas were selected throughout the study basins in a semi-random way (2), with size and configuration similar to the gully systems under study, but occupying unincised locations (i.e. locations where the actual topography is known). All selected methods were then successively used to interpolate the surface corresponding to these test areas (3). The model surfaces were validated by comparison with the actual topography, extracted from contour maps at a 1:10000 scale (4). In a subsequent and final stage (5), we strived to optimize the methods yielding the best results by applying them once more, but introducing experimental changes in the input parameters.
Figure 3 A schematic outline of the adopted methodology
This final stage ensured that an optimal method was ultimately selected. All stages are described in detail below.

3.1 Preliminary method selection

Numerous spatial interpolation methods are available in different GIS software packages, implying different sets of assumptions about the modelled phenomena. Associated algorithms normally include user-defined parameters, thus extending their advantages and applicability to a greater number of situations.
The methodology proposed begins with a preliminary analysis of the different methods, leading to the selection of one or more that are clearly more suitable for the purpose at hand. It is fundamental to consider that there is not a universal ‘best interpolation method’, each one presenting advantages and disadvantages in different situations (Weibel and Heller, 1991).
A brief synthesis of the most common methods is presented in Table 2. All produce grid-based information, although in the case of linear interpolation a (vector-based) TIN can be first created through Delaunay triangulation and then converted to grid. Only exact methods were considered, i.e. methods producing surfaces in which the values initially known are maintained.
Setting apart inexact interpolators (e.g. polynomial interpolation), the set of remaining methods can be differentiated according to three criteria: (1) the smoothing effect; (2) the proximity effect, and (3) the existence of geostatistical assumptions (Eastman, 2006; Heywood et al., 2006; Hengl and Evans, 2009).
The smoothing effect refers to the adjustment between resulting values and the original data, and has expression in two aspects: on the one hand, the difference between the original values and the generated ones. In this sense, smoothing determines the position of each method across a spectrum having in one extreme an exact result (all original values are maintained, i.e. no smoothing) and on the other a very inexact one (resulting values are very different from the original ones, i.e. elevated smoothing) (Neteler and Mitasova, 2008).
On the other hand, smoothing refers to the differences between the range of variation of resulting values in relation to that of the original ones. No smoothing results in both amplitudes being the same (e.g. linear interpolation), whereas methods with smoothing can produce values beyond the amplitude of the original values (e.g. Splines).
Table 2 General properties of the most common commercially available exact spatial interpolation methods
Method General features Smoothing Proximity Geostatistical assumptions
Linear interpolation May be based on a previous Delaunay triangulation, with the value for each cell being defined by the linear surface of the triangle it overlays (e.g. Surfer 101, ArcGIS 9.1). In other cases, estimations are obtained simply as a function of the nearest known values and the respective distances (e.g. IDRISI Andes2: Eastman, 2006) None Local No
Inverse Distance Weighted Interpolated values are a function of the values of the nearest points (quantity is user-defined), with the weight of each in the result being a function of distance. None Local to Global No
Splines Generated surface results from fitting a polynomial to a quantity of user-defined known values, subjected to two constraints: (1) surface passes exactly through the known data points; (2) curvature of generated surface is minimized. Has problems representing discrete transitions (e.g. limits of flood plains, slope breaks), sometimes ‘overshooting’ the true surface (Hengl and Evans, 2009). Elevated Local to Global No
Topo to Raster Similar to Spline, but modified in order to produce a hydrologically correct surface and incorporate slope breaks. Conceived to use points, lines and polygons as input. Elevated Local to Global3 No
Ordinary Kriging Based on preliminary analysis and statistical modelling of the variation of differences between all known values with spatial distance and/or direction. For each location, the functions thus defined are used to estimate values from surrounding data points of known value. Medium Local to Global Yes
Natural neighbour Based on the construction of a network of Voronoy polygons incorporating all known data points. Each point to be estimated is inserted on the network, and the latter is modified in order to incorporate it. Each estimated value is the average of all known surrounding points of known value, weighted by the proportion of the new Voronoy polygon overlaying each of the initial polygons. None Local No

1 Golden Software;2 Clark Labs;3 In its ArcGIS 9.1 implementation, the algorithm uses a maximum of four input points, thus being local.

The proximity effect describes the area around each value to be estimated where all known values influence the result. It varies between a global extreme (all known values have influence) and a local one (only the more immediate values have influence).
Finally, the existence of geostatistical assumptions differentiates methods in which a preliminary analysis of spatial correlation between known values determines the importance of each one on the estimation process (e.g. ordinary Kriging) from the so-called ‘deterministic’ methods, in which no assumptions are made and the same mathematical formulation is indiscriminately used across all the area to be estimated.
For the present case, Spline interpolation was selected because it allows a variable degree of smoothing, allowing for realistic representations of the topography (Hengl and Evans, 2009). On the contrary, Inverse Distance Weighted and Natural Neighbour interpolations were eliminated for being incapable of generating curved surfaces in areas with no data, and Kriging for frequently exaggerating terrain smoothness and being sensible to values statistically different from the rest of the known population (i.e. hotspots), which makes it little suited to topographic modelling (op.cit.).
The Topo to Raster algorithm, which constitutes an adaptation of the ANUDEM method (Hutchinson, 1988; 1989) implemented in ArcGIS 9.1 (ESRI) was also selected. Although it includes constraints to the interpolation process designed to force a hydrologically correct topography and include abrupt changes in the land surface (e.g. ridges and streams), it is in essence a discretized thin-plate Spline technique (ESRI, 2005), being considered by several authors as especially adequate for modelling topography with morphometric purposes (Hengl and Evans, 2009; Reuter and Nelson, 2009). It was used with all input parameters set as the default in ArcGIS 9.1.
Finally, linear interpolation was also selected in order for its results to be compared with those of other methods, in accordance with the objectives defined in the Introduction.
It is important to mention at this point that the consideration of available interpolation methods included only those available in commercial GIS software packages, i.e. those methods which application to a given study area does not require programming knowledge on the part of the user. Apart from the commercially available methods, the HASM (High Accuracy Surface Modelling) method developed by Yue et al. (2010) should be mentioned, as it allowed its authors to reproduce known values with a significantly higher accuracy than the IDW, Spline, TIN and Kriging methods.
All interpolations were performed using ArcGIS 9.1. This software package includes two main types of Spline methods: Regularized and Tension. Although both include the same two constraints (see Table 1), the Regularized method produces smoother surfaces that vary more gradually. The user-defined parameter weight (w) describes the weight attached to the third derivative of the resulting surface during the curvature minimization procedure. Higher w values will produce a smoother surface (i.e. less close to the range of variation of input values).
In opposition, the Tension method produces surfaces that are less smooth, but closer to the range of variation of known values. The curvature minimization procedure incorporates first derivative terms, whose attached weight is described by the weight parameter. Higher values will produce a coarser surface closely conforming to the input points (ESRI, 2005).
The different methods adopted and the corresponding parameterizations (in a total of 13) are presented in Table 3.
The Topo to Raster method was applied using separately as inputs points and contour lines. This was done not only to evaluate differences in results, but also to explore the capacity of this algorithm to run on contour data, a distinguishing feature in face of all other methods discussed.
Finally, each of the two Spline types was run with five distinct values of w (Table 3), selected taking into consideration the reference values given in the software documentation (ESRI, 2005).
Table 3 The interpolation methods and parameterizations adopted
Method (different parameter sets used) Parameters
Linear interpolation (1) Obtained through triangulation and conversion of a TIN model (points as input)
Topo to Raster (2) Two parameterizations: points as input and contours as input. Further parameters were set as default.
Spline (10) Spline Regularized: w = 0; 0.001; 0.01; 0.1; 0.5
Spline Tension: w = 0, 1, 4, 7, 10.

3.2 Generation of test areas

A set of 50 points was randomly generated using the Create Random Points tool (ArcGIS 9.1). Some of these were then relocated manually across the shortest possible distance in order not to be situated over flat valley bottoms or plateau areas (where no erosive features occur). The same was done for points generated in the unincised headmost sectors of both study basins. Although these changes necessarily compromised the initial randomness of the location of the original point dataset, they centred attention on the topographic contexts where gullies and gully complexes effectively occur, i.e. on the hillslopes.
In order to construct test areas with similar size and configuration to the erosive features under study (as opposed to ‘unnatural’ shapes as circles or rectangles), the 90 photo-interpreted polygons were ordered with respect to the area occupied by each and analysed in order to define the features over 33% and 66% of the distribution. Copies of these (with 854.7 m2 and 1873.5 m2, respectively) were then alternately positioned over each of the 50 points and rotated so as to be oriented along the direction of maximum slope, therefore reproducing the characteristic position of gullies and gully systems in the study area. Two features were selected as models so as to better represent the existing variability in term of shape.

3.3 Surface reconstruction for test areas

Contours at a 1:10000 scale were adopted to represent real topography. In order to generate empty areas for reconstruction, the contour segments coincident with the 50 test area polygons were removed using the Simmetric Difference tool in ARCGis 9.1.
All remaining contours were converted to points using the Features to Points tool, in order to be used as input for the different interpolations.
The set of removed contour segments coincident with the test areas was also converted into points. This point dataset, comprising a total of 2344 points, was later used for validating reconstruction results.
Given the enormous quantity of points generated by the contour conversion process (3084979 in both basins), a slight simplification of the initial lines was performed using the ArcGIS 9.1 simplify line tool. This was done with the purpose of reducing data redundancy as well as computational demands. The simplification algorithm was constrained in order that the distance between simplified and original contours was never greater than 0.25 m, an insignificant variation in comparison to the 2 m resolution adopted for all interpolation results. In order to further minimize computational demands, seven contiguous zones were defined in order to contain all 50 test areas, thus allowing interpolations to be run separately for each zone. Care was taken to place zone limits away from test area polygons, thus ensuring that enough data points were available for interpolation (12 points were defined as input for each interpolated cell following the program defaults except for Topo to Raster, which only uses four points). The seven zones are shown in Figure 4, together with the 50 test areas.
Figure 4 Zones and test areas (represented as points). The limits of the two studied basins are represented with a dashed line
Finally, all the interpolation methods selected (see 3.1 and Table 3) were applied. All used points as input, with the sole exception of Topo to Raster, which was also executed using the simplified contours.
In the case of linear interpolation, a TIN was created and then converted to raster (2 m resolution) by linearly estimating the altitude of each cell from the vertices of the triangle it overlays.

3.4 Validation

Validation was performed by overlaying the set of 2344 points extracted from the original contour data and corresponding to the test areas on the generated surfaces, and extracting the values of the corresponding cells in each to a table using the Extract Values to Points tool (ArcGIS 9.1). These tables were then exported from ArcGiS into Microsoft Excel XP and used to compare real values to interpolated values, and thus measure the absolute error produced by the different interpolation methods. Results are presented in Table 4.
Values show that the Topo to Raster method running over contours produced the lowest Mean Absolute Error (0.752 m), with the second lowest value (higher only by 1.5 cm) being obtained with the Spline Regularized method and w = 0.01.

3.5 Parameter optimization

Setting these two best methods/parameterizations as reference, a second set of interpolations was undertaken with the purpose of optimizing the values of user-defined entry parameters and thus assessing to which degree it was possible to further reduce the errors initially obtained.
Table 4 Characteristics of the distributions of absolute error obtained for each interpolation method (i.e. square root of the square of the difference between real and interpolated values). Methods are ordered by ascending mean absolute error (MAE). Spline Reg and Spline Ten respectively identify the Regularized and Tension methods; w = weight parameter; P50 and P80 are the 50th and 80th percentiles; SD - standard deviation. All values are in metres.
Method MAE Min Max P50 P80 SD
Topo to Raster (contours) 0.752 0.000 5.399 0.440 1.203 0.872
Spline Reg w=0.01 0.767 0.000 5.001 0.463 1.176 0.889
Spline Reg w=0.1 0.771 0.000 5.395 0.470 1.218 0.913
Spline Reg w=0.001 0.810 0.000 5.338 0.493 1.239 0.904
Spline Reg w=0.5 0.813 0.000 5.827 0.473 1.242 1.000
Spline Reg w=0 0.834 0.000 5.456 0.526 1.288 0.912
Spline Ten w=1 0.887 0.000 5.375 0.540 1.496 0.939
Spline Ten w=4 0.954 0.000 5.354 0.587 1.654 0.976
Spline Ten w=7 0.998 0.000 5.349 0.619 1.758 1.004
Spline Ten w=10 1.034 0.000 5.350 0.639 1.828 1.028
Linear 1.214 0.000 6.908 0.815 1.962 1.239
Topo to Raster (points) 1.589 0.000 12.773 1.094 2.383 1.795
Spline Ten w=0 3.463 0.001 72.642 1.084 3.494 7.500
For the Regularized Spline method, a set of six additional values were selected for the parameter w. These were defined at equal intervals between the best initial value (0.01) and the two closest values initially tested (0.001 and 0.01), and are presented in Table 5.
For the Topo to Raster method, the roughness penalty parameter R (initially set as 0, as is the default in ArcGIS 9.1) was defined as 0.1, 0.2, 0.3, 0.4 and 0.5, thus varying between the extreme reference values recommended in the software documentation (0-0.5). This parameter describes the integrated squared second derivative as a measure of roughness, minimized by the algorithm during interpolation (ESRI, 2005). All remaining user-defined input parameters were set as default (e.g. maximum number of iterations, discrete error factor).
Table 5 Characteristics of the distributions of absolute error obtained during parameter optimization. Methods are ordered by ascending mean absolute error (MAE).; w = weight parameter in the regularized (Reg) Spline method; R = roughness penalty in the Topo to Raster method; P50 and P80 are the 50th and 80th percentiles; SD - standard deviation. All values are in metres.
Method Mean Min Max P50 P80 SD
Spline Reg w=0.033 0.758 0.001 5.206 0.458 1.160 0.890
Spline Reg w=0.055 0.762 0.001 5.291 0.468 1.191 0.896
Spline Reg w=0.078 0.766 0.000 5.349 0.466 1.190 0.905
Spline Reg w=0.008 0.771 0.000 5.042 0.463 1.171 0.890
Spline Reg w=0.006 0.776 0.000 5.095 0.466 1.175 0.892
Spline Reg w=0.003 0.791 0.000 5.210 0.480 1.191 0.897
Topo to Raster R=0.4 0.825 0.000 5.499 0.496 1.318 0.922
Topo to Raster R=0.3 0.871 0.000 4.920 0.547 1.442 0.923
Topo to Raster R=0.2 0.938 0.000 5.585 0.618 1.519 0.980
Topo to Raster R=0.1 0.987 0.000 5.634 0.664 1.599 1.007
Topo to Raster R=0.5 1.020 0.000 5.716 0.701 1.635 1.029
Validation results for the resulting surfaces, as well as the parameterizations adopted, are shown in Table 5.
The best result was obtained using the Regularized Spline method with w = 0.033. Despite its proximity to the lowest MAE initially obtained, this earlier value is still higher by 0.006 m, thus defining Topo to Raster with the default parameterization (R = 0) and contour lines as input as the optimal spatial interpolation method for reconstructing pre-erosion topography in the study area (see Table 4).

3.6 Final reconstruction of pre-erosion topography

Following the same methodology already described in 3.3 and in accordance with the objectives stated, the interpolation method defined as optimal was used to reconstruct the topography prior to the formation of the 90 gully systems considered. The segments corresponding to gullied areas were removed from the contour data describing present topography (using the Simmetric Difference tool), and the gaps were interpolated using the Topo to Raster method with the Roughness Penalty parameter R set as 0 in ArcGIS 9.1. In order to better illustrate the result, contours with 5 m equidistance (i.e. same as the equidistance on the original 1:10000 information) were automatically generated with the Contour tool in ArcGIS 9.1, being contrasted to the original topography in Figure 5. The three examples presented show these contours to be strikingly similar in terms of curvature and configuration to the original surrounding topography.
Figure 5 Three examples of surface reconstructions using the optimal interpolation method and parameterization (Topo to Raster with roughness penalty R = 0 and contours as input). (a), (c), (e) - original topography; (b), (d), (f) - reconstructed surfaces

4 Discussion

Tables 4 and 5 show that the different interpolation methods produced contrasting results, with Mean Absolute Errors varying from 0.752 to 3.463 m. Generally, all the remaining distribution parameters represented (to the exception of the minimum error) increase with MAE, showing that higher means are accompanied by a higher dispersion of error values. This increase with the mean is not so clear for the maximum error obtained in each distribution, but is plain for P50, P80 and SD values.
Considering the Spline parameterizations used, the Regularized method consistently produced lower errors than the Tension method, marking it as clearly more adequate to the context of the study area. In fact, all the parameterizations adopted for the Regularized method produced better results than all other methods/parameterizations, with the sole exception of Topo to Raster with roughness penalty R = 0 and running over contours, and this with a MAE lower by only 1.5 cm.
It is worthy of notice that the Topo to Raster method produced much better results using contours as input than using points, showing that for point only datasets (e.g. GPS or LiDAR), the Regularized Spline technique should be preferred. This is in accordance with the original purpose of the Topo to Raster procedure, as it is the only ArcGIS interpolator designed to work intelligently with contour data (ESRI, 2005).
Linear interpolation, herein performed using a preliminary Delaunay triangulation, was clearly shown by the results to be a much weaker method than Splines for reconstructing topography in the geomorphic context considered, with a MAE 61.4% higher than the lowest value obtained, and the third highest among all the 24 methods/parameterizations tested. As curvature is a very frequent quality of natural relief, it is expectable that similar results will be obtained in the majority of hillslope contexts, in accordance with the theoretical considerations made in the Introduction.
The variability in the results obtained underlines the importance of comparing both different spatial interpolation methods and different input parameter values for the same method, as small changes in one single input may significantly influence results. Examples are the contrasting MAE values obtained with Topo to Raster and Roughness Penalty R = 0 (0.752m), and the same method with R = 0.1 (0.987), or Tension Spline with weight w = 1 (0.887 m) and w = 0 (3.463 m).
Finally, it should be remembered that the methodology described had the purpose of selecting the most adequate interpolation method and respective parameterization from a limited set of methods available in commercial GIS software packages. Two points must therefore be made explicit: (1) we did not intend to be exhaustive as to the types of interpolations compared. Other GIS software packages offer interpolation methods differing from the ones considered here. Our emphasis is on the proposed methodology, which may be applied to any method. (2) By considering only commercially available tools, we necessarily left aside a number of methods that may produce significantly better results. As an example, the HASM (High Resolution Spatial Modelling) technique developed by Yue et al. (2010b) and further enhanced by Yue et al. (2010a) produced smaller errors than the TIN and Spline methods in reproducing experimental topography, and was further tested as a void filling technique for SRTM data by Yue et al. (2012). Significantly, it produced better results in three contrasting topographic contexts than all the former methods as well as the ANUDEM technique, the algorithm behind the Topo to Raster tool in ArcGIS 9.1. herein selected as optimal for our study area.

5 Conclusions

In accordance with the objectives originally defined, the proposed methodology allowed to reconstruct the topography preceding the development of a set of 90 large gullies and gully complexes using an optimal interpolation method, and thus minimizing errors.
Apart from the validation-based approach that constitutes the core of the methodology, allowing any set of potential methods to be compared, the stages identified in Figure 3 provide a general outline that should be adapted to different situations. For example, although ArcGIS was used, numerous other GIS applications exist that could be applied, either alone or in combination (e.g. Surfer, GRASS GIS, ANUSPLIN for some examples of applications offering Spline interpolation). At the same time, the exact number of methods/parameterizations first adopted for comparison and afterwards refined in the parameter optimization stage is arbitrary, and will depend on available time, software and computational capacity.
The contrasting results obtained underline the importance of comparing both different methods and different parameterizations in each situation, as small changes in the user-defined input values can significantly influence error.
Results also evidenced the limitations of linear interpolation methods for reconstructing pre-erosion topography in the studied area. Although it is expectable that similar limitations will be extended to most geomorphic settings, a linear approach may eventually prove more adequate than Splines in contexts of flat topography, markedly rectilinear hillslopes or breaks in slope. The preliminary analysis as to the potential applicability of each method that was defined as the first stage of the methodology (Figure 3) should therefore not be disregarded in future applications.
The proposed methodology can easily be applied to different topographic contexts, enabling better calculations of eroded volumes and denudation rates. At the same time, the interpolated surfaces can be used to investigate controls by antecedent topographic properties (e.g. profile and planform curvature) over erosive processes.

The authors have declared that no competing interests exist.

References

Publishing order | Descend order by publishing year | Descend order by cited within
1
Alexander R W, Calvo-Cases A, Arnau-Rosalén Eet al., 2008. Erosion and stabilization sequences in relation to base level changes in the El Cautivo badlands, SE Spain.Geomorphology, 100: 83-90.

2
Anderson M G, Burt T P, 1978. The role of topography in controlling through flow generation.Earth Surface Processes, 3: 331-344.

3
Bergonse R, Reis E, 2011. Theoretical constraints to gully erosion research: Time for a re-evaluation of concepts and assumptions?Earth Surface Processes and Landforms, 36: 1554-1557.

4
Betts H D, Trustrum N A, De Rose R C, 2003. Geomorphic changes in a complex gully system measured from sequential digital elevation models, and implications for management.Earth Surface Processes and Landforms, 28: 1043-1058.

5
Bou Kheir R, Wilson J, Deng Y, 2007. Use of terrain variables for mapping gully erosion susceptibility in Lebanon.Earth Surface Processes and Landforms, 32: 1770-1782.

6
Buccolini M, Coco L, 2010. The role of the hillside in determining the morphometric characteristics of “calanchi”: The example of Adriatic central Italy.Geomorphology, 123: 200-210.

7
Buccolini M, Coco L, Cappadonia Cet al., 2012. Relationships between a new slope morphometric index and calanchi erosion in northern Sicily, Italy.Geomorphology, 149/150: 41-48.

8
Daba S, Rieger W, Strauss P, 2003. Assessment of gully erosion in eastern Ethiopia using photogrammetric techniques.Catena, 50: 273-291.

9
De Rose R C, Gomez B, Marden Met al., 1998. Gully erosion in Mangatu Forest, New Zealand, estimated from digital elevation models.Earth Surface Processes and Landforms, 23: 1045-1053.

10
Dunne T, Moore T R, Taylor C H, 1975. Recognition and prediction of runoff-producing zones in humid regions.Hydrological Sciences Bulletin, XX(3): 305-327.

11
Eastman J R, 2006. IDRISI Andres Guide to GIS and Image Processing. Clark University.

12
ESRI (Environmental Systems Research Institute), 2005

13
Hengl T, Evans I S, 2009. Mathematical and digital models of the land surface. In: Hengl T, Reuter H I. Geomorphometry: Concepts, Software, Applications. Elsevier, 31-63.

14
Heywood I, Cornelius S, Carver S, 2006. An Introduction to Geographic Information Systems. 3rd ed. Pearson Education.

15
Hutchinson M F, 1988. Calculation of hydrologically sound digital elevation models. Paper presented at the Third International Symposium on Spatial Data Handling. Sydney, Australia.

16
Hutchinson M F, 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits.Journal of Hydrology, 106: 211-232.

17
Kakembo V, Xanga W W, Rowntree K, 2009. Topographic thresholds in gully development on the hillslopes of communal areas in Ngqushwa Local Municipality, Eastern Cape, South Africa.Geomorphology, 110: 188-194.

18
Martínez-Casasnovas J A, Ramos M C, Garcia-Hernández D, 2009. Effects of land-use changes in vegetation cover and sidewall erosion in a gully head of the Penedès region (northeast Spain).Earth Surface Processes and Landforms, 34: 1927-1937.

19
Neteler M, Mitasova H, 2008. Open Source GIS: A GRASS GIS Approach. Springer.

20
Olaya V, 2009. Basic land-surface parameters. In: Hengl T, Reuter H (eds.). Geomorphometry: Concepts, Software, Applications. Elsevier, 141-169.

21
Parkner T, Page M, Marutani Tet al., 2006. Development and controlling factors of gullies and gully complexes. East Coast, New Zealand.Earth Surface Processes and Landforms, 31: 187-199.

22
Perroy R L, Bookhagen B, Asner G Pet al., 2010. Comparison of gully erosion estimates using airborne and ground based LiDAR on Santa Cruz Island, California.Geomorphology, 118: 288-300.

23
Poesen J, Nachtergaele J, Verstraeten Get al., 2003. Gully erosion and environmental change: Importance and research needs.Catena, 50: 91-133.

24
Reuter H I, Nelson A, 2009. Geomorphometry in ESRI packages. In: Hengl T, Reuter H I (eds.). Geomorphometry: Concepts, Software, Applications. Elsevier, 270-291.

25
Summerfield M, 1991. Global Geomorphology. Pearson Education.

26
Weibel R, Heller M, 1991. Digital terrain modelling. In: Maguire D J, Goodchild M F, Rhind D W (eds.). Geographical Information Systems: Principles and Applications: 269-297.

27
Wells N, Andriamihaja B, Rakotovololona H F, 1991. Patterns of development of lavaka, Madagascar’s unusual gullies.Earth Surface Processes and Landforms, 16: 189-206.

28
Wells S G, Gutierrez A A, 1982. Quaternary evolution of badlands in the south-eastern Colorado plateau, USA. In: Bryan R, Yair A (eds.). Badland Geomorphology and Piping. GeoBooks.

29
Wood A L, Simon A, Downs P Wet al., 2001. Bank-toe processes in incised channels: the role of apparent cohesion in the entrainment of failed bank materials.Hydrological Processes, 15: 39-61.

30
Yao C, Lei T, Elliot W Jet al., 2008. Critical conditions for rill initiation.Transactions of the ASABE, 51(1): 107-114.

31
Yue Tianxiang, Chen Chuanfa, Li Bailian, 2010a. An adaptive method of high accuracy surface modeling and its application to simulating elevation surfaces.Transactions in GIS, 14(5): 615-630.

32
Yue Tianxiang, Chen Chuanfa, Li Bailian, 2012. A high-accuracy method for filling voids and its verification.International Journal of Remote Sensing, 33(9): 2815-2830.

33
Yue Tianxiang, Song Dunjiang, Du Zhengpinget al., 2010b. High-accuracy surface modelling and its application to DEM generation.International Journal of Remote Sensing, 31(8): 2205-2226.

34
Zucca C, Canu A, Della Peruta R, 2006. Effects of land use and landscape on spatial distribution and morphological features of gullies in an agropastoral area in Sardinia (Italy).Catena, 68: 87-95.

Options
Outlines

/

〈 〉