Accuracy Assessment of Digital Terrain Model Dataset Sources for Hydrogeomorphological Modelling in Small Mediterranean Catchments

Abstract

Digital terrain models (DTMs) are a fundamental source of information in Earth sciences. DTM-based studies, however, can contain remarkable biases if limitations and inaccuracies in these models are disregarded. In this work, four freely available datasets, including Shuttle Radar Topography Mission C-Band Synthetic Aperture Radar (SRTM C-SAR V3 DEM), Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Map (ASTER GDEM V2), and two nationwide airborne light detection and ranging (LiDAR)-derived DTMs (at 5-m and 1-m spatial resolution, respectively) were analysed in three geomorphologically contrasting, small (3–5 km²) catchments located in Mediterranean landscapes under intensive human influence (Mallorca Island, Spain). Vertical accuracy as well as the influence of each dataset’s characteristics on hydrological and geomorphological modelling applicability were assessed by using ground-truth data, classic geometric and morphometric parameters, and a recently proposed index of sediment connectivity. Overall vertical accuracy—expressed as the root mean squared error (RMSE) and normalised median deviation (NMAD)—revealed the highest accuracy for the 1-m (RMSE = 1.55 m; NMAD = 0.44 m) and 5-m LiDAR DTMs (RMSE = 1.73 m; NMAD = 0.84 m). Vertical accuracy of the SRTM data was lower (RMSE = 6.98 m; NMAD = 5.27 m), but considerably higher than for the ASTER data (RMSE = 16.10 m; NMAD = 11.23 m). All datasets were affected by systematic distortions. Propagation of these errors and coarse horizontal resolution caused negative impacts on flow routing, stream network, and catchment delineation, and to a lower extent, on the distribution of slope values. These limitations should be carefully considered when applying DTMs for catchment hydrogeomorphological modelling.

1. Introduction

Digital terrain models (DTM) provide a continuous mathematical representation of the Earth’s bare surface [1]. In catchment hydrogeomorphological modelling, DTMs are the most important data inputs, as the accurate replication of landscape elements highly influences model accuracy [2,3,4,5]. Adequate representation of landscape elements is often limited by intrinsic errors related to the data acquisition process as well as preprocessing steps used to generate the DTM [6,7]. In fact, DTM dataset characteristics such as spatial resolution and vertical accuracy, as well as the accurate filtering of vegetation and human-made structures, may importantly affect terrain representation and subsequent modelling [8,9,10,11,12,13,14].

Most of the DTMs used in hydrogeomorphological applications are satellite-derived datasets, such as the SRTM DEM [15] and the ASTER GDEM [16], due to their nearly global topographic coverage and unrestricted availability. However, their spatial resolution is coarse (~30 m) and their targeted vertical accuracy is ~16 m [17] and ~17 m [18], respectively. Airborne light detection and ranging (LiDAR)-derived DTMs are recognised to overcome these restrictions in spatial resolution and vertical accuracy [19,20,21,22], but their spatial coverage is limited due to their cost-intensive acquisition and processing. Nationwide, LiDAR DTMs are becoming increasingly available for some regions (e.g., several countries in the European Union). However, the use of these datasets is frequently limited to small-scale landscape modelling, such as landslide and hillslope failure assessment [23,24,25]. A broader range of LiDAR DTM potential applications in catchment hydrology and landscape geomorphology has been identified [26,27].

At the catchment scale, studies focusing on vertical accuracy and DTM effects on hydrogeomorphological modelling have intensively applied satellite-derived SRTM and ASTER DTM data. These studies have been conducted in a wide range of geomorphologically contrasting study areas, including low-relief terrains [28,29,30] and steep-gradient areas [31,32,33], as well as different geographic regions, such as mid-latitude [34,35,36], tropical, and subtropical areas [37,38,39]. However, Mediterranean environments are clearly underrepresented in such validation and comparison studies, despite the few relatively old studies [40,41] and several more recent investigations about the detection of karst features in satellite-derived DTMs [42]. This underrepresentation must be seen as a critical point, as the Mediterranean region is arguably one of the most human-influenced areas on Earth. For instance, a major characteristic of this region is the massive presence of traditional soil and water conservation structures constructed since Roman times [43]. However, changes in traditional land use systems, such as urban expansion and increases in irrigated agriculture [44], have altered complex geoecosystems, thereby favouring the progress of land abandonment and further land degradation [45,46] as well as profound changes in hydrological systems [47].

In order to model current and future changes in the hydrogeomorphological behaviour of Mediterranean environments, the use of reliable and accurate DTM is fundamental. Accordingly, the objectives of this study were twofold. First, ground-truth assessment of the vertical accuracy of SRTM C-SAR DEM V3, ASTER GDEM V2 (spatial resolution 1 arcsecond ~ 30 m at the equator), and two airborne LiDAR DTM datasets with 5- and 1-m spatial resolution, respectively, was applied in the island of Mallorca (Spain) over three small catchments representative of characteristic Mediterranean landscapes with morphologies and land uses differently affected by a long-term history of human changes. Second, evaluation of the studied DTM data sources and dataset characteristics on hydrogeomorphological modelling applications was carried out over the analysed catchments by calculating widely used statistics and descriptors (e.g., hypsometry, slope–area relationships, stream network patterns, and distribution of the hillslope length and steepness factor values). Last, representative plots were selected within the catchments to investigate the role of characteristic landscape features, including human-made objects (e.g., agricultural terraces) and dominant geomorphological features (e.g., major rock outcrops), on the sediment connectivity index proposed by Borselli et al. [48] and improved by Cavalli et al. [49], providing critical assessment of small-detail patterns of water and sediment transference between different landscape compartments of the catchment. The results obtained by this study aim to expand the existing knowledge on the influence of DTM data sources and dataset characteristics on terrain representation for Mediterranean landscapes, while also providing valuable information about the application of remote-sensing-derived DTMs in small-catchment hydrogeomorphological modelling in and beyond Mediterranean environments.

2. Study Area

Three small contrasting catchments (Sa Font de la Vila, Es Fangar, and Es Telègraf; ranging from 3 to 5 km²) located in the Island of Mallorca (Spain) were selected for this study (Figure 1a). All three catchments are located in the alpine Tramuntana Range, although each of them is characterised by a unique underlying morphology with decreasing terrain complexity and land use intensity from the mountainous Es Telègraf and over the mixed Es Fangar to the terraced Sa Font de la Vila catchment. This gradient in terrain complexity and land use intensity is caused by the underlying geological settings of the Tramuntana Range, which is aligned in a NE to SW direction with several folds and thrusts with increasing relief energy and hillslope steepness to the northern part (Figure 1b; [50]).

The southernmost catchment is Sa Font de la Vila (4.8 km²). This catchment is located in the southwestern part of Mallorca (2°24′50″E, 39°35′20″N; Figure 1c). The relief of the catchment is complex due to the interaction of soft and hard lithology and different tectonics [50]. The catchment is characterised by the massive presence of traditional agricultural terraces, covering 37% of the whole catchment area and supported by 147 km of dry-stone walls [51]. Since the mid-twentieth century, the abandonment of traditional agriculture in these marginal areas caused an afforestation process with increased fuel availability [52]. As a result, two severe wildfires affected the catchment in 1994 and 2013. Current land use is dominated mainly by agriculture (44%), sparsely vegetated terrains (38%), and forests (18%).

The catchment with the highest relief energy and steepest hillslopes is the mountainous Es Telègraf catchment (2°51′0″E, 39°48′30″N; Figure 1c), which is located in the highest section of the Tramuntana Range (Figure 1b). The main geomorphological traits of the mountainous Es Telègraf catchment (ca. 3 km²) are aligned along a tectonic structure dominated by a NW-directed thrust system, including several cliffs and steep slopes. Vegetation shows a sharp contrast between the headwaters (25% of total area), dominated by sparse vegetation (culminal Balearic stage; [53]), and the lower parts (75% of total area), mostly covered by forest.

The mixed Es Fangar catchment (3.2 km²; 2°60′00″E, 39°50′10″N; Figure 1b) is characterised by a combination of thrust and normal faults and synclinal–anticlinal structures, as well as by the presence of different land use classes including forests (31%), sparsely vegetated areas (20%), and agricultural fields (49%), which are also affected by check-dam terraces and surface-leveled fields. Thus, the catchment shows a mosaic of different land uses and management practices as a result of a complex interaction between human and natural conditions [46].

Three small 5-ha plots, including characteristic landscape elements of the three catchments (Figure 1c), were selected to assess how accurate DTM representation of typical Mediterranean landscape elements—induced either by natural geomorphological conditions or human activities—can affect small-detail hydrogeomorphological modelling by using a morphometric hydrosedimentological connectivity index (see Section 3.3). The selected plot for the Sa Font de la Vila catchment is located in a terraced hillslope (Figure 1d). For the mountainous Es Telègraf catchment, the plot is selected in an area of high mountainous relief with bedrock outcrops and steep hillslopes of ca. ≥60% gradient slope (Figure 1e). This plot is centered on a massive bedrock outcrop formed by erosion-resistant limestone, while softer materials dominate the surrounding area. In the mixed Es Fangar catchment, the plot encompasses a flat agricultural field covered by rainfed herbaceous crops that are affected by traditional drainage systems, combining human-made artificial channels and subsurface tile drains [54] (Figure 1f).

3. Materials and Methods

3.1. DTM Datasets

The present study compares the widely used satellite-derived SRTM DEM V3 and ASTER GDEM V2 datasets as well as the nationwide, Spanish airborne LiDAR DTM datasets over the study catchments in the Island of Mallorca. SRTM data was acquired during February 2000 during an 11-day campaign of the National Aeronautics and Space Administration’s (NASA) space shuttle Endeavour, carrying both C-band (5.6 cm) and X-band (3.1 cm) synthetic aperture RADAR (SAR) devices as payload [15]. A first release of the SRTM C-SAR data was made available fpr public use by the government of the United States in 2003, providing nearly global topographic coverage with an initial horizontal resolution of 3 arcseconds (~90 m at the equator). Since 2014, the SRTM data has been available with the increased spatial resolution of 1 arcsecond (~27 m at Mallorca). The SRTM mission requirements have a targeted vertical accuracy of 16 m [55,56], but validation studies reported the actual vertical accuracy ranging between 3 and 12 m [29,30,31,32,33,35,39,57,58]. For use in this study, 1-arcsecond SRTM C-SAR version 3 data was downloaded from the USGS EROS, Earth Resources and Science Center at the United States Geological Survey.

Photogrammetry was employed in a collaborative effort of the NASA Jet Propulsion Laboratory (NASA-JPL) and the Japanese Ministry of Trade and Industry (METI) to produce the ASTER GDEM product from stereographic pairs of images in 1-arcsecond spatial resolution. A first release of the ASTER GDEM data was published in 2009 with an absolute vertical accuracy of 18.3 m. The second release of the dataset (ASTER GDEM V2) was published in 2011, with an improved vertical accuracy of 17 m [16]. Higher vertical accuracy was achieved by the inclusion of an additional number of 260,000 ASTER scenes and enhanced processing algorithms [18]. The reported vertical accuracy of the second release of ASTER GDEM was confirmed in the literature [38,39,58]. ASTER GDEM V2 was obtained for this study from the NASA tool. The number of stereographic pairs of ASTER images involved in DTM generation (stack number) for our study areas was between 11 and 27 in Sa Font de la Vila, between 2 and 7 in Es Telègraf, and between 5 and 9 in Es Fangar.

Airborne LiDAR data was acquired during an airborne campaign carried out by the Spanish National Institute of Geography (IGN) in 2014, covering the whole island of Mallorca [59]. The LiDAR dataset for the study areas is provided free of charge (i) as an already-processed DTM model with a spatial resolution of 5 m and as (ii) a preclassified point cloud by the Balearic Islands Autonomic Government and IGN. These nationwide IGN LiDAR datasets are just available for Spain, but are representative of most European airborne topographic LiDAR datasets, as evidenced by the European Union’s contribution to GMES (Global Monitoring for Environment and Security) services [60]. Minimum density of the point cloud is reported to be 0.5 points per m² [59]. In this study, an average density of 1.65 points per m² was calculated from the available LiDAR point cloud data of the catchments. Both 5-m DTM and point-cloud data was obtained from the IGN Geoinformatics webpage. Vertical accuracy of the 5-m DTM data is reported to be better than 0.2 m [59]. In order to generate a second LiDAR DTM model with 1-m spatial resolution, the preclassified point data was interpolated to a regularly spaced grid using multilevel B grid spline interpolation [61] with a threshold level of 0.0001 and a maximum level of 11 iterations. The resulting 1-m terrain model was smoothed by a multidirectional Lee filtering algorithm [62] to reduce the amount of noise.

Accordingly, four different datasets were analysed in this study: 1-arcsecond SRTM C-SAR DEM V3; 1-arcsecond ASTER GDEM V2; 5-m LIDAR IGN; and 1-m LIDAR IGN (hereafter referred to as SRTM, ASTER, IGN 5 m, and IGN 1 m, respectively). The datasets differ notably in their data source, spatial resolution, and vertical accuracy, so they provide an excellent study framework for investigating their applicability for catchment hydrogeomorphological modelling.

3.2. Vertical Accuracy Assesment

Vertical accuracy was assessed by collecting ground control points (GCPs) using a Real-Time Kinematic Global Positioning System (RTK-GPS, Leica GPS Series 1200) with a base station placed at central locations of the study catchments and absolute accuracy in the range of centimetres. Before surveying in September 2017, the initially established random spatial distribution of the GCPs was checked against the CORINE land use classes and 2012 cover data available from the COPERNICUS Land Monitoring Service and high-resolution (0.25 m/pixel) aerial imagery facilitated by the Government of the Balearic Islands. GCPs were adjusted manually later to ensure a representative proportion of densely vegetated (i.e., forest cover) and open field areas (i.e., sparse vegetation and agricultural land use cover classes) in the survey. At Sa Font de la Vila, 32 GCPs that were previously surveyed in January 2016 over terraced hillslopes using the same equipment and procedures were also included in this study. As a result, a total of 140 GCPs were included in our DTM accuracy assessments (Figure 1c).

After field acquisition, GPS measurements were rectified from ellipsoidal to orthometric height to ensure comparability with the DTM datasets that refer to the EGM1996 geoid (ASTER and SRTM) and the EGM2008-REDNAP geodetic vertical data (IGN data). The robustness of the accuracy assessment in terms of statistic measurements was ensured by applying a threshold of a maximum absolute vertical error of 7 cm in the GPS data, which is nearly 3 times more accurate than the accuracy given by IGN for the 5-m LiDAR data [59]. Under this condition, the derived DTM accuracy had an error of 5%, which is considered as tolerable [63].

Vertical accuracy of the DTM datasets was expressed as error statistics. In order to determine if the underlying distribution of elevation errors equals a normal distribution, histogram plots were examined visually. Robust measures such as the normalised median deviation (NMAD) and the 68.3% and 95.0% sample quantiles of the error between GPS- and DTM-derived elevation were computed. All these measures were reported by Höhle and Höhle [63] and Höhle [64] to be reliable even for non-normal distributions. NMAD is computed according to:

$$NMAD=1.4826\ast median(|∆-\widetilde{∆}|)$$

(1)

where $\widetilde{∆}$ denotes the median in elevation error $∆$, computed from the $i\le n$ differences between GCP- and DTM-derived elevation values. For comparability with other accuracy assessments, the root mean squared error (RMSE) is also provided:

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left({∆}_i\right)2}$$

(2)

RMSE and NMAD values were provided at the levels of the DTM source (SRTM, ASTER, IGN 5 m, and 1 m), catchment (Sa Font de la Vila, Es Telègraf, and Es Fangar), and land use type (densely vegetated forest cover versus open-field agricultural plus sparse vegetation cover classes) for analysis.

3.3. Quality Assessment of DTMs for Hydrogeomorphological Modelling

All the different DTM datasets were projected to UTM Zone 31 N (ETRS89 ellipsoid) to ensure comparable geolocation. The original vertical data of the datasets remained unchanged as the difference in height between EGM1996 and EGM2008-REDNAP were <1 m in the study areas, which can be considered small if compared to the vertical resolution and precision of the satellite-derived DTMs. In order to correct the data for hydrogeomorphological purposes, the Planchon and Darboux [65] surface filtering algorithm was applied, imposing a minimum slope gradient between grid cells of 0.1%. The algorithm adds a virtual layer of water to the DTM data, filling all depressions up to the level of their outflow point. This results in a surface that ensures drainage even for flat areas. Comparison of elevation statistics before and after pit-filling was carried out to assess the influence of surface filtering.

A broad variety of hydrogeomorphological statistics and descriptors were analysed to explore the feasibility of the studied DTMs for hydrogeomorphological modelling in the three study catchments. These indicators encompassed basic terrain characteristics, widely used catchment geomorphometric parameters and relationships, and stream network and flowpath properties, as well as small-detail arrangement patterns of water and sediment fluxes (

Table 1. Summary of statistics and descriptors applied to analyse the feasibility of the studied DTMs for catchment hydrogeomorphological modelling, and relevant references for the methods.

Hydrogeomorphological Statistics and Descriptors	References
Basic terrain characteristics
Minimum, maximum, mean elevation, and total relief
Mean and SD of slope gradient
Mean and SD (D8*) of flowlength	O’Callaghan and Marks [66]
Total catchment area	O’Callaghan and Marks [66]
Catchment geomorphometric parameters and relationships
Terrain hypsometry	Strahler [67]
Slope–area relationship	Hancock et al. [68], Willgoose [69]
Cumulative area distribution	Rodríguez-Iturbe et al. [70]
Mean and SD of LS** factor	Wischmeier and Smith [76], Desmet and Govers [77]
Stream network and flowpath properties
Stream network patterns	O’Callaghan and Marks [66]
Cumulative distribution function of flowpath lengths	Moreno-de-las-Heras et al. [28]
Small-detail water/sediment flow arrangement patterns
Surface connectivity index (IC)	Borselli et al. [48], Cavalli et al. [49]
D8 = D8 flow direction algorithm [66]; ** LS = hillslope length and steepness factor [76]

). Extraction of basic terrain characteristics, such as general elevation statistics, slope gradient and flow direction, upslope contributing area, and catchment area were computed using the ESRI ArcGIS (Version 10.3) spatial analyst toolbox. General flow accumulation for analysis was obtained by applying the deterministic nondispersive D8 algorithm [66] that simulates water flow in the direction of the strongest slope gradient.

The hypsometric curve [67], slope–area relationship [68,69], cumulative area distribution function [70], and empirical hillslope length and slope steepness factor (LS factor [76,77]) were computed, providing geomorphometric evaluation criteria for analysis. The hypsometric curve relates relative catchment area to relative elevation and allows the estimation of dominant erosion processes [71] and landform maturity [72]. Thus, it is a simplistic measure of the mass and energy stored within a landscape. The slope–area relationship—defined by catchment mapping of the local slope gradient to the contributing area—and the cumulative area distribution—defined as the proportion of the catchment that has a contributing area greater than or equal to a specified drainage area—provide information about characteristic fluvial processes [69] and flow aggregation structures [70], respectively. In more detail, the analysis of slope–area plots generally reveals the emergence of an inflection point that separates relatively small catchment areas dominated by the action of diffusive erosion from larger catchment areas characterised by fluvial erosion and transport [69]. In theory, the shape of the slope–area plots of catchments is a “boomerang”, with the roll-over (inflection) point corresponding to the threshold in drainage area where hillslopes transition into channels [73]. Differently, the cumulative area relationship focuses on flow aggregation structures of hillslope elements and the stream network [74,75]. Herein, the leftmost section of the distribution curve is associated with flow aggregation at hillslopes, whereas the middle section describes channelised flow as a log–log linear straight line, and the rightmost section is related with boundary effects. Such an organisation of flow structures can be determined for all DTM datasets and catchments [74].

The empirical hillslope length and slope steepness factor (LS factor) provides critical information on the terrain morphological susceptibility to surface soil erosion. In its original formulation, as part of the universal soil loss equation (USLE) [76], the LS factor is determined as:

$$LS={(\frac{l}{a_0})}^p\ast \left(0.065+4.56\ast \sin \left(\theta \right)+65.41\ast \sin 2\left(\theta \right)\right)$$

(3)

where $l$ accounts for the hillslope length (usually calculated from flowlength), $a_0$ is the standardised plot length applied in the USLE experiment (~22.1 m), Ɵ is the slope gradient (rad), and $p$ is the slope length exponent, which ranges between 0 and 1 and approaches 0 when the ratio of rill-to-interrill erosion is close to 0. Desmet and Govers [77] extended the classical LS approach to large-scale modelling applications in a Geographic Information System (GIS) by replacing hillslope flowlength by upslope contributing area. We applied the improved Desmet and Govers’ LS approach, assuming a rill-to-interrill erosivity ratio of 1 for LS calculation in SAGA GIS [78].

Differently from the above-described geomorphometric properties that describe the dominant erosion processes, landform maturity, flow aggregation structures, and terrain morphological susceptibility to surface erosion, stream network properties provide important clues for understanding the large-scale spatial arrangement of water and sediment fluxes. Flow routing differences between the different DTM models were analysed for each study catchment by exploring variations in stream network representation and flowpath length distribution. Stream networks were delineated from the derived D8 flow-accumulation raster maps of the catchments by imposing a minimum size of upslope contributing area of 30 ha, which produced the most realistic results when compared to field observations. In addition, the empirical inverse cumulative distribution function (iCDF) of derived flowlength values for each studied catchment and DTM were determined for analysis, following the approach of Moreno-de-las-Heras et al. [28].

Last, the connectivity index (IC) proposed by Borselli et al. [48] and modified by Cavalli et al. [49] was applied for exploring small-detail patterns of water and sediment transference between different landscape compartments of the study catchments, considering that a dynamic assessment of landscape connectivity can help to incorporate aspects of the process linkages that essentially drive water and sediment fluxes [79]. IC is a dynamic property of the catchment that indicates the probability of a particle at a certain location to reach a defined target area, which in this study, was established at the catchment outlet, with its effects being analysed at representative small-scale (5-ha) plots from each study catchment (Figure 1d–f). IC is calculated as:

$$IC={\log }_{10}\left(\frac{D_{up}}{D_{dn}}\right)={\log }_{10}\left(\frac{\overline{W}\overline{S}\sqrt{A}}{{\sum }_i\frac{di}{WiSi}}\right)$$

(4)

where $D_{up}$ and $D_{dn}$ are the up- and downslope components of the connectivity index, respectively, and $\overline{S}$ represents the average percentage slope, A the size of the upslope contributing area, $\overline{W}$ an averaged weighting factor represented by the terrain roughness, and d_i the length of the flowpath from the ith cell to the targeted outlet along the steepest downslope direction. IC was determined for this study using the freely available SedInConnect (Version 2.3) software developed by Crema and Cavalli [80].

4. Results

4.1. DTM Vertical Accuracy

Significant errors in elevation for all the analysed datasets were observed by comparing DTM-derived elevation values against those elevation values obtained in our dGPS-derived GCPs (

Table 2. Results of the vertical accuracy assessment expressed as root mean squared error (RMSE) and normalised median deviation (NMAD) between the GPS-measured elevation values and those values derived from the different DTM datasets (in metres).

	SRTM DEM	ASTER GDEM	IGN 5 m	IGN 1 m
All sites (n = 140)
RMSE	6.98	16.10	1.73	1.55
NMAD	5.27	11.23	0.84	0.44
All open terrain sites (agricultural + sparse vegetation cover, n = 87)
RMSE	7.38	16.26	1.59	1.41
NMAD	5.22	10.82	0.93	0.61
All densely vegetated sites (forest cover, n = 53)
RMSE	6.28	15.84	1.94	1.76
NMAD	5.31	11.46	0.73	0.33
Sa Font de la Vila catchment (n = 53)
RMSE	8.28	9.62	2.09	2.03
NMAD	5.59	7.63	0.98	0.62
Es Telègraf catchment (n = 40)
RMSE	7.76	26.77	1.59	1.25
NMAD	5.76	12.17	0.89	0.49
Es Fangar catchment (n = 47)
RMSE	4.06	7.62	1.35	1.10
NMAD	2.98	7.72	0.71	0.26

). Accuracy increased from ASTER over SRTM to the 5-m and 1-m IGN LiDAR models. A systematic overestimation of elevation values was observed for all datasets, as indicated by the relative frequency distribution of (dGPS minus DTM) vertical errors (i.e., the error distributions are skewed to the left for all DTMs as shown in Figure 2).

The SRTM data showed an overall accuracy of 6.98 m RMSE and 5.27 m NMAD (Table 2), also evidencing a strong tendency to overestimate surface elevation (up to 15–25 m; Figure 2). Among the three studied catchments, SRTM reached the lowest accuracy in terms of RMSE and NMAD (over 7.70 and 5.50 m, respectively) for the terraced Sa Font de la Vila catchment and the mountainous Es Telègraf catchment. The highest accuracy of the SRTM data was observed in the mixed Es Fangar catchment (RMSE = 4.06 m, NMAD = 2.98 m; Table 2). No clear influence of vegetation was observed in the SRTM errors, since the elevation bias showed the same relative distribution patterns (Figure 2) as well as similar RMSE and NMAD values (Table 2) for forest and open (agricultural and sparse vegetation cover) areas.

ASTER showed the lowest overall accuracy among all datasets, according to the obtained RMSE and NMAD values (16.10 and 11.23 m, respectively; Table 2). The frequency distribution of the elevation errors showed that ASTER values largely overestimated elevation (by up to 30 m), although significant underestimations were also identified for the dataset (Figure 2). The lowest accuracy was observed at the mountainous Es Telègraf catchment (RMSE = 26.77 m, NMAD = 12.17 m; Table 2). The mixed Es Fangar catchment showed the best accuracy in terms of RMSE (7.62 m), whilst the terraced Sa Font de la Vila catchment displayed the best NMAD performance (7.63 m). Comparing open and densely vegetated terrains, dense vegetation considerably increased the overestimation of elevation values (i.e., the histogram green bars are more skewed to the left than red bars in Figure 2), although both RMSE and NMAD values were similar for forest and open-terrain points (Table 2).

The vertical accuracy of the LiDAR-based IGN DTMs was considerably higher than for the SRTM and ASTER data. The IGN 5 m model showed a vertical accuracy of 1.73 and 0.84 m, expressed as RMSE and NMAD, respectively (Table 2). Differently to the SRTM and ASTER data, the vertical accuracy for the IGN 5 m increased from the terraced Sa Font de la Vila catchment (RMSE = 2.09 m, NMAD = 0.98 m) to the mountainous Es Telègraf (RMSE = 1.59, NMAD = 0.89) and the mixed Es Fangar (RMSE = 1.35, NMAD = 0.71 m) catchments. The IGN 1 m model, which was developed from the same LiDAR data cloud as the 5-m DTM, showed the best model performance in terms of vertical accuracy (RMSE = 1.55 m, NMAD = 0.44 m), performing slightly more accurately than the 5-m data (Table 2). Similarly to the IGN 5-m model, the results for the 1-m model showed the lowest accuracy in the terraced Sa Font de la Vila catchment (RMSE = 2.03 m, NMAD = 0.62 m), while reaching the best model performance in the mixed Es Fangar catchment (RMSE = 1.10 m, NMAD = 0.26 m). Both the IGN 5-m and 1-m models showed a slight tendency to overestimate elevation values, with densely vegetated areas showing slightly higher RMSE than open terrain (Table 2), but similar relative frequency distribution patterns of elevation errors for both (forest and open) land use classes (Figure 2).

4.2. DTM Hydrogeomorphic Modelling Performance

4.2.1. Basic Terrain Statistics

Descriptive statistics of minimum, maximum, and mean elevation values as well as average slope gradient values before and after applying surface filtering (i.e., pit-filling) are summarised in

Table 3. Descriptive statistics of the used DTM datasets. “Uncorr” refers to results before surface filtering and “Corr” after surface filtering “-” indicates not applicable measures.

		SRTM DEM		ASTER GDEM		IGN 5 m		IGN 1 m
		Uncorr	Corr	Uncorr	Corr	Uncorr	Corr	Uncorr	Corr
Sa Font de la Vila	Minimum Elevation (m)	71.5	71.5	71.0	71.0	66.6	66.6	66.0	66.3
	Maximum Elevation (m)	470.0	470.0	505.0	505.0	516.1	516.1	516.4	516.4
	Mean Elevation (m)	252.7	252.7	257.8	257.8	256.7	256.7	256.4	256.7
	Mean Slope (%)	30.9	30.8	33.4	33.1	39.7	39.7	40.5	40.5
	Relief (m)	398.5	398.5	434.0	434.0	449.5	449.5	449.8	450.1
	Catchment Area (km2)	-	5.06	-	4.83	-	4.82	-	4.83
	Filtered Area (%)	-	0.2	-	0.6	-	<0.1	-	<0.1
Es Telègraf	Minimum Elevation (m)	638.5	638.5	632.0	639.0	625.0	625.0	624.6	624.6
	Maximum Elevation (m)	1337.5	1337.5	1350.0	1350.0	1349.5	1349.5	1351.0	1351.0
	Mean Elevation (m)	947.4	947.4	921.4	921.4	912.7	912.7	911.6	911.6
	Mean Slope (%)	44.7	44.7	45.4	45.2	54.2	54.2	55.6	55.6
	Relief (m)	699.0	699.0	718.0	711.0	724.5	724.5	726.4	726.4
	Catchment Area (km2)	-	3.55	-	3.17	-	2.73	-	2.72
	Filtered Area (%)	-	0.1	-	0.7	-	<0.1	-	<0.1
Es Fangar	Minimum Elevation (m)	74.3	74.9	70.0	74.0	73.0	73.4	73.2	73.9
	Maximum Elevation (m)	369.9	369.9	405.0	405.0	403.8	403.8	403.9	403.9
	Mean Elevation (m)	157.8	157.8	163.6	163.7	163.0	163.0	162.7	162.7
	Mean Slope (%)	18.8	18.8	22.3	21.9	23.8	23.8	24.3	24.3
	Relief (m)	295.6	295.0	335.0	331.0	330.8	330.4	330.7	330.0
	Catchment Area (km2)	-	3.04	-	3.26	-	3.42	-	3.39
	Filtered Area (%)	-	0.3	-	1.7	-	<0.1	-	<0.1

. Surface filtering did not significantly affect the terrain models, as the percentage of pit-filled areas was between 0.1 and 0.3% for SRTM, 0.6 and 1.7% for ASTER, and <0.1% for the IGN LiDAR models. Since the extension of the DTM-derived catchment areas was partly different (Table 3), the absolute values of relief (difference between maximum and minimum elevation) cannot be directly compared. For the mountainous Es Telègraf catchment, SRTM and ASTER data showed larger catchment areas (i.e., 3.55 and 3.17 km², respectively) than the LiDAR models (i.e., 2.73 and 2.72 km² for the 5-m and 1-m models, respectively). Conversely, the ASTER-derived catchment area for the terraced Sa Font de la Vila and the mixed Es Fangar catchment coincided well with the LiDAR-derived values. Among all datasets, SRTM showed the largest catchment extent for Sa Font de la Vila (5.06 km²) and the smallest catchment area for Es Fangar (3.04 km²). The differences between the two LiDAR models were negligible in all the studied catchments. In addition to the observed discrepancies regarding catchment area, a general increase in average slope gradient was detected from the SRTM model and over the ASTER data to the 5-m and 1-m LiDAR models (Table 3).

4.2.2. Geomorphometric Parameters and Relationships

The values of both the hypsometric integral and LS factor decreased from the mountainous Es Telègraf catchment and over the terraced Sa Font de la Vila to the mixed Es Fangar catchment (

Table 4. Geomorphological and hydrological characteristics of the studied datasets and catchments, including the hypsometric integral, mean flowlength, and LS factor and their corresponding standard deviation (SD). “-“ mean that the measure has no physical unit.

		SRTM DEM	ASTER GDEM	IGN 5 m	IGN 1 m
Sa Font de la Vila	Hypsometric Integral (-)	0.46	0.44	0.43	0.43
	Mean Flowlength (m)	2304.3	2344.4	2520.2	2566.9
	SD Flowlength (m)	1020.1	969.3	1011.5	1031.0
	Mean LS factor (-)	5.4	4.8	6.7	12.1
	SD LS factor (-)	2.8	3.2	11.1	24.0
Es Telègraf	Hypsometric Integral (-)	0.45	0.41	0.41	0.41
	Mean Flowlength (m)	1760.3	1820.6	1821.3	1877
	SD Flowlength (m)	814.4	902.7	950.2	991.8
	Mean LS factor (-)	7.4	5.3	9.7	18.1
	SD LS factor (-)	9.5	4.5	33.1	37.3
Es Fangar	Hypsometric Integral (-)	0.27	0.28	0.28	0.28
	Mean Flowlength (m)	1554.3	1789.3	1668.3	1731.0
	SD Flowlength (m)	711.0	867.3	791.8	826.1
	Mean LS factor (-)	3.3	3.6	4.5	9.1
	SD LS factor (-)	2.7	2.8	7.9	20.6

). Differently, mean flowlength was largest at the Sa Font de la Vila catchment (Table 4), which also showed the largest catchment areas (Table 3). The hypsometric integral was nearly equivalent for all datasets, except for the mountainous Es Telègraf catchment, where SRTM showed a slightly higher integral value than the other models (Table 4). Differently, clear differences among the DTM datasets were detected for the (mean and SD) flowlength values of the catchments (Table 4). SRTM showed the smallest flowlength values in all catchments, followed by the ASTER data in Es Telègraf and Sa Font de la Vila, but not in the Es Fangar catchment. The IGN LiDAR models obtained the largest flowlength values in all catchments, except for the mixed Es Fangar catchment, where the ASTER model generated the largest value. A slight flowlength difference in the order of decametres was observed between the two IGN LiDAR models, with the 1-m dataset showing the highest (mean and SD) values (Table 4). Considering the LS factor, the LiDAR-based IGN data generated larger (mean and SD) values than the SRTM and ASTER datasets (Table 4). The LS values obtained by the IGN 1-m dataset in the studied catchments were (twofold) larger than those reached by the 5-m dataset (Table 4).

In agreement with the general patterns obtained using the values of the hypsometric integral, the hypsometric curves of the four elevation datasets showed a similar shape and trend for each of the three analysed catchments (Figure 3a). Differently, the slope–area and cumulative area distribution plots of the catchments showed important differences between DTMs. In fact, the slope–area plot showed the classic boomerang-like shape trend for the two (LiDAR-based) IGN models, separating the diffusive and fluvial-dominated section of the catchments (i.e., slope gradient slightly increased or kept constant values with contributing area in the 1 to 10³ m² section and then decreased linearly in the log–log plot), but only showed the fluvial-dominated (negative log–log trend) section for the SRTM and ASTER datasets (Figure 3b). Similarly, the cumulative area distribution plots for the IGN datasets accurately separated the hillslope and stream-dominated areas in the study catchments (depicted by the change in the slope of the cumulative distribution at 10²–10³ m² contributing area; Figure 3c). However, for the ASTER and SRTM datasets, a very small number of cells (leftmost section of the cumulative area plots) described the hillslope region of the catchments (Figure 3c).

4.2.3. Stream Networks

At the mountainous Es Telègraf catchment, SRTM produced straight-line stream features that did not reflect the actual terrain topography (Figure 4a). Furthermore, the satellite-derived SRTM and ASTER models indicated the existence of second-order streams that were not represented by the IGN LiDAR data for this catchment. It should be noted that Es Telègraf was also the catchment with the largest differences in delimited catchment area between the datasets (Table 2), where the SRTM and ASTER models largely failed to determine the actual boundary of the catchment, particularly in the southernmost contributing area region. In contrast, a better consistency between the models was obtained for the terraced Sa Font de la Vila catchment, despite some important shifts in the position of second-order streams at the headwaters and at the outlet point being detected (Figure 4a). The small area of interest represented in Figure 4a illustrates for this catchment how the streams derived from the IGN LiDAR data (green and yellow lines for the 5-m and 1-m IGN models, respectively) tended to reproduce the natural flow patterns of the landscape, following the valley floor as indicated by the 10-m contour lines, while the SRTM and ASTER streams (red and blue lines, respectively) partly deviated from those pathways. In the mixed Es Fangar catchment, the position of the outlet was identical for all DTMs, but again, remarkable differences were observed in the second-order streams, with the ASTER model showing the largest deviations (Figure 4a). Complementarily, the empirical inverse cumulative distribution function of flowlength probability (flowpath length iCDF; Figure 4b) revealed important differences among the datasets for Es Telègraf (i.e., shorter extreme flowpath lengths for the SRTM data) and Es Fangar (i.e., longer extreme flowpath lengths for the ASTER data).

4.2.4. Small-Detail Water and Sediment Connectivity

The coarser spatial resolution of the satellite-derived SRTM and ASTER models (~27 m) strongly limited the applicability of the analysed index of water/sediment connectivity (IC) for revealing small-detail surface flow patterns in the studied landscapes (Figure 5). Differently, the IC showed reliable surface flow patterns for the study sites using the LiDAR-based IGN DTMs, although the 1-m data considerably improved the details provided by the 5-m data (Figure 5). In accordance with these general results, a variety of differences characterised the IC performance of the SRTM, ASTER, and IGN LiDAR data for the small-detail 5-ha study plots. For example, for the mountainous Es Telègraf catchment, both the SRTM and LiDAR datasets showed a highly connected NW–SE flowpath that was enforced by the presence of a massive outcrop in the study site, while this dominant surface flow pathway was absent in the ASTER-derived IC map (Figure 5). The IC representation for this plot of the SRTM model was, however, limited as the high connectivity area appeared blurred and not as sharp as in the LiDAR-based IGN data. In Sa Font de la Vila, where the 5-ha small-detail plot was centred on a terraced hillslope, the SRTM and ASTER models failed to represent the effects caused by the presence of bench terraces on the surface flow patterns (Figure 5). The dendritic-like, highly connected network represented along the N–S direction by the IGN-derived maps of IC for this terraced site was replaced by a nearly uniform IC area in the SRTM and ASTER data. Furthermore, the effects of terracing on the IC patterns of this plot were better established by the 1-m data than the coarser 5-m LiDAR-based IGN model. Last, the IC results obtained in the 5-ha plot of Es Fangar again highlighted the severe difficulties of the ASTER model to reveal actual water and sediment pathways that are controlled by small-detail landscape objects. In this site, the small-detail plot was located in an agricultural field delimited by a road and a parallel artificial channel running along the SE corner of the landscape (Figure 5). Both the SRTM- and LiDAR-based IGN data showed a highly connected pathway running from NE to SW at the eastern border of the plot, reproducing the natural direction of the flow. Contrarily, the ASTER-derived IC values were highest along the northern part of the plot and did not show any sort of linear features associated to the actual water and sediment pathways of the landscape. Further comparisons of the SRTM- and LiDAR-derived IC patterns for this plot revealed that SRTM could not resolve minor features of the agricultural field (i.e., preferential flow patterns), indicating high connectivity crossing the plot from N to S. It should be also noted that even for the high-resolution 1-m LiDAR dataset, the position of the highest IC values did not exactly reproduce the actual position of the artificial channel in the landscape.

5. Discussion

5.1. Vertical Accuracy

5.1.1. General Dataset Evaluation

The vertical accuracy results of the studied (SRTM and ASTER) satellite- and (LiDAR-based IGN) airborne-derived models (Table 2) were in good agreement with other similar studies in hilly and mountainous areas [31,32,81], despite the fact that they were mainly focused in regions characterised by very different climate, vegetation patterns, and land use systems in Mediterranean areas. SRTM was less affected by systematic errors than ASTER, as indicated by the corresponding histogram plots, and the errors tended to show minor deviations from the normal error distribution (Figure 2). The SRTM data clearly outperformed the ASTER datasets, showing lower RMSE and NMAD values for all catchments and land use cover classes (Table 2), which, overall, agreed with previous analyses of these global datasets [39,58].

The vertical accuracy of the (5-m and 1-m) IGN LiDAR models (about 1.6 m RMSE and <1 m NMAD) was considerably higher than that obtained from the SRTM and ASTER datasets (accounting for 7–16 m and 5–11 m RMSE and NMAD, respectively). However, our analysis indicated that the actual vertical accuracy of the IGN DTMs is considerably lower than the (0.2 m) reference vertical accuracy that is officially posted by the IGN [59] for these datasets. The error magnitude (RMSE and NMAD) of the studied LiDAR models was in agreement with Estornell et al. [82] and Simpson et al. [22], showing that the 1-m dataset had the lowest RMSE and NMAD for all study sites and land use/cover classes. The LiDAR data showed significant shifts of the distribution in elevation errors towards negative values (Figure 2), although the slight overestimation in elevation caused by these errors was notably lower than the large overestimations produced by the satellite-derived SRTM and ASTER models. The increase in spatial resolution from 5 to 1 m in the LiDAR-derived models did not involve remarkable changes in the frequency distribution and magnitude of the errors. This may be related to the fact that both LiDAR DTMs were generated from the same 3D point-cloud and predefined classification scheme to separate ground from nonground returns; thus, identical classification errors are likely to affect both models [83,84].

5.1.2. Assessing the Effects of Vegetation on the Datasets

Classification of elevation errors according to land use/cover type revealed further information about the structure of the DTM datasets. The presence of densely vegetated forest cover exacerbated the tendency of the ASTER dataset to overestimate elevation values, reaching deviations of up to 30 m (Figure 2). In this context, Li et al. [90] defined the ASTER data as a “first-return system”. In fact, ASTER, as an optical remote-sensing system, captures the top of the vegetation canopy. Consequently, a large offset to the true surface elevation is expected to be obtained for the ASTER data in densely covered areas. Interestingly, in spite of exacerbating the tendency to overestimate elevation values in the ASTER data, the absolute vertical accuracy of this dataset was very similar for both densely vegetated (forest) and open (agricultural fields and sparse vegetation areas) surfaces. The poor elevation accuracy of the ASTER data over the studied catchments (16 and 11 m RMSE and NMAD values, respectively; Table 3) could be more influenced by the low amount of stereographic pairs (or stack number) applied to generate the elevation data in these sites, which ranged between 2 and 27 pairs. The use of low stack numbers—particularly where applying fewer than 10 stereographic pairs—in the production of ASTER GDEM data is reported to very remarkably reduce the vertical accuracy of the elevation values, as the amount of residual artefacts in the data can be very large [18,28,93]. In fact, the inherent, large vertical errors that affect the analysed ASTER data over the study catchments very likely surpass the influence of vegetation in measured absolute vertical accuracy. For example, although some relatively tall (up to 20 m) pine forest patches can be found in the studied areas, forest cover for the analysed catchments is more widely dominated by relatively open xeric formations of holm oak trees of up to 6–8 m of tree height, which are largely below the obtained ~16 m RMSE for the analysed ASTER data.

No clear effects of vegetation cover were detected in the absolute vertical accuracy and the relative distribution of the vertical errors for the SRTM data. Kellndorfer et al. [87] and Ludwig and Schneider [35] reported that thick forest vegetation can generate negative bias (i.e., overestimation) in SRTM elevations when comparing model-derived elevation to in situ-derived values. This can be explained by the relatively short wavelength (i.e., 5.6 cm) of the C-band SAR used for DTM generation, which does not penetrate the tree canopy completely as it is affected by (back-) scattering by branches and leaves within the vegetation canopy [80,88,89]. Our SRTM results, however, did not show any clear increased tendency to overestimate elevation values under forest cover as compared to open (agricultural and sparse vegetation) areas (Figure 2). This may be due to the low tree height and relatively open structure that characterise the Mediterranean dry formations that dominate the forest cover in the studied catchments. Similarly to the analysed ASTER data, any elevation bias produced by the Mediterranean dry forest formations distributed within the studied catchments in the SRTM data may be largely surpassed by the large (~7 m RMSE) intrinsic vertical errors of this dataset.

For the development of LiDAR-based elevation models, last returns are usually taken from discrete waveform LiDAR cloud data to ensure that backscattering occurs as close as possible to the bare surface, maximizing DTM vertical accuracy. However, severely filtered LiDAR returns are reported to affect DTM accuracy [83]. In this study, the density of preclassified ground returns of the LiDAR data underlying the analysed (5-m and 1-m) IGN DTMs was 1.65 points per m² in average for the three studied catchments, with the flat, open areas reaching up to 10 points per m² and the steep, vegetated hillslopes reaching the lowest values (generally below 1 point per m² on average). Data gap-filling techniques, such as the employed multilevel B-spline interpolation algorithm, are likely to produce artefacts in the production of high-resolution DTMs in the case that filling data gaps becomes necessary. Our results indicate that LiDAR DTM RMSE values slightly increased (by about 50 cm) for vegetated areas with low density of ground returns, as compared to flat, open areas (Table 2). Accordingly, previous studies indicated that gap filling can affect vegetated areas with low-density LiDAR ground returns, resulting in non-negligible errors in the resulting DTMs [85,86]. The analysed IGN LiDAR models showed a very similar relative distribution of elevation errors for both bare and vegetated areas (Figure 2). These results contrasted with Harding et al. [90], who observed higher underestimations in areas with dense vegetation compared to open terrain due to the dispersion of LiDAR pulses within the canopy. The results of the present work are more consistent with Clark et al. [83], who observed that the precision of LiDAR-based elevation models might vary where different patterns and vegetation types are combined with different slope gradients and morphological characteristics. Furthermore, Estornell et al. [82] reported that the effect of slope on the LiDAR DTM accuracy could largely surpass the direct effects caused by vegetation. In the studied Mediterranean catchments, the effect of vegetation on LiDAR-based DTM accuracy may be relatively small as compared to other error sources related to complex morphological features (e.g., terraces or artificial channels) and derived technical limitations related to filtering or interpolation processes in these complex areas.

5.1.3. Evaluating Terrain Morphology on DTM Accuracy

SRTM showed a low vertical accuracy at the mountainous Es Telègraf and the terraced Sa Font de la Vila catchments (Table 2), probably due to the presence of geometric distortions in the underlying C-band SAR imagery. Such distortions are caused by shadowing effects under extreme viewing conditions in very steep and fragmented areas [91,92]. These findings can also explain the observed straight-line stream features generated by the SRTM data in the mountainous Es Telègraf catchment (Figure 5), which are also likely to be related to severe intrinsic errors of the SRTM data in highly rugged, mountainous areas. Conversely, optical imagery is less affected by such relief distortions, as viewing angles are usually small. Nevertheless, the vertical accuracy of the ASTER model was extremely poor at the Telègraf catchment. Such poor performance can be directly related to the very low number of stack layers used for DTM generation, which was between 2 and 7 for the specific case of Es Telègraf. As previously discussed, the presence of large residual artefacts in the elevation data very importantly increases for areas where fewer than 10 stack layers were applied for ASTER GDEM production, inducing large errors and distortions in both steep [93] and low-gradient terrain [28]. Therefore, the insufficient amount of data available for the DTM generation largely conditioned the vertical accuracy of the ASTER data for this study, while the intrinsic errors of the underlying remotely sensed C-band SAR data represented for the SRTM dataset the main responsible factors affecting the poor performance of this synthetic aperture RADAR-based model in rugged terrain.

As mentioned above, a low density of ground returns for complex surfaces in the LiDAR data may be responsible for the observed small deviations in the vertical accuracy of the IGN DTMs. Model performance was relatively good for the mixed Es Fangar and the mountainous Es Telègraf catchment, although the data showed a poorer performance for the terraced Sa Font de la Vila catchment (Table 2). The relatively poor results obtained in Sa Font de la Vila are probably related to the highly fragmented topography that characterises this terraced catchment [94], which favours data oversampling in the flat areas of the terraces, but large data undersampling for the small-detail abrupt elevation changes aligned along the dry-stone walls of the terraces. This effect may be significantly increased by the presence of erroneously classified LiDAR returns in the datasets [82]. For the Es Telègraf and Es Fangar catchments, where the terrain is also complex, but less fragmented by human-made structures, these effects are also likely to be present, but to a lower extent.

5.2. Assessing Hydrogeomorphological Modelling Reliability

Hydrogeomorphological modelling reliability strongly depends on the ability of digital topographic models to accurately represent landscape morphology. In this study, very remarkable differences were identified between the descriptive statistics and hydrogeomorphological characteristics of the four assessed datasets.

5.2.1. Basic Terrain Attributes

Hypsometry did not show a high sensitivity to either DTM grid size or data source (Figure 4 and Table 4). Therefore, the ability of the DTMs to replicate general patterns of mass and energy stored within a landscape appears to be nearly independent of the horizontal resolution and dataset characteristics. The identified independencies are induced by the application of relative values of elevation and area that surpass the absolute differences [95]. Our results also showed consistent smoothing of the slope gradient values for the SRTM and ASTER datasets in all catchments (Table 3). The coarse spatial resolution of these models (i.e., 1 arcsecond) caused an oversimplification of the landscape elements, favouring the loss of hydrological and geomorphological details, which is in accordance with similar results obtained in other geographic regions [10,11,96]. Such coarse DTM resolution did not only simplified the landscape elements, but also affected subsequent catchment hydrogeomorphological modelling in terms of physical measures. In more detail, Zhang and Montgomery [8] reported the large implications of altered slope-gradient distributions on process-based hydrological modelling, as the laws of transport obtained from coarse-grained DTMs did not match those obtained from field studies, especially for small catchments.

Important differences in the distribution of the slope-gradient values were observed between datasets, showing a clear nonlinear relationship between the raster horizontal resolution and the derived terrain representation of hillslopes. Whilst the slope-gradient differences observed between the (1-m and 5-m) LiDAR models were tightly coupled to their differences in spatial resolution, slope-gradient differences identified between the (1-arcsecond) ASTER and SRTM datasets cannot be attributed to differences in horizontal grid resolution (Table 3). Thus, residual artefacts and errors associated with the intrinsic characteristics of these coarse-grained datasets were probably relevant not only in terms of vertical accuracy, but also in terms of landscape representation, such as slope gradient, stream network patterns, and catchment area delineation. In fact, these results were also supported by the large observed differences in flow routing characteristics (e.g., mean and SD flowlength) and catchment delineation of the coarse DTMs (Table 2 and Table 3).

5.2.2. Geomorphic Parameters and Relationships

In order to check the ability of the investigated DTMs to address hydrogeomorphological processes, different geomorphometric parameters and relationships were analysed, including the slope–area relationship and the cumulative area distribution (Figure 4) as well as the LS factor (Table 4). The obtained slope–area relationships and the cumulative area distributions for the satellite-derived SRTM and ASTER models reflected a poor landscape representation performance at the hillslope scale, and hence, a poor ability to assess soil erosion and flow aggregation processes. Contrarily, the LiDAR-based IGN models showed a good capacity to represent hillslope areas within the domain of diffusive and interrill erosion. This domain characterises effective catchment areas that contribute to the catchment conveyor belt [97]. Complete representation of hydrogeomorphological processes for meaningful assessment of catchment physical structure requires accurate DTM representation of dominant processes at different spatial and temporal scales. Thus, the SRTM and ASTER datasets provided insufficient information for hydrogeomorphological modelling in terms of water and sediment transfer processes in the studied small Mediterranean catchments. Both the SRTM and ASTER datasets showed very similar patterns in slope–area relationship and cumulative area distribution, whereas the 1-m IGN DTM induced changes in the position of critical (erosion-domain and flow aggregation) breaking points towards smaller values of contributing area (Figure 3). These results emphasise that the effect of the spatial resolution on the representation of different landscape features and erosion forms within the catchments was probably stronger than the effects caused by the intrinsic errors of the DTMs in the obtained geomorphological relationships. Differently to the slope–area relationship and the cumulative area distribution, the LS factor does not differentiate between regions or domains of erosion and flow aggregation within the catchments. A strong nonlinear dependency on raster resolution and derived slope gradient values was observed for our LS estimations, which largely agrees with the results obtained by Wu et al. [98], who assessed the influence of grid size on soil loss modelling in a mountainous catchment in Virginia (USA).

5.2.3. Stream Network Organisation

The residual artefacts and errors related to the intrinsic characteristics of the analysed datasets showed a strong effect on the organisation of stream networks (Figure 4). While differences in spatial resolution generally represent minor changes in the position of first-order streams [99], the residual artefacts and intrinsic errors of the studied satellite-derived SRTM and ASTER datasets induced the emergence of unrealistic fluvial networks (e.g., SRTM stream delineation for Es Telègraf; Figure 4a) and large uncertainties in stream network organisation and flowpath distribution (e.g., large ASTER deviations for Es Fangar; Figure 4b). Therefore, stream network delineation was directly affected by the propagation of errors in the input DTMs. Similar error propagation effects were reported by Oksanen and Sarjakoski [100], whereas Chaubey et al. [101] found that DTM grid size can significantly affect stream network and catchment delineation without invoking intrinsic elevation errors. Our analysis suggest that error propagation was more important than the grid size in this study, since stream network discrepancies between the LiDAR-based DTMs and the satellite-derived SRTM and ASTER models were relatively small for some of the study areas (e.g., Sa Font de la Vila catchment; Figure 4a).

5.2.4. Small-Detail Arrangement Patterns of Water and Sediment Fluxes

A relatively good consistency in the spatial pattern of modelled surface water and the sediment connectivity index (IC) was observed between the SRTM data and the two LiDAR-based models, while the patterns obtained from the ASTER data showed serious pitfalls for interpreting small-detail patterns of surface water and sediment flow. Since the ASTER imagery was collected with a ground-sampling distance of 15 m, small detail features remained within the subpixel range. Furthermore, the presence of large elevation inaccuracies and artefacts in the studied ASTER data strongly hindered the application of this dataset for detailed hydrogeomorphological analysis. Although more consistent within the SRTM data, the surface flow patterns obtained using this SAR-based DTM were considerably blurred, due to the coarse spatial resolution and lack of small-detail features of this dataset. In fact, the flow arrangement effects caused by microtopographic features such as concrete roads and preferential flow pathways were only captured within the 1-m LiDAR-based data (e.g., small detail plot at the mixed Es Fangar catchment; Figure 5), thereby highlighting the need for high-resolution (i.e., ≤5 m) DTMs in the application of small-detail hydrogeomorphological modelling.

Anthropogenic landscape features, such as small artificial drainage channels and ditches in agricultural fields, can significantly alter the hydrosedimentological functioning of traditional Mediterranean landscapes [102]. These artificial (metre- to submetre-wide) features were very likely underrepresented even for the 1-m LiDAR model. The coarser spatial resolution of the other analysed models made the detection of such traditional drainage systems largely unfeasible. Consequently, the model-derived stream networks and pathways of water and sediment connectivity obtained using the studied DTM sources approached a drainage pattern close to the “original” state of the system, but were unable to accurately reproduce the changes induced by human-made structures. This was probably true even for the analysed high-resolution 1-m LiDAR model. Therefore, the use of LiDAR-based models with very high resolution is recommended only if the underlying point density is high enough to reveal these small-detail features; otherwise, the consequent increases in computing time and required hardware resources cannot be justified in terms of increased accuracy and hydrogeomorphological modelling reliability.

6. Conclusions

This study examined the vertical accuracy of four DTMs (the 1-arcsecond SRTM C-SAR V3 and ASTER GDEM V2 datasets, as well as the 5- and 1-m LiDAR-based IGN models) with regard to data source, dataset characteristics, and terrain morphology in three small Mediterranean catchments. The reliability of these DTMs for subsequent hydrogeomorphological modelling purposes was also assessed. The results of this study can be summarised in the following keystone conclusions:

• The analysed LiDAR-based models and—to a lower extent—SRTM provided reliable sources for most of the discussed hydrological and geomorphological modelling aspects. However, SRTM notably failed to produce reliable results for highly rugged, mountainous areas due to intrinsic errors associated with topographic RADAR shadowing effects. For the analysed LiDAR-based models, attention should be also paid to the influence of data processing steps such as grid interpolation and point-cloud classification.
• ASTER showed the lowest vertical accuracy and considerable residual artefacts producing strong, non-normally distributed elevation errors that largely constrained the reliability of the ASTER elevation data. The presence of forest vegetation exacerbated the tendency of the ASTER dataset to overestimate elevation values (accounting for up to 30-m deviations), although the inherent, large vertical errors that affected this dataset largely surpassed the influence of Mediterranean dry forest vegetation in measured absolute vertical accuracy.
• Intrinsic errors and scarcity of the underlying DTM production data, vegetation patterns, and complex terrain morphology as well as relief fragmentation (especially for Mediterranean landscapes with traditional terraced structures) influenced the analysed datasets to different extents, resulting in significant deviations of elevation values.
• Both the vertical accuracy and horizontal resolution of the datasets were found to influence catchment hydrogeomorphological modelling in the studied sites. Error propagation impacted flow routing, stream network, and catchment delineation, and to a lower extent, the distribution of slope gradient values. Coarse horizontal raster resolution was found to reduce the degree of hydrological and geomorphological detail available from the DTMs and their reliability in representing processes at different spatial scales within the catchment.

The results presented in this study are transferable to other geographic regions dominated by fluvial processes and are not restricted to Mediterranean environments. Current trends in LiDAR-derived DTMs such as the use of unmanned aerial vehicles (UAV) and full-waveform LiDAR systems, as well as their possible gains in hydrogeomorphological modelling, should be addressed in future research.