Suitability assessment of global, continental and national digital elevation models for geomorphological analyses in Italy

Digital elevation models (DEMs) represent a fundamental resource in geomorphological analysis. The increasing availability of open‐access DEMs over wide areas is advantageous, but requires an evaluation of DEM quality and errors. This work applies a hierarchical assessment of global, continental and national DEMs in Italy in order to explore the differences and analyze the vertical accuracy, spatial error distribution and agreement of morphometric measurements. The selected DEMs are compared with local reference data as a ground points dataset, extracted from the national geodetic network, and regional DEMs at high spatial resolution, through both a qualitative and a quantitative approach. The results identify limits and potentialities of the selected DEMs, showing accuracy and errors in height representation, also affected by the topographic characteristics of the surface, such as steep slope in mountain zones and some defects in hydromorphological derivatives that could condition the geomorphological applications.

Here we focus on the exploitation of DEMs in geomorphological analysis based on a quantitative approach (i.e., geomorphometry). In order to identify and monitor landforms and the geomorphological processes that shape them, geomorphometry applies a numerical characterization of the surface and the detection of morphological variations. Therefore, a DEM represents a fundamental resource to observe Earth surface changes (Lane, Westaway, & Hicks, 2003), extract morphologic features and measure descriptors and parameters (Dixon, 1998;Rinaldi, Surian, Comiti, & Bussettini, 2016), and integrate topographic structures in the computation of geomorphic indices (Cavalli, Trevisani, Comiti, & Marchi, 2013;Rigon et al., 2013). During the last decades, the development of remote sensing led to the realization of DEMs at different resolutions. Mostly, the greater acquisition of satellite data determined the production of global and continental DEMs with ever-higher spatial resolution and coverage, greatly increasing the availability of these datasets over many areas of the Earth's surface (Baiocchi, Crespi, De Vendictis, Mazzoni, & Salerni, 2005;Capolongo, Marangu, Albanese, & Pennetta, 2008). Obviously, this aspect constitutes a considerable advancement in geomorphological analysis because it enhances the possibility to explore the land surface, ensuring improved access to digital models over wider areas. However, the use of DEMs in geomorphometry strictly depends on: (1) their accuracy in surface elevation; and (2) their sensitivity to topographic variability. In fact, the applicability of a DEM to derive geomorphological parameters needs to be verified according not only to altimetric accuracy, but also capacity to reproduce Earth surface complexity without the latter, in turn, affecting data acquisition techniques. Indeed, the various technologies used for surface detection (interferometric synthetic aperture radar-InSAR, photogrammetry, etc.) may suffer from biases and noise that depend on the physical structure (e.g., slope), cover (e.g., vegetation), and require additional processing to be removed (Schumann & Bates, 2018). Several studies evaluated global DEM accuracy in test areas (Alganci, Besol, & Sertel, 2018;Florinsky, Skrypitsyna, Trevisani, & Romaikin, 2019;Sertel, 2010;Uuemaa, Ahi, Montibeller, Muru, & Kmoch, 2020) or at nearly global scale (Kulp & Strauss, 2019;Rodriguez, Morris, & Belz, 2006;Wecklich et al., 2015;Wessel et al., 2018) through a comparison with ground truth and local reference data (ground control points-GCPs, GPS points, LiDAR, and local DEMs) by applying a statistical assessment and a spatial analysis of errors. These previous studies found different DEMs behavior related to the different surface characteristics such as topography, morphology, and land cover, highlighting that the available DEMs show discrepancies predominantly influenced by the spatial variability of Earth surface attributes such as the direction of slopes, steep and high mountains, terrain roughness, and dense vegetation (Miles, 2013;Mukherjee et al., 2013;Padova, 2013;Santillan & Makinano-Santillan, 2016;Zhao et al., 2011). This aspect becomes even more relevant when DEMs are used for geomorphological applications at regional and national scale. In Italy, elevation data and morphometric indices from open-source DEMs are integrated and correlated in quantitative analyses for the characterization and evolution of landscapes (Coco, Cestrone, & Buccolini, 2015;Martino, Nico, & Schiattarella, 2009), and in risk assessment studies for the realization of landslide susceptibility maps (Marchesini, Ardizzone, Alvioli, Rossi, & Guzzetti, 2014) and the definition of more exposed coastal areas to sea-level rise (Antonioli et al., 2017;Satta, Puddu, Venturini, & Giupponi, 2017).
It follows that the geomorphological analysis of large regions necessarily requires an assessment of the available DEMs to identify the advantages and drawbacks in geomorphic applications. Specifically, a focus on the quality of the opensource DEMs and their reliability in reproducing the characteristics and variability of the land surface in Italy is needed. Global DEM (ASTER); the Shuttle Radar Topography Mission 1 DEM (SRTM1); the European DEM v.1.1 (EUDEM); and the TINITALY DEM (TINITALY). DEMs are evaluated all over the Italian region, which is a very large area with considerable topographic and morphologic variability, not yet investigated in its overall territoriality. The main aim is to provide a suitability assessment of the elevation models for geomorphological analyses in such a large area that it requires the use of DEMs with wide spatial coverage (global, continental, or national DEMs). The study applies both qualitative and quantitative approaches to the evaluation steps. First, exploring the differences among the selected DEMs for the basic elevation information (i.e., vertical errors) and for derived surface attributes (i.e., slope angle and aspect). Second, comparing the selected DEMs with local reference data such as mark ground points from the national geodetic network and local high-resolution DEMs for three assessment levels: (1) statistical vertical accuracy; (2) spatial characterization of errors; and (3) agreement with morphometric measurements.
This study contributes to determine the reliability of global, continental, and national DEMs in accordance with the topographic characteristics and morphological features of Italy for addressing the most suitable use of models for geomorphological applications on the overall national territory.

| Italy
The Italian territory covers an area of 302,068 km 2 and comprises the northern region of the Alps, which connects it to the European continent, and the peninsular and island regions, which extend into the center of the Mediterranean Sea basin (Figure 1). Seas, lakes, lagoons, and rivers contribute to shape the physical landscape characterized by various coastal landforms (and a long coastline of about 7,500 km), glacial and fluvial valleys, mountainous reliefs and plains, depressions and caves, in heterogeneous sequences of morphological ups and downs (Bartolini, 2010).
The topographic variability distinguishes the Italian territory, with a prevalence of hill zones (42%) over mountain (35%) and flat (23%) areas. The elevation range varies from −3.44 (height of the lowest site in Le Contane locality, northeastern Italy) to 4,810 m a.s.l. (height of the highest site in Monte Bianco, northwestern Italy).
Different land cover classes define the surface of Italy, which overall can be divided into vegetated areas (forest, woodlands, agricultural trees, herbaceous, shrub, and pastures), crop fields, urban zones, and water bodies (Falcucci, Maiorano, & Boitani, 2006).
The Italian territory is divided into 20 administrative regions that identify distinct areas for geographical and topographic character, besides cultural, ethnographic, and socio-political aspects. This variability allows us to consider the regional unit as an additional test area to focus the analysis and evaluation of DEMs.

| Regional test areas
Three regions located in northern Italy (Piemonte), central/northern Italy (Emilia-Romagna), and southern Italy (Basilicata) (red areas and boxes on the right in Figure 1) are identified as test areas according to two main criteria: (1) the variability of the topographic characteristics that make these regional territories representative of the spatial and morphological complexity of Italy; and (2) the availability of open-source data.
The Piemonte region (top right box in Figure 1) has an average elevation of 421 m a.s.l. and a wide surface area of almost 25,300 km 2 . The mountainous area belongs to the Alpine range and it occupies 43% of the region, compared to 31% hill area and 26% flat area. A dense hydrographic network consists of creeks and minor and main channels, all tributaries of the Po River, which is the longest river in Italy (Carraro et al., 1995).
The Emilia-Romagna region (center right box in Figure 1) has an average elevation of 211 m a.s.l. and a surface area of almost 22,400 km 2 . The region is flat over 47.8% of the total area, due to the presence of Pianura Padana, the large floodplain of the Po River. The mountainous Apennines occupy 25.1% of the region, and the hill area extends over 27.1% of the territory. Emilia-Romagna is characterized by a low sandy coast and rivers that downslope the Apennines from inland, flowing into the Po River with a south-west/north-east direction (Castiglioni, Bondesan, Bondesan, Cavallin, & Gasperi, 1998). Figure 1) has a surface area of almost 10,000 km 2 and an average elevation of 633 m a.s.l. Mountain and hill areas predominantly occupy the territory, covering 47% and 45%, respectively; the flat zone represents 8% of the region with only one large plain, the Metaponto Plain, that extends along the Ionian coast in southern Basilicata. The main rivers, Bradano, Basento, Agri, Sinni, and Cavone (from north to south), flow from the Apennines range in a north-west/south-east direction (Lanorte et al., 2019).

| Global, continental, and national DEMs
Global, continental, and national DEMs generated by different types of data were selected for the suitability assessment, taking into account open-source models accessible for the whole area of the Italian territory. The F I G U R E 1 Study area and reference local dataset. Left: Italy with regional test areas (red areas) and GRPs dataset (blue points). Base map: ESRI Physical map. Right: DEM (hillshade) of regional test areas-Piemonte (top), Emilia-Romagna (center), Basilicata (bottom). GRPs, ground reference points rationale of the choice is to consider the resources available to apply geomorphological analysis in large Italian areas and/or over all the Italian region. The elevation models considered for this study, hereinafter named "selected DEMs," are digital surface models (DSMs).
AW3D30 was developed by the Japan Aerospace Exploration Agency (JAXA) through the photogrammetric technique applied on the optical stereo images (2.5 m resolution) acquired from the Panchromatic Remote-Sensing Instrument for Stereo Mapping (PRISM) on board the ALOS spacecraft in (JAXA, 2021. The spatial resolution of the dataset is approximately 30 m, which corresponds to a 1″ grid pixel. The elevation values are referred to the Earth Gravitational Model 1996 (EGM96, vertical datum) and the height accuracy of 4.40 m (root mean square error-RMSE) is confirmed by the quality assessment applied using GCPs distributed worldwide (Tadono et al., 2016). The 3.1 and 3.2 latest versions of AW3D30 (with global spatial coverage) can be downloaded in geographic World Geodetic System (WGS84) coordinates from the JAXA website free-of-charge (https://www. eorc.jaxa.jp/ALOS/en/aw3d3 0/). reference vertical datum is EGM96. According to Fujisada, Urai, and Iwasaki (2012), the height accuracy is 20 m at 95% confidence. The dataset is available in WGS84 coordinates (from 83° north to 83° south latitude) at no charge from Japan Space Systems (https://gdemdl.aster.jspac esyst ems.or.jp/index_en.html) and NASA Earthdata Search (https://search.earth data.nasa.gov/search) in its version 3, which has a decrease in elevation void area and an improvement in water body identification that led to an additional water-body dataset release. SRTM1 was developed jointly by NASA and the National Geospatial-Intelligence Agency (NGA). The DEM was generated through single-pass interferometry using data from different radar antennas (the C-Band Spaceborne Imaging Radar and the X-Band Synthetic Aperture Radar-X-SAR) aboard the space shuttle in 2000 (USGS, 2021). SRTM1 has approximate spatial resolution of 30 m, which corresponds to a 1″ grid posted. The heights refer to EGM96 geoid and the absolute height error is better than 9 m (Farr et al., 2007).
EarthExplorer of the U.S. Geological Survey (USGS) can be used for access to the collection data available (between 60° north and 56° south latitude) in the geographic WGS84 coordinates for the download (https:// earth explo rer.usgs.gov/).

EUDEM was developed by the European Environmental Agency (EEA) and managed by the European
Commission and the Directorate-General Enterprise and Industry (DG-ENTR). EUDEM was derived by fusing SRTM1 and ASTER data through a weighted averaging approach, including the results of the Copernicus program (Bashfield & Keim, 2011  TINITALY covers the whole Italian territory. It was developed by the National Institute of Geophysics and Volcanology (INGV) through the fusion of different data sources (contour lines and spot heights from technical regional cartography, LiDAR, and GPS data, orthophotos) (INGV, 2021). The DEM has a spatial resolution of 10 m and a mean RMSE value of 4.3 m (Tarquini et al., 2007(Tarquini et al., , 2012. The dataset is freely available in the UTM32N-WGS84 projection system as GeoTIFF tiles from the INGV website (http://tinit aly.pi.ingv.it/Downl oad_Area2.html).

| Reference local data
Different local reference data were used for the various assessment levels of the selected global and continental DEMs (see Section 4.2): the national geodetic network of the Italian Military Geographic Institute (IGM) and regional DEMs relative to the administrative territory of Piemonte, Emilia-Romagna, and Basilicata test areas (see Section 2.2).
The Italian national geodetic network was completed in 1995 and named the IGM95 network. It includes 5,000 ground points measured by GPS differential correction techniques through the Network Real Time Kinematic (NRTK) stations (https://www.igmi.org). The network is connected to the other geodetic national networks (triangulation and leveling networks). Each IGM95 point (here ground reference points-GRPs) is described by a local name, id number, geographic (from various coordinate systems) and topographic localization, altimetry with elevation derived from the ITALGEO 2005 Italian geoid (Albertella, Barzaghi, Carrion, & Maggi, 2008;Barzaghi, Borghi, Carrion, & Sona, 2007), and other information. The vertical accuracy of GRPs is about 5 cm (https://www.igmi. org/en/descr izion e-prodo tti/eleme nti-geode tici-1/rete-igm95). A set of 515 GRPs ( Figure 1) is extracted from the geodetic network, ensuring a homogeneous distribution over the Italian region, a representative allocation in each topographic zone, and a reasonable spatial coverage of the national territory. The set of GRPs is used as reference data for the comparison of point elevation with corresponding cell elevation from each selected DEM (see Section 4.2). A map visualization and an overall report about GRPs are available at the cadastral cartography website (http://carto grafia.globo gis.it/fiduc iali_gfmap let/).
The regional DEMs, realized by the corresponding regional authorities, have the same spatial resolution (5 m) and are open access.
The regional DEM of the Piemonte territory was obtained from LiDAR data acquired in 2009-2011 with vertical accuracy of about 0.3-0.6 m (Piemonte, 2021). The dataset can be downloaded as tiles in GeoTIFF format in the UTM32N-WGS84 projection datum from the Geo-Piemonte website (http://www.geopo rtale.piemo nte.it/ geoca talog orp/?sezio ne=catalogo).
The regional DEM of the Emilia-Romagna territory was developed by merging the altimetric data extracted from the technical regional cartography and LiDAR data (acquired in 2009) with a vertical accuracy of about 1 m (Emilia-Romagna, 2021). The DEM is available in the UTM32N-WGS84 projection datum from the Geo-website of the Emilia-Romagna region (https://geopo rtale.regio ne.emili a-romag na.it/downl oad/downl oad-data?type=raster).
The regional DEM of the Basilicata territory was realized using ground data from the photogrammetry and technical regional cartography with a vertical accuracy of 1 m (Basilicata, 2021). Through composing tiles in ASCII format, the dataset can be downloaded in the UTM33N-ETRS89 projection system from the regional WebGIS site (http://rsdi.regio ne.basil icata.it/webGi s2/gisDe vel.jsp?viewS tate=yes&proje ct=CCBB2 CC7-4D39-B4B4-3A07-72E3D E7F0443).

| Analysis of DEM differences
We analyzed the performance of the various available DEMs over the Italian region in order to distinguish the potential differences in the elevation data and in the main derived parameters, such as surface slope and slope direction (aspect) (Figure 2). For this evaluation step and the following comparison levels, the selected DEMs (AW3D30, ASTER, SRTM1, EUDEM, TINITALY) have to be considered in a common geodetic reference system (Alganci et al., 2018;Florinsky et al., 2019), which here is based on the UTM32N-WGS84 projected reference system and ~30 m spatial resolution. The rationale of this choice is to minimize the pre-processing of the dataset (i.e., geocoding, alignment, and resampling) in order to avoid overly affecting the original data, thus maximally exploiting the common spatial properties. Noting that Italy is divided into three UTM zones (i.e., 32N, 33N, and 34N), the use of the UTM32N projection is suggested by the common practice of adopting this UTM zone for the data that cover the whole national territory (national institutes such as INGV, Istituto Superiore per la Protezione e la Ricerca dell'Ambiente-ISPRA, and Istituto Nazionale di Statitica-ISTAT). The projection and resampling of all original raster data are performed by utilizing the GDAL geospatial software library. In particular, the gdal.Warp utility is applied for image warping using geometric transformation functions and nearest-neighbor resampling methods (GDAL/OGR Contributors, 2021).
After the pre-processing applied over the Italian region, the histograms of the selected DEMs are computed relative to three different elevation zones defined according to the homogeneous altimetric and geomorphological characteristics: flat (for areas with heights <300 m); hill (for areas with heights 300-800 m); and mountain (for areas with heights >800 m). Furthermore, a per-cell statistical analysis is applied by computing the value range corresponding to the difference between the largest and smallest height value in the stack of DEMs. After deriving the slope and aspect from each selected DEM, the same cell statistics from the multiple raster are applied in order to observe the potential discrepancies of DEMs in the morphometric products. The slope tool and the aspect tool are used to derive the spatial parameters from the raster surface (Spatial Analyst ArcToolbox; Burrough & McDonnell, 1998). Slope identifies the steepest downhill descent from the cell by computing the maximum change in elevation over the distance between the cell and its eight neighbors; it is measured in degrees. Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors; it is measured in degrees.

F I G U R E 2
Methodological approach for suitability assessment of DEMs. The selected AW3D30, ASTER GDEM, SRTM1 DEM, EUDEM, TINITALY DEM, and relative derived products are subjected to two investigation steps: analysis of differences (dashed box) and suitability assessment for geomorphological applications (dashed and dotted box). The latter is composed of three levels, each of them based on a different local reference dataset (circles), applied in various test areas (lozenges), from overall to specific accuracy evaluation target (square box on the right). See text for DEM acronyms

| DEM evaluation through three assessment levels
The suitability of the selected DEMs for geomorphological application is evaluated by comparing the DEMs with local reference data (i.e., assumed ground truth) through three assessment levels, which estimate the statistical vertical accuracy (level 1), the spatial characterization of errors (level 2), and the agreement of the morphometric products (level 3). This exploration by different levels is needed to move from an overall accuracy evaluation toward an increasingly specific investigation of the reliability of DEMs in morphometric applications ( Figure 2). Level 1. The set of GRPs (see Section 3.2) is used to estimate the vertical errors of the selected DEMs in Italy.
The GRP coordinates are converted into the same spatial system of DEMs (see Section 4.1) in order to extract the height value from the cell in which the reference point falls. Then the reference GRP value is subtracted from the model value: where E is the error, Z m is the height value of the selected DEM, and Z GRP is the height value of the GRP. The calculus is applied for each selected DEM.
For the calculus of the metrics (RMSE, standard deviation-SD, mean error, median, kurtosis, and skewness), the voids and outliers are excluded from the GRPs set by obtaining a common sample for the statistical analysis. Voids mean "return no value" from the cell in which the reference point falls. Outliers mean "differ significantly" from other error values. In order to identify outliers, the interquartile range (IQR) method is used (Tukey, 1977). Furthermore, the error distribution in 10 predefined classes of height difference is observed for all DEMs in order to examine overestimation and underestimation (positive and negative classes, respectively) of the height values. Error frequency histograms and relative cumulative frequency histograms are derived from this computation.
Level 2. The same analytical approach is applied by taking into account the different elevation zones (described in Section 4.1). The observation of the errors in areas with distinct topographic aspects (flat, hill, and mountain zones) contributes to correlate the performance of each DEM with the topographic conditions of Italy.
Level 3. Both qualitative and quantitative analyses are applied in regional test areas to compare the selected DEMs with local reference DEMs (see Sections 2.2 and 3.2). In order to make the comparison, the regional DEMs are subjected to spatial resolution reduction (from 5 to ~30 m) and spatial coordinate conversion in order to utilize the same reference system of the selected DEMs. First, topographic profiles are extracted from cross-sections in each elevation zone for all DEMs; the cross-sections are defined by setting a distance of 3,000 m and keeping the same metric pitch (about 30 km) passing through different elevation zones. The overlap of the height curves relative to each DEM (the selected DEMs and the regional DEM) allows a visual interpretation of the different height trends among the DEMs, identifying the topographic cross-sections along which the curves match and those along which the curves mostly deviate. Second, a quantitative approach is applied to compare DEM-derived products that are mostly used in geomorphological analysis, such as slope, aspect, and hydrographic network. After obtaining the raster of slope and aspect from all (selected and regional) DEMs in each test area (see Section 4.1 for slope and aspect computation method), the differences between values of regional derived products and values of products derived from each selected DEM are computed per cell: where P m is the (slope, aspect) value derived from the selected DEMs and P reg is the (slope, aspect) value from the regional DEM. The application of the formula generates a raster of differences of slope and aspect for each selected DEM (i.e., slope from AW3D30 minus slope from regional DEM, slope from ASTER minus slope from regional DEM, etc.). Statistics and plots of these rasters of differences are computed in order to observe the error trend by individualizing the elevation model that differs least from the regional DEM in the estimation of surface slope and aspect.
Furthermore, a fluvial basin is selected in each regional test area and the corresponding hydrographic network is derived from all DEMs. Rasters of flow direction and flow accumulation are obtained from DEMs (after sink filling); then the hydrographic network is defined by setting a flow accumulation threshold in order to extract channels (here a threshold of 1,000 identifies as stream channel cells those that have more than 1,000 pixels flowing into them).
Through the overlay of the hydrographic networks, the discrepancies in channel cell recognition and river flow mismatch are shown by identifying the most accurate elevation model in the description of the drainage network.
For the processing and analysis computations, the above cited GIS and calculus software is used (ArcMap™, QGIS®, Matlab®).

| DEM differences
The histograms of the selected DEMs show a similar frequency of height values over the Italian region. However, the overlap of the histograms in elevation zones ( Figure 3) shows a different behavior of ASTER and TINITALY (blue and orange line, respectively), deviating from the other DEMs. In particular, ASTER's height curve describes a distinct trend in the flat zone, with a higher pixel frequency value around 10 m. Also, spikes at constant elevation intervals characterize TINITALY, with greater evidence in hill and mountain areas.
The computation of cell statistics brings out a more relevant discrepancy among all DEMs. The range raster, which quantifies this value dissimilarity, shows that the hill and mountain zones are characterized by the greatest height differences among the DEMs, >10 to >30 m (green, orange, and red areas in Figure 4a).
The slope range (Figure 4b) visualizes moderate gradient discrepancy in hill and mountain zones (between about 5° and 30°); in the aspect range (Figure 4c), all elevation zones present moderate discrepancy in the gradient direction that shows greater differences in hill and flat zones (from 90° to 315°). It is specified that a range greater than 315° does not indicate high discrepancy in the gradient direction, considering that both 337° and 22° represent the north direction. This explains the first-class values set in the aspect range representation (Figure 4c).

| DEM assessment levels
For level 1 assessment, the initial dataset of 515 GRPs was used to compute the height differences [Equation (1)] of each selected DEM. A common sample of 435 error points obtained, excluding voids and outliers (see Section 4.2), was applied to calculate the metrics (Table 1). The results show that TINITALY has the lowest error values, with RMSE = 5.9 m, followed by AW3D30 and SRTM1 with RMSE = 6.5 m. Instead, the ASTER and EUDEM In the level 2 assessment, the statistical analysis of the same sample was applied in different elevation zones, with almost similar numbers of GRPs in mountain (133), hill (130), and flat (172) areas. The metrics reveal that all selected DEMs show lower accuracy in the mountain zone and higher accuracy in the flat zone (Table 2). In particular, RMSE values of TINITALY, AW3D30, SRTM1, ASTER, and EUDEM in the mountain zone correspond to 6.2, 7.4, 7.9, 10.3, and 12.8, respectively. Instead, in the flat zone, RMSE corresponds to 5.1 for SRTM1, 5.6 for TINITALY and AW3D30, and 6.9 for ASTER and EUDEM. This error trend is visible in error frequency histograms in which the error distribution in classes of height differences is similar in mountain and hill zones (Figures 5b,c), keeping the general overestimation and underestimation trends observed in the overall Italian region. Histograms relative to the flat zone (Figure 5d) show an improved performance of DEMs with the errors that are distributed in

RMSE (m)
SD ( The Level 3 assessment compared the selected DEMs with the regional DEMs in test areas (see Section 2.2). Figures 6, 7, and 8 show topographic profiles extracted from cross-sections in each elevation zone in Piemonte,

F I G U R E 5
Graphs of error distribution of the selected DEMs in predefined classes of height differences: histograms of error frequency (left) and corresponding error cumulative frequency (right). (a) Histograms relative to the overall Italian region (assessment level 1). (b-d) Histograms relative to the different elevation zonesmountain, hill, and flat, respectively (assessment level 2) Emilia-Romagna, and Basilicata, respectively. There are differences among the profiles of DEMs relative to the cross-section of the flat zone in the Piemonte region; in particular, the height plots of EUDEM and ASTER have a different trend that becomes reversed compared to the profile of the regional DEM in some traits of the transect (detailed box of profile 1 in Figure 6). In hill and mountain zones, the profiles of the selected DEMs are close to the reference DEM, but shift in height along the cross-sections (see detailed box in Figure 6). This overall correspondence with the regional DEM in hill and mountain zones can also be observed in the Emilia-Romagna region (Figure 7), where the profiles relative to the cross-section of the flat zone show an evident mismatch of ASTER compared to the others (profile 1 and detailed box in Figure 7). A focus on the profiles in the Emilia-Romagna region highlights that TINITALY deviates in some points along the cross-sections (detailed boxes in Figure 7).
The same aspect can be recognized by observing the profiles in the Basilicata region (Figure 8), where the trend of TINITALY and EUDEM curves does not coincide with the trend of the regional DEM in flat and hill zones.
Moreover, the profile of SRTM1 in the hill zone and the profile of ASTER in the mountain zone show a difference in height values compared to the regional DEM.
The correspondence between the products (slope and aspect) derived from the selected DEMs and those derived from the regional DEMs was evaluated by applying Equation (2) in each regional test area. Rasters of differences were analyzed by computing statistics (Tables 3, 4, and 5 corresponding to Piemonte, Emilia-Romagna, and Basilicata, respectively) and by overlapping the relative plots (Figures 9a,b, 10a,b, 11a EUDEM and TINITALY for slope differences in this region. In general, histogram curves describe differences less than ±5° for slope and less than ±50° for aspect, with a minimum tendency to overestimation in slope computation detectable in Piemonte and Emilia-Romagna, and underestimation detectable in Basilicata.

F I G U R E 6
Topographic profiles of the regional and selected DEMs in flat (Profile 1), hill (Profile 2), and mountain (Profile 3) zones of Piemonte test area (below). Square boxes correspond to detailed view of profiles The hydrographic network of Tanaro River in Piemonte (Figure 9c), Taro River in Emilia Romagna (Figure 10c), and Basento River in Basilicata (Figure 11c) is derived from regional and selected DEMs. The overlay of the networks highlights great differences among DEMs. While in upslope areas there is a minimum mismatch more visible F I G U R E 7 Topographic profiles of the regional and selected DEMs in flat (Profile 1), hill (Profile 2), and mountain (Profile 3) zones of Emilia-Romagna test area (below). Square boxes correspond to detailed view of profiles in the Piemonte region, in downslope areas all hydrographic networks deviate greatly from the network derived from the regional DEM. Besides TINITALY, which seems not to recognize channel cells in downslope areas of the Taro basin (Emilia-Romagna), SRTM1 and EUDEM also show differences in drainage cell estimation that determine F I G U R E 8 Topographic profiles of the regional and selected DEMs in flat (Profile 1), hill (Profile 2), and mountain (Profile 3) zones of Basilicata test area (below). Square boxes correspond to detailed view of profiles a negligible shift of channels in upslope areas and considerable direction changes of channels in downslope areas in all fluvial basins. A detailed observation suggests that AW3D30 defines similar hydrographic networks compared to the regional DEMs in all test areas.
The results relative to the third assessment level indicate that the extraction of some geomorphological features from the selected DEMs can generate different products, which do not always agree with local reference data.

| D ISCUSS I ON
Global, continental, and national DEMs were analyzed on the overall Italian territory and compared with local reference datasets to observe differences among the elevation models that derived from various types of TA B L E 3 Statistics corresponding to the differences of slope and aspect from the selected DEMs with slope and aspect from Piemonte DEM

Difference of derived products from regional DEM a SD (°) Mean (°) Median (°) Kurtosis Skewness
Slope AW3D30  sources, and to assess their reliability starting from the vertical accuracy, even considering the topographic variability, to obtain the basic parameters usually extracted from DEMs and useful for geomorphological analysis.
AW3D30, ASTER, SRTM1, EUDEM, and TINITALY express different height values over the Italian region, and in particular in mountain areas where the height difference range reaches and exceeds 30 m (Figure 4). This dissimilarity in surface elevation is also shown by the asymmetrical trend of ASTER's height curve in the flat zone, where the shift in value peak from 0 to around 10 m of ASTER's curve with respect to the other DEM curves is likely localizable in coastal plains, as already observed in the previous studies (Grohmann, 2018;Sertel, 2010).
The different method of DEM generation could explain the distinct reproduction of surface elevation, further conditioned by the geomorphology of detected areas. For example, shadowing effects caused by steep slopes can influence interferometry processing, and homogeneous surfaces can affect the photogrammetric technique due to more difficult feature identification. Certainly, the characteristics of source data also play a role in surface modeling, thus producing DEMs with different behavior-albeit derived from the same methodology and referred to similar topographic and morphologic surfaces. Therefore, the higher spatial resolution of optical stereo images could differentiate AW3D30 from ASTER in coastal plains detection, without height value shift occurrence. Also, the use of contour lines among input data of TINITALY could determine artifacts visible in spikes at constant distance in all elevation zones (Capolongo et al., 2008).
Different RMSE values obtained in assessment levels 1 and 2 (Tables 1 and 2) confirm the dissimilarity among DEMs and distinguish TINITALY, AW3D30, and SRTM1 as more accurate elevation models using the GRPs dataset as benchmark (though the RMSE is greater than the RMSE declared by the producers, except for SRTM1). This result is explainable by considering that: (1) TINITALY is a national DEM derived from the fusion of local data at high spatial resolution (see Section 3.1), strongly controlled and supervised in order to correct errors and avoid anomalies (Tarquini et al., 2007); and (2) AW3D30 and SRTM1 are global datasets with a fair vertical accuracy, as described by relatively small residual errors and supported by the previous studies (Alganci et al., 2018;Uuemaa et al., 2020). In contrast, EUDEM (with RMSE that disagrees with the officially stated worldwide accuracy) and ASTER show vertical errors that could condition their quality and potential applicability, above all in the mountain zones (Florinsky et al., 2019;Reuter, Neison, Strobl, Mehl, & Jarvis, 2009). Furthermore, the metrics, and kurtosis values in particular, suggest the tendency to produce more extreme errors for SRTM1 and EUDEM than AW3D30, ASTER, and TINITALY, which approximate the normal distribution of errors.
TA B L E 5 Statistics corresponding to the differences in slope and aspect from the selected DEMs with slope and aspect from Basilicata DEM By observing the spatial characterization of errors in level 2, the vertical quality of all selected DEMs in Italy appears affected by the topography (Table 2 and Figure 5). In fact, frequent errors are attested in areas with rough terrain and high relief; the greatest accuracy in height representation is denoted in flat zones. The influence of variations in topography on the vertical quality of DEMs resulted from other areas in previous studies (Falorni, Teles, Vivoni, Bras, & Amaratunga, 2005;Miles, 2013;Mukherjee et al., 2013). Generally, the trend of error distribution in classes of height differences describes an overall overestimation of height values by TINITALY, AW3D30, and SRTM1 in all analyzed elevation zones of Italy, an overall underestimation F I G U R E 9 Comparison of derived products (slope, aspect, hydrographic network) from the selected DEMs and regional Piemonte DEM. (a) Histogram of the difference relative to the slope. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±30° for a better visualization of the results. (b) Histogram of the difference relative to the aspect. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±300° for a better visualization of the results. (c) Overlay of hydrographic networks in upslope (yellow box) and downslope (green box) areas of Tanaro River catchment by EUDEM except in the flat zone, and both overestimation and underestimation by ASTER. The positive and negative errors of DEMs could be influenced not only by the topographic characteristics, but also by the surface cover (Padova, 2013;Santillan & Makinano-Santillan, 2016). The overestimation, for example, may be due to high vegetation, such as forest and woodland, in the mountain zone and/or the presence of buildings. In order to better investigate this aspect, we applied an analysis of the largest errors by spotting the location and corresponding land cover with height difference equal to or greater than ±10 m. The results confirm the majority of the largest errors in mountain and hill zones, but do not show an evident correlation between positive/ negative errors and land cover classes. In fact, the largest DEM errors can be found in urban areas, along the F I G U R E 1 0 Comparison of derived products (slope, aspect, hydrographic network) from the selected DEMs and regional Emilia-Romagna DEM. (a) Histogram of the difference relative to the slope. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±30° for a better visualization of the results. (b) Histogram of the difference relative to the aspect. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±300° for a better visualization of the results. (c) Overlay of hydrographic networks in upslope (yellow box) and downslope (green box) areas of Taro River catchment roads, in forest and woodlands, in areas with pastures and herbaceous vegetation. Furthermore, the analyses of lower errors with height difference equal to or less than ±3 m fall in the same land cover classes, with the only exception in crop fields that seem to characterize surfaces well detected by DEMs. Therefore, we could argue that the topographic variability plays a significant role compared to the land cover type in influencing the elevation errors of DEMs.
Although this analysis gives some indication of the vertical accuracy of DEMs, it is based on spot height differences that strongly depend on the spatial location of reference points and the horizontal accuracy of DEMs.

F I G U R E 11
Comparison of derived products (slope, aspect, hydrographic network) from the selected DEMs and regional Basilicata DEM. (a) Histogram of the difference relative to the slope. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±30° for a better visualization of the results. (b) Histogram of the difference relative to the aspect. Each color curve describes the difference of the selected DEM minus the regional DEM. The range of x-axis values was fixed to ±300° for a better visualization of the results. (c) Overlay of hydrographic networks in upslope (yellow box) and downslope (green box) areas of Basento River catchment Regarding this aspect, we specify that the assessment of levels 1 and 2 was also applied in the 33N zone of the UTM projection, obtaining results with a negligible difference.
The visual observation of the topographic profiles in regional test areas confirms some results and introduces new evaluation elements. In fact, if the profiles of EUDEM and ASTER mostly, and SRTM1 to a lesser extent, show small height oscillations along some cross-sections with respect to the reference DEM, deviations of TINITALY in the profiles of the Emilia-Romagna and Basilicata regions suggest a potential deficiency of this DEM in some areas (Figures 6-8), visible also by deriving the hydrographic network in the Tanaro (Piemonte) and Taro (Emilia-Romagna) basins (Figures 9c and 10c, respectively). Probably, this behavior is connected to the basic interpolation technique of TINITALY, which can detect the density and spatial resolution of source data, and improves as density and spatial resolution increase.
The quantitative assessment of derived products in regional test areas allows us to formulate some evaluations about slope and hydrographic network extraction. Generally, slope computation from DEMs at low spatial resolution introduces inclination errors that are visible by observing the curves of differences with regional DEMs in all test areas (major frequency corresponding to ±5°; Figures 9a, 10a, and 11a). These errors do not inhibit the slope map production (as described by the trend of curves), but can determine inaccuracies, more detectable in flat zones (Grohmann, 2018;Mouratidis, Karadimou, & Ampatzidis, 2017;Zhao et al., 2011), as proved by the histogram relative to the Emilia-Romagna region ( Figure 10a). Here, where the flat areas occupy 47.8% of the total territory (see Section 2.2), the curve of slope differences shows more oscillations (see Section 5.2), identifiable with artifacts. Clearly, the latter are associated with the DEM production process and depend on the spatial distribution of DEM height errors (Albani & Klinkenberg, 2003), less traceable through spot height assessment (levels 1 and 2), and more evident through derivatives analysis (level 3). Obviously, inaccuracies in the slope gradient affect the definition of flow direction and thus the identification of the hydrographic network, causing stream shift and channel deviations in downslope areas more than in upslope areas. The procedure of drainage network extraction from the selected DEMs in the Tanaro, Taro, and Basento basins (Figures 9c, 10c, and 11c) leads to this assessment scheme, by underlying the extreme impact of applying DEMs in hydrogeomorphic analysis and further focusing on the need to investigate more deeply DEMs' suitability for geomorphological applications. The present study contributes to observing limits and advantages in surface modeling of the global, continental, and national DEMs, by considering the potential effects of the topographic and morphologic characteristics of the overall Italian territory and regional test areas on DEMs' quality and applicability.
AW3D30 and SRTM1 demonstrate a good trend in surface representation, with minimal height accuracy in mountain zones and some defects in slope gradient estimation for plains. The high spatial resolution of the source data (optical stereo images in panchromatic band) and the high accuracy of the source technical approach (SAR interferometry), respectively, advantageously influence the reproduction of surface attributes, though the lower spatial resolution of the final product (~30 m) surely conditions the quality of corresponding derivatives.
ASTER and EUDEM contain height errors that produce some anomalies in both height accuracy and surface shape representation. In particular, ASTER presents atypical behavior along coastal flat areas and a general overestimation and underestimation of elevation values, and EUDEM proves not to enhance the quality aspect of ASTER and SRTM1, from which it derives. In fact, this hybrid product seems to smooth the discrepancies of the two global DEMs on the one hand, and keep the defects and artifacts on the other hand, in Italy as well as in other regions (Florinsky et al., 2019;Mouratidis et al., 2017). However, it should be considered that in the definition of the common spatial system for the analysis and comparison of the present work, EUDEM is the most affected DEM in its original spatial characteristics (reference system and resolution).
TINITALY configures as a national DEM available at relatively high spatial resolution (~10 m) that depends on the high spatial accuracy of local input data used for DEM generation. This advantageous aspect is retained in the height accuracy of the DEM in all elevation zones; however, the DEM production technique limits the quality of surface modeling, determining some anomalies and discrepancies in derived surface attributes. However, it should be considered that the definition of the common spatial system for the analysis and comparison of the present work might have influenced the results obtained for this DEM with respect to its initial spatial resolution (from 10 to 30 m).
This assessment could be useful in applying global, continental, and national DEMs in the overall Italian region for geomorphological analyses, in order to exploit the potential of open-source DEMs in such a large area, by improving the awareness of a more suitable use.

| CON CLUS IONS
Our analysis explored the use of available DEMs in a large area such as Italy in order to take into account the effects of DEM application in geomorphological analyses. The results show that in the case of the Italian territory: 1. TINITALY, AW3D30, and SRTM1 appear more accurate in height values than ASTER and EUDEM.
2. The topography affects the elevation representation in the DEMs that generally show different spatial pattern of errors, with the worst trend in mountain zones.
3. The morphologic characteristics of surfaces may not be well reproduced in the derivatives (such as slope, aspect, and hydrographic network) because of some defects and artifacts of the DEMs, which can be used with no particular restrictions but with a further investigation of conditions and effects on the geomorphological applications.
The assessment of DEMs over such a large area was made difficult in this work because of the complexity in terms of spatial characteristics and data availability. These elements represent limitations of the present study that force choices such as the selection of DEMs, the definition of the common spatial system for the comparison of different DEMs, the research of national local data for the vertical accuracy evaluation, and the identification of the regional test areas partly conditioned by the open-source local data.

ACK N OWLED G EM ENTS
The contribution of M.L.S., R.C, R.R. and P.P. was supported by the funding of PhD course in Geosciences.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study, derived from resources available in the public domain (list of open source DEMs and URLs can be found in the text), are available from the authors upon request.