Global digital elevation models for terrain morphology analysis in mountain environments: insights on Copernicus GLO-30 and ALOS AW3D30 for a large Alpine area

This study focuses on the quality evaluation of two of the best 1 arc-second public global digital elevation models (DEMs), Copernicus GLO-30 DEM and ALOS AW3D30 DSM, from the perspective of their capability to represent the terrain fine-scale morphology of a complex alpine landscape, located in the Italian Trentino Province. The analysis is performed on an area of 6210 km2, considering a reference DEM derived from a high resolution and accurate airborne Lidar survey. The quality assessment goes beyond a conventional approach based on elevation differences statistics, computed on a pixels-by-pixel basis. An ad hoc approach for evaluating the capability to represent fine-scale morphology, including surface roughness, is adopted. Moreover, the quality analysis is performed considering the influence of local morphology and of the different land covers. The findings show that although the two global DEMs have comparable overall quality, their relative performances change according to local landscape characteristics. Copernicus DEM performance is on average better than ALOS in correspondence of urbanized areas as well as in areas without vegetation cover, with gentle slopes and relatively low short-range roughness. Meanwhile, ALOS DEM performance is slightly better than Copernicus in rougher terrain and steeper slopes. In general, both DEMs have poor performances in steep slopes, with a limited capability to describe fine-scale morphology. The adoption of these global DEMs for terrain analysis and modelling of earth surface processes should be performed carefully, considering the impact of different land covers and of local morphology, including surface roughness.


Introduction
Global or quasi-global open-access digital elevation models (DEMs) produced from satellite data are marked by continuous improvements in quality and accuracy (Gesch 2018;Kakavas et al. 2020;Mudd 2020;. The most common and popular global DEMs include versions of SRTM DEM (Farr et al. 2007), ASTER GDEM (Abrams et al. 2020), ALOS AW3D30 DSM , MERIT DEM (Yamazaki et al. 2017), and Copernicus GLO-30 DEM (Airbus 2020;Strobl 2020). 1 These global DEMs can be further processed to improve their characteristics such as FABDEM (Hawker et al. 2022), which is derived from post-processing Copernicus GLO-30 DEM to filter out the effect of vegetation and buildings, adopting machine learning and various input features.
Global DEMs (Guth and Geoffroy 2021), generally more similar to digital surface models (DSMs) than to digital terrain models (DTMs), are frequently used in 1 3 198 Page 2 of 22 geoenvironmental and geoengineering studies. These models have been adopted for multiple research purposes, including the analysis and modeling of geoenvironmental processes related to terrain morphology more than to the overall surface morphology, which also includes anthropic structures and vegetation . Terrain morphology is fundamental in the context of geostructural, geomorphological, and geoenvironmental interpretation. In many applications, morphometric variables derived from DEMs are used as input features in supervised learning approaches or as proxies of other environmental variables (Wilson and Gallant 2000;Hengl and Reuter 2009;Florinsky 2016Minár et al. 2016;Trevisani and Florinsky 2021;. Examples of DEM-based studies include, but are not limited to, landslides susceptibility models (Meena and Nachappa 2019;Titti et al. 2021), derivation of shear wave seismic velocity of shallow subsoil (Wald and Allen 2007;Thompson et al. 2014;Heath et al. 2020), derivation of drainage networks and basins with a possible adoption of ad-hoc conditioning procedures of the drainage network in anthropized areas (Cavalli et al. 2013;Boulton and Stokes 2018). Another area of DEM-based research concerns geodiversity studies, where morphometric parameters contribute to derive indexes of geodiversity (Melelli 2014;Chrobak et al. 2021). Geodiversity studies covering large areas often employ global DEMs, although they focus on terrain morphology. DEMs are also fundamental for many processes concerning remote sensing methodologies such as SAR (Synthetic-Aperture Radar) interferometry technology, applied to the study of ground displacements (Bayer et al. 2017).
Given the importance of global DEMs (Schumann and Bates 2018), it is not surprising the large amount of studies concerning the quality and accuracy of these digital products (Purinton and Bookhagen 2017;Rizzoli et al. 2017;Boulton and Stokes 2018;Caglar et al. 2018;Florinsky et al. 2018;Gonzalez and Rizzoli 2018;Grohmann 2018;Jain et al. 2018;Li and Zhao 2018;Florinsky et al. 2019;Liu et al. 2019;González-Moradas and Viveen 2020;Kakavas et al. 2020;Vassilaki and Stamos 2020; Mesa-Mingorance and Ariza-López 2020; Guth and Geoffroy 2021;Hui et al. 2022). It is worth mentioning the recent digital elevation model intercomparison experiment (DEMIX, see Strobl et al. 2021), which aims at the development of a participative and standardized approach for global DEMs quality evaluation. The majority of studies focus on the evaluation of the inherent quality of DEMs considering a reference DEM or a set of ground control points ). The analysis is mostly performed on a pixel-by-pixels basis, that takes into account elevation differences (hereinafter referred to as EDs) between a global DEM and a reference, and sometimes takes into account geomorphometric derivatives, such as slope (Guth and Geoffroy 2021). The most advanced studies evaluate the quality of global DEMs considering different land covers and morphological settings. In fact, in areas with dense vegetation and anthropic structures the capability of global DEMs to reproduce terrain morphology is mainly influenced by the land cover characteristics, while in other areas it is mainly affected by the characteristics of local morphology, such as slope and surface roughness. However, the evaluation of global DEMs performance in reproducing finescale morphology is generally conducted by a visual analysis of the shaded relief, while quantitative approaches are rarely adopted (e.g., Florinsky et al. 2019;Polidori and El Hage 2020;Guth and Geoffroy 2021).
The capability of a DEM to reproduce terrain fine-scale morphology cannot be fully captured by the statistical analysis of EDs computed on a pixel-by-pixel basis, for example adopting conventional metrics such as root mean square error (RMSE, with the EDs considered as errors), mean/ median ED and mean/median absolute ED (e.g., Florinsky et al. 2019;Crema et al. 2020;Polidori and El Hage 2020). Otherwise, the local spatial-statistical variability structure (Isaaks and Srivastava 1989) of DEMs should be considered in the quantitative evaluation of DEMs quality, for example including geostatistical metrics linked to the fine-scale morphology, including surface roughness (e.g., Trevisani et al. 2012). The necessity to move beyond the evaluation of error statistics on a pixel-by-pixels basis is advocated also in the analogous context of image analysis. From this perspective, the structural similarity index (Wang et al. 2004;Crema et al. 2020) could be mentioned, since it considers the spatial variability structure in evaluating the similarity between a sample image and a reference one.
Accordingly, the main aim of this research is the quality evaluation of two of the best 1 arc-second public global digital elevation models, Copernicus GLO-30 version 2021_1 (hereinafter referred to as COP) (Airbus 2020;Strobl 2020) and ALOS AW3D30 version 3.2 (hereinafter referred to as ALOS) (JAXA 2021a; b), from the perspective of their capability to represent terrain fine-scale morphology of a complex alpine landscape, located in the Trentino Province (Italy). This study represents the development of a previous research (Florinsky et al. 2019) that analyzes the quality evaluation of several global DEMs in the same area, considering EDs with respect to a reference DEM (herewith referred to as reference), derived from the upscaling of an accurate and high-resolution airborne Lidar DTM (Cavalli et al. 2013;Portale geocartografico 2016). The cited study highlighted the necessity to evaluate quantitatively the capability to reproduce finescale morphology, for a more complete quality analysis of global DEMs. For this reason, in the current research a geostatistical approach is adopted to evaluate the performances of the selected global DEMs in reproducing fine-scale morphology. The study area is well suited for an analysis of this kind for several reasons: (1) large extent, approximately 6210 km 2 ; (2) high heterogeneity in morphology and land cover; (3) demanding morphological conditions for the remote sensing technologies adopted for the derivation of the global DEMs; (4) public availability of high-quality reference DEMs with creative common license (Portale geocartografico 2016).
The two global DEMs are seen among the best freely available global DEMs at 1 arc-second resolution. On the one hand, COP version 2021_1 seems particularly promising (Guth and Geoffroy 2021) with remarkable improvements respect to the older EU-DEM version 1.1 (Bashfield and Keim 2011) and other popular global DEMs. On the other hand, ALOS resulted the one with the best performances in the previous study conducted in the Alpine area (Florinsky et al. 2019) as well as in other regions (Caglar et al. 2018;Florinsky et al. 2018;Jain et al. 2018;Li and Zhao 2018;Liu et al. 2019;González-Moradas and Viveen 2020).
The quantitative analysis of the DEMs quality is performed both in terms of conventional EDs statistics, computed on a pixel-by-pixel basis, as well as from the perspective of the capability to represent fine-scale morphology (e.g., roughness). Moreover, the analysis is conducted taking into account specific land cover types and morphological settings.

Study area
The study area is situated in the Province of Trento, in the northeastern Italian Alps, with an extent of 6210 km 2 (Fig. 1). The morphology of the province is heterogeneous, also in relation to the complexity of the geostructural setting (Castellarin et al. 2005), with elevations ranging from 65 to 3760 m above sea level, with approximately the 70% of the area lying above 1000 m. The median slope is 26.5°, with 20% of the territory with slopes steeper than 37° (Table 1). Forests cover approximately 60% of the area, while rocky outcrops occupy about 7% of the territory. Flat and gentle slopes are intensely anthropized, including urban areas and agricultural land, that cover approximately the 19% of the territory. Lakes and glaciers cover minor proportions of the landscape (Cavalli et al. 2013).

Global DEMs
For the study area ( Fig. 1), ALOS AW3D30 DSM version 3.2 (JAXA 2021a) and Copernicus GLO-30 DEM version 2021_1 (ESA 2021) have been downloaded and cropped according to the study area extent.
The ALOS AW3D30 DSM is a photogrammetric DEM derived from visible-band stereo images with a 2.5 m resolution acquired by the ALOS spacecraft in 2006-2011 Takaku et al. 2016;Takaku and Tadono 2017). ALOS AW3D30 DSM versions 1.1 and 2.1 became respectively available in 2016 and 2018 (JAXA 2021a) at a resolution of 1 arc-second (approximately 30 m at the equator, and 26 m at the considered latitude) and a reported vertical accuracy (RMSE) of 5 m (Takaku and Tadono 2017). Further versions have been then released with quality improvements in northern and southern regions. In the high-latitude regions, different pixel spacing was adopted for each latitude zone. Coastlines were updated using new auxiliary data and removing incorrect sea masks over land. The anomaly detection algorithm has been improved and the anomalies were partially corrected with new additional data (JAXA 2021b).
The Copernicus GLO-30 DEM at a resolution of 1 arcsecond was obtained by resampling the original 12 m-gridded WorldDEM™ (Airbus 2020;Strobl 2020). The World-DEM™ has been derived from radar interferometry data (X-band) obtained during the TanDEM-X mission in 2011-2015 (Rizzoli et al. 2017), with a first release in 2019. The declared absolute vertical accuracy is approximately 4 m (90% linear error), the relative vertical accuracy is approximately 2 m for slopes ≤ 20%, and 4 m for slopes > 20%. The WorldDEM™ has been edited applying various filters and filling gaps using data from other global DEMs (Airbus 2020). In the studied spatial domain, the Copernicus DEM is characterized by several void filled areas resolved with a wide set of procedures and alternative DEMs (Fig. 2).
COP and ALOS DEMs have been projected in WGS 84 UTM 32N at a resolution of 25 m (the projected grid spacing along E-W and N-S directions is approximately 21.5 m and 30.9 m). Bilinear interpolation is used for the alignment (using Esri-ArcMap 10.8.1. GIS) according to a grid geometry based on the reference DEM. The decision to perform the analysis in a projected grid system and not in the original one of the global DEMs is dictated by the fact that geo-environmental applications, at the local and regional level (e.g., landslide susceptibility mapping, flood analysis, landslide runout modeling, etc.), are commonly conducted in projected coordinate systems and not in geographical ones. Accordingly, it is of interest to evaluate the quality of the global DEMs after the projection and alignment on a common grid system. Moreover, the alignment on the reference grid is convenient, considering that ALOS and COP have a different grid coding system, that induces a shift of half pixel between the two .
It is worth noting that the projection and alignment cause some degradation of the original quality of global DEMs; however, the error induced by the process is not marked and the impact on the statistical distribution of EDs is low. An alternative procedure preserving global DEMs grid geometry has been tested for COP, in order to derive an empirical evaluation of the impact of the followed procedure on EDs statistics. A subset of the study area, with an extent of approximately 12 km × 10 km, has been considered for the test. The alternative procedure maintains unaltered the COP coordinate system and grid geometry, hence preserving the original pixel values and accuracy. In fact, the reference DEM is derived according to the COP grid system, upscaling the original reference DEM at 2 m resolution. The upscaling is performed computing the average elevation of the 2 m pixels comprised in each pixel at 1 arc-second resolution. The two procedures provide a comparable statistical distribution in EDs between COP and the reference. The approach adopted in this study with respect to the alternative one leads to an increase of 5% for the median absolute ED and of 1.5% for the RMSE.
The geoid EGM2008 has been selected as the common potential gravimetric surface for the DEMs examined in this study. The elevation of COP is already expressed according to the geoid EGM2008 (Pavlis et al. 2012), which resembles the local geoid adopted in the Trentino Province. The elevations of ALOS are expressed according to EGM96 , less accurate and with lower spatial resolution than EGM2008. The transformation of ALOS elevations from EGM96 to EGM2008 has been conducted deploying the EGM96 deviations, derived from the 360 spherical harmonics components (code and data are available from https:// earth-info. nga. mil/), and the EGM2008 geoid grid at 1 arc-minutes resolution.

Reference DEM
The reference DEM at 25 m resolution ( Fig. 1) has been derived from high resolution DTMs obtained by means of airborne Lidar scanning (Cavalli et al. 2013). The Lidar survey has been conducted between 2006 and 2008 covering the Trentino Province (Portale geocartografico 2016). The DTMs have a resolution of 1 m and a reported vertical accuracy of 15 cm in valleys and urbanized areas; in the remaining areas the resolution is of 2 m with a reported vertical accuracy of 30 cm. The reference DEM at a resolution of 25 m has been generated by upscaling the original DEM resolution by cell aggregation using the mean as estimator.
To perform this operation the function "Aggregate" of the Arcmap10.8.1 Spatial Analyst extension has been deployed. Accordingly, the pixels values of the reference DEM represent elevation values with a spatial support of 25 m × 25 m.
As reported in the previous study (Florinsky et al. 2019), an indicative assessment of the accuracy of the upscaled reference DEM has been performed considering 111 ground control points, pertaining to a 4th-order regional geodetic network, stationed with satellite receivers using static or rapid-static techniques (Servizio catasto di Trento 2011;Chistè et al. 2013). The resulting mean error is of − 2.14 m with a standard deviation of 3.5 m; the error can be regarded relatively low, if one takes into account that the ground control points have a punctual spatial support (often located on peaks), and the elevations of the reference DEM are representative of the mean elevation of a 25 m × 25 m area. The reference DEM elevations are expressed according to a local geoid (Di Girolamo 2008) adopted for the topographic and cadastral services of the Trentino Province. The elevations have been transformed to the common geoid EGM2008, by means of the local geoid grid and the EGM2008 at 1 arc-minute resolution. From the tests performed on thousands of ground control points (Chistè et al. 2013), located mainly on structures, the interpolated local geoid is quite accurate, with a standard deviation of 0.2 m. Moreover, the differences between EGM2008 and the local geoid are not marked, with a mean difference of 0.47 m and a standard deviation of 0.07 m.

Fig. 3
Reclassified Corine land cover (CLC2006) for the study area. White rectangles show the approximate location of sites (S1-S9) represented in figures of "Materials and methods" and "Results and discussion"

Data analysis
The statistical and spatial-statistical analysis of the EDs between the global DEMs and the reference have been conducted using various software, that includes proprietary geographic information systems (ESRI ArcMap 10.8.1 and MapInfo PRO) and an open-source solution based on the Terra package (Hijmans et al. 2022) in the R statistical computing environment (R Foundation 2022). ArcMap 10.8.1 has been deployed for visual analysis, creation of maps, and cross-checking some of the results obtained in R. The basic statistical functions of R and the possibility to manipulate raster data provided by the Terra package have been deployed to compute the statistical metrics and statistical graphs.
Standard statistical indexes have been used to characterize the statistical distribution of EDs and absolute EDs. They have been computed pixel-by-pixel between the global DEMs and the reference: mean, median, standard deviation, RMSE, and quantiles. The Corine 2006 (CLC2006) land cover map (EEA 2020) has been adopted to stratify the analysis according to land cover typologies, being derived in the same period of the survey of the reference DEM. It should be mentioned that CLC2006 has a minimum mapping unit of 25 ha and, accordingly, some noise due to misclassification should be expected because of the higher resolution of the DEMs and the chronological discrepancies. Nevertheless, as confirmed also by means of the visual analysis of orthomosaics (pixel size of 1 m), the accuracy and resolution of CLC2006 permit to differentiate appropriately, for the intended tasks of this study, the different land covers in the area. Given the characteristics of the area and the study aims, the CLC2006 has been reclassified in a lower number of classes: 'Agriculture' (areas dedicated to agriculture), 'Forest' (dense vegetation cover), 'Forest (Low)' (transitional woodland-shrub), 'Grasslands' (mainly pastures and natural grasslands), 'Outcrop' (bare rock and shallow bedrock), 'Urban' (anthropized areas). The proportion of the selected land covers in the examined area are: 'Agriculture' 10.2%, 'Forest' 55.8%, 'Forest (low)' 6%, 'Grassland' 17.4%, 'Outcrop' 7%, 'Urban' 2.8%.
Metrics outlining the capability to reproduce fine-scale morphology have been considered to evaluate the global DEMs quality, in addition to the usual EDs statistics. The calculations have been mainly performed using focal functions of the Terra package and the standard statistical functions of R; ArcMap algorithms deploying Spatial Analyst focal functions (Trevisani and Rocca 2015) have been adopted for cross-checking the results obtained in R. Spatialstatistical indexes capable to provide information on local surface spatial-variability have been adopted (Trevisani et al. 2012). The methodology relies on a geostatistical-based approach for surface roughness analysis (Trevisani and Rocca 2015;Trevisani and Cavalli 2016).
In this study, the residual DEMs and a simple isotropic short-range roughness index are considered to characterize key aspects of fine-scale morphology. The residual DEM, analogous to the topographic position index (TPI, Guisan et al. 1999), is obtained by means of a high pass filter (e.g., Wilson and Gallant 2000;Hiller and Smith 2008;Florinsky 2016), removing large-scale variations (i.e., a smoothed DEM) from the original DEM. From a geostatistical perspective, the residual DEM can be interpreted as a detrended surface (e.g., Isaaks and Srivastava 1989), that highlights fine-scale morphology and permits to compute local roughness indexes, not influenced by local slope. In this study, the residual DEMs, for COP, ALOS and the reference, have been derived subtracting from the original corresponding DEM a smoothed DEM obtained via a simple single-pass moving average approach. A circular moving window with a radius of 3 pixels (i.e., covering 29 pixels) has been adopted. The size of the search window has been selected as a tradeoff between two contrasting requirements: on the one hand, the need to keep the window small to detect local details and, on the other hand, the necessity for a window large enough to contain a sufficient number of data (see Trevisani and Rocca 2015). The residual DEM highlights fine-scale morphologies including potential artifacts (Fig. 4), and it is fundamental for the calculation of the roughness index considered (Fig. 4c), even if solutions bypassing DEM detrending are available (e.g., Trevisani et al. 2023). Moreover, the analysis of the correlation (Pearson) between the residual global DEMs and the residual reference provides an immediate quantitative index that reports the capability to reproduce fine-scale morphology.
The study focuses on a specific aspect of surface roughness, namely the short-range isotropic roughness, computed by means of the geostatistical algorithm based on MAD (Median Absolute Differences) estimator (Eq. 1). MAD has been proposed by Trevisani and Rocca (2015) for the analysis of geomorphometric data as an improvement of the variogram (Eq. 2) or its generalization (Eq. 3) (e.g., Isaaks and Srivastava 1989): . 4 Shaded relief (a) and residual DEM (b) highlighting fine-scale morphology for site S1 (Fig. 3). The residual DEM is the input to calculate local short-range isotropic roughness (lag distance of 2 pixels) represented in c with Δ(h) = z u − z u + h In Eqs. (1-3), h is the separation vector (lag) between two locations (u), z(u) is the value of the variable of interest in the location u (e.g., residual elevation), and N(h) is the number of point pairs separated by h found in the search window. While the variogram (Eq. 2, or Eq. 3 with p = 2) and the madogram (Eq. 3 with p = 1) use as estimator of spatial variability the mean, MAD adopts the median, resulting less sensitive to data contamination. For a given search window, MAD reports the median absolute difference in residual elevation between pixels separated by a given lag distance and for a given direction. Calculating the indexes for different values of h it is possible to derive a multiscale and directional characterization of surface roughness.
In this study, a simple short-range isotropic (i.e. omnidirectional) index has been calculated according to Eq. (1), selecting a lag distance of 2 pixels in any direction and a circular search window with a radius of 3 pixels (Fig. 4c). The shortrange isotropic roughness index considered (for brevity, further referred as roughness), focusing on short-range spatial variability, provides an intuitive index of the capability of global DEMs to reproduce fine-scale terrain roughness. This roughness index was chosen for its simplicity and high informativeness; however, it is worth noting that other estimators can be adopted, for example tailored to the detection of artifacts in DEMs. From this perspective, variogram and madogram (Eqs. 2 and 3) would permit to highlight anomalous data, being these estimators more sensitive to outliers; however, the general description of surface roughness would appear blurred, given the heterogeneity in roughness of the studied terrain (see Trevisani and Rocca 2015). Moreover, other aspects and scales of surface roughness could be analyzed (e.g., roughness anisotropy, Trevisani et al. 2012). It is worth noting that the selection of scales and characteristics of roughness to be described as well as the inclusion of other geomorphometric derivatives (e.g., curvature indexes) is subjective and dependent on the objectives of the study.

Elevation differences in relation to land cover
The statistical distribution of EDs between the global DEMs and the reference ( Fig. 5 and Table 2), computed with all land cover typologies, reveals a higher positive bias and a higher dispersion of differences for COP than for ALOS. Even regarding absolute EDs, COP is characterized by higher median ED (7.1 m vs. 5.45 m), mean ED (10.06 m vs 8.34 m), and RMSE (16.09 m vs 14.01 m). However, when focusing on areas without vegetation and anthropic structures, i.e., "Outcrop" and "Grassland" land covers, the situation reverts, and COP is better than ALOS in terms of median ED and median absolute ED. For instance, in correspondence of "Outcrop" land cover, the median absolute ED for COP and ALOS is, respectively, 2.45 m and 4.2 m. The analysis of EDs for "Outcrop" and "Grassland" land covers provides information on the inherent quality of global DEMs, at least when adopting robust statistics such as the median, in order to reduce the impact of noise due to land cover misclassifications. It is evident ( Fig. 5 and Table 2) that in these land covers, COP is better than ALOS both in terms of median ED as well as of median absolute ED. ALOS shows a wider variability (Fig. 5) and in general a negative bias (e.g., median ED for grassland − 0.89 m). This is particularly evident in "Low Forest" and "Agriculture" land covers (Fig. 5), with negative values of the first quartile, that is unexpected given the presence of vegetation and considering that ALOS is a photogrammetric product. In general, COP is characterized by larger positive median differences for forested areas in comparison to ALOS. In "Urban" land cover the median differences are similar, but ALOS has a higher dispersion of differences (Fig. 5). As expected, the EDs between the global DEMs and the reference are strongly influenced by the land cover type (Fig. 5), with reflections also on the capability to represent fine-scale terrain morphology, as described in the "Quality in reproducing fine-scale terrain morphology". Given the statistics, the use of the global DEMs in forested areas should be performed carefully, especially when the purpose is to model surface flow processes or the derivation of terrain morphometric variables. The analysis of EDs in bare earth conditions reveals the negative bias of ALOS, that can justify its better performance in terms of EDs statistics, when looking at all land covers.

Elevation differences in relation to local morphology
Limiting the analysis to "Grassland" and "Outcrop" land covers (Figs. 6, 7), mostly representative of bare earth conditions, it is possible to highlight the influence of the local morphological setting on the EDs. The local morphology is characterized in terms of local slope and short-range roughness (Figs. 6, 7, 8), both computed on the reference DEM. The statistical metrics of EDs have been calculated stratifying the data according to classes of slope (class interval of 2.5°) and isotropic short-range roughness (class interval of 1 m). The chosen intervals permit a fine description still maintaining enough samples in each class. For each combination of slope-roughness classes, the following statistics have been calculated, considering a minimum of 50 data: number of data (i.e., pixels), mean slope of the reference, mean roughness of the reference, RMSE for COP, RMSE for ALOS, correlation between the residual global DEMs and the residual reference.
As general tendency (Figs. 6, 7), COP outperforms ALOS in areas with slopes lower than 40° and a relatively smooth topography. It should be noted that considering all land covers (Fig. 8) this behavior is not observed, because the statistical distribution of EDs is mostly influenced by land covers.
Limiting the analysis to "Grassland" and "Outcrop" land covers, it is interesting to analyze in more details the EDs between global DEMs and the reference, looking at areas with gentle slopes and low roughness. These settings are comparable to those adopted in the usual statistical performance tests of the global DEMs, which examine mainly flat areas (e.g., Caglar et al. 2018). Further, if the analysis is narrowed to slopes lower than 5° (for a statistical distribution of slopes see Table 1) and in any condition of local roughness, the RMSE of COP and ALOS is 4.54 m and 5.11 m, respectively. On the other hand, in areas with a roughness less than 1.3 m, corresponding to the median value of roughness in that slope range, the RMSE of COP and ALOS is 1.86 m and 2.62 m, respectively. It is worth noting that the impact of local roughness on EDs statistics is significant also in correspondence with mild and steep slopes. Up to slopes of approximately 40°, COP outperforms ALOS in areas with relatively smooth topography. For example, considering slopes between 30° and 35° in terrains with roughness values less than the median value of roughness (3.72 m) for that slope range, the RMSE of COP and ALOS is 4.86 m and 7.13 m, respectively. Above the roughness threshold, the RMSE of COP and ALOS is 15.85 m and 14.85 m, respectively.
At steeper slopes the accuracy of the ALOS and COP are comparable with slightly better performances of the former. For instance, with slopes between 50° and 55° in terrains with roughness values less than the median value of roughness (9.98 m), the RMSE of COP and ALOS is 22.5 m and 20.22 m, respectively. Above that roughness threshold, the RMSE of COP and ALOS is 48.44 m and 46.09 m, respectively. Anyway, at these very steep slopes the differences in accuracy between the two global DEMs are erratic, and it is not feasible to individuate a clear pattern, partly in relation to the low amount of data in this slope range (Table 1) and due to the interpolation error.  6 RMSE (m) for COP and ALOS in relation to slope and shortrange isotropic roughness, considering "Grassland" and "Outcrop" land covers. Contour lines report the number of data used for the cal-culation of the RMSE in each class of slope (class interval of 2.5°) and roughness (class interval of 1 m). The values have been interpolated on a regular grid for improving the appearance

Spatial distribution of elevation differences
A further appreciation of the EDs between the global DEMs and the reference can be provided analyzing their spatial distribution in areas with different land covers and morphology (Fig. 9). As confirmed by the statistics, in areas without vegetation cover and anthropic structures the differences between COP and the reference are not biased, while ALOS tends to have a negative bias, except for some areas with positive errors. A similar trend was observed for other European regions (Florinsky et al. 2018). In "Grassland" and "Outcrops" land covers, the variability of errors for COP seems lower than for ALOS. In urban areas, the patterns of positive EDs for COP resemble the Fig. 7 Ratio of the RMSE of COP to the RMSE of ALOS in relation to slope and isotropic short-range roughness, considering "Grassland" and "Outcrop" land covers. Contour lines report the number of data used for the calculation of the RMSE ratio in each class of slope (class interval of 2.5°) and roughness (class interval of 1 m). The values have been interpolated on a regular grid for improving the appearance Fig. 8 Ratio of the RMSE of COP to the RMSE of ALOS in relation to slope and isotropic short-range roughness, considering all land covers. Contour lines report the number of data used for the calcula-tion of the RMSE ratio in each class of slope (class interval of 2.5°) and roughness (class interval of 1 m). The values have been interpolated on a regular grid for improving the appearance Page 13 of 22 198 spatial distribution of anthropic structures and agriculture fields; differently, for ALOS the spatial distribution of EDs seems less spatially structured. These characteristics are then confirmed by the analysis of the correlation between the residual DEMs and by the capability to reproduce short-range roughness, discussed in the next sections.

Residual DEMs analysis
The visual comparison of the residual global DEMs (Fig. 10) with the residual reference highlights the poor quality in terms of fine-scale morphology representation and the presence of various artifacts.
A quantitative evaluation of the capability to represent fine-scale morphology is achieved calculating the correlation coefficient (Pearson) between the residual global DEMs and the residual reference (Table 3 and Fig. 11). In general, regarding all morphological settings and land covers, ALOS tends to show higher correlation than COP (0.72 vs. 0.65), even if it is characterized by the presence of various artifacts. Looking at the different land covers (Table 3), marked differences between the two global DEMs are in correspondence of "Outcrops", "Grassland" and "Forest" land covers. The fact that ALOS shows higher correlations in forested areas respect to COP (0.76 versus 0.68) is unexpected, considering that the methodology used for deriving COP should permit some penetration into vegetation (Guth and Geoffroy 2021) and consequently a better representation of underlying morphology. In contrast, the performance is comparable for "Agriculture" and "Urban" land covers, with the last one characterized by a slightly better performance of COP (0.66 vs. 0.65).
Focusing on "Outcrop" and "Grassland" land covers, it is possible to analyze as the correlation between the residual Global DEMs and the residual reference (Fig. 11) changes according to the local morphology, considering both slope and roughness, as reported in Sect. "Elevation differences in relation to local morphology". It is evident (Fig. 11) how the performances are influenced by the morphological setting, as already seen for the RMSE (Figs. 6, 7). COP outperforms ALOS in areas with low/moderate slopes (e.g., slope less than 35°-40°) and a relatively smooth surface (e.g., short-range roughness less than 8 m). ALOS has better performances with steeper slopes (up to 55°-60°) and higher short-range roughness.  Fig. 3). COP tends to be smoother than ALOS, however ALOS has various artifacts

Short-range roughness
The quality of fine-scale morphology represented by global DEMs can be further investigated evaluating how COP and ALOS reproduce short-range surface roughness (Figs. 12,13,14). The metric adopted is the ratio of the roughness computed from the Global DEM to the roughness computed from the reference. Both COP and ALOS smooth severely surface roughness, and the smoothing changes in function of the different land covers and of the local slope (Fig. 12).
This suggests that the real spatial support of measurement (Isaaks and Srivastava 1989;Burrough and McDonnell 1998) is much larger than 1 arc-second. COP, apart from "Urban" land cover, tends to smooth with all land covers, with an overall median reduction of short-range roughness up to 29%; ALOS is characterized by a lower level of smoothing, with a median reduction of 17%. However, the general lower smoothing observed for ALOS should be interpreted carefully, as it is related to multiple factors including land cover, morphological setting, and the Fig. 11 Correlation coefficients between the residual global DEMs and the residual reference DEM in relation to slope and short-range isotropic roughness (lag two pixels) for "Grassland" and "Outcrop" land covers. Contour lines report the number of data used for the calculation of the correlation in each class of slope (class interval of 2.5°) and roughness (class interval of 1 m). The values have been interpolated on a regular grid for improving the appearance presence of artifacts. If we focus on "Grassland" and "Outcrops" land covers (Fig. 12), it is evident a drop in the roughness ratios in correspondence of slopes steeper than 45°. In areas (Fig. 13) of high elevation and rugged topography, with a prevalence of "Outcrop" and "Grassland" land covers, the smoothing of roughness is marked, apart from localized sites where artifacts may increase the roughness.
Respect to COP, ALOS presents a higher frequency of areas with marked higher roughness than the reference (e.g., in the valley of Fig. 14) in "Agriculture" and "Urban" land covers. The smoothing of ALOS seems tendentially less severe than for COP, but this is partly due to the presence of artifacts (e.g., red boxes in Figs. 13,14,15) that contribute to increase local roughness.
The differences between COP and ALOS in their capability to reproduce terrain fine-scale morphology in the different morphological and land cover settings are partly justified by the different technology and methodology adopted for deriving these global DEMs. In fact, the specifics of signal acquisition and processing are of critical importance in the creation of the surface model. Even though attempts have been made to correct implausible values, high-altitude peaks in rugged terrain and/or with snow cover are still a challenge for terrain modeling (e.g., red box and arrows in Figs. 13, 15), both with the data obtained by radar imaging as well as Fig. 12 Box plots of the roughness ratios (ratio of global DEM roughness to reference roughness) for COP and ALOS. Top: roughness ratios versus land covers. Bottom: roughness ratios versus slope for "Grassland" and "Outcrop" land covers by optical range imaging. The radar surveying on high, rugged, and snowy terrain is often characterized by information gaps and void filling algorithms are applied. In the same terrains, photogrammetric processing of images generates artifacts caused by correlation errors. In urban areas (Fig. 14), the lower performance of ALOS respect to COP are likely caused by the interpolation of the photogrammetric point cloud without prior filtering of structures and vegetation. Finally, the stripy artifacts of ALOS (Figs. 14,15) are also related to the integer representation of elevation values. Regarding ALOS (Fig. 15), depressions and bumps are particularly evident and are likely related to correlation errors due to perspective distortions of the slopes, their varying conditions of illumination and, in some circumstances, presence of snow. The relevance of artifacts in the deterioration of DEM quality (Polidori and El Hage 2020) suggests the need to implement new quantitative indexes capable to characterize the different typologies of noise.

Conclusions
In the examined large Alpine area it is difficult to assess which of the two global DEMs, COP or ALOS, is the best from the perspective of terrain morphology representation. In fact, the relative performances are dependent on land covers and local morphological settings; accordingly, neither of the two DEMs is the best in all conditions. COP has on average better performances with low/mild slopes (up to 35°-40°) and relatively low short-range roughness, Fig. 13 Comparison of the global DEMs fine-scale morphology and short-range roughness with the reference DEM (site 3, Fig. 3). First row: shaded relief; second row: residual DEMs; third row: roughness for the reference DEM and roughness ratios (ratio of global DEMs roughness to reference roughness). Blue boxes: areas with different levels of smoothing of surface morphology. Red boxes: areas with artifacts and void filling processes especially in urbanized valleys and areas with a prevalence of agriculture and grasslands. ALOS shows slightly better performances in more rugged terrain and steeper slopes (up to 55°-60°); in regard to fine-scale morphology, ALOS has also better performances in correspondence of forested areas. In general, with steep slopes and rugged topography both DEMs have poor performances in reproducing finescale terrain morphology. In the studied area, especially in correspondence of gentle slopes with "Grassland" and "Agriculture" land covers, COP can be preferred to ALOS, even if it is smoother.
The analysis of residual topography and of short-range roughness permitted to highlight peculiar aspects on the quality of these DEMs related to local fine-scale morphology. Accordingly, this approach represents an additional tool to quantify specific aspects of DEM quality. Surface roughness is a relevant morphological factor, in conjunction with slope, influencing the accuracy of the DEMs. In addition, the analysis of residual DEMs and of the short-range roughness suggest that the true spatial resolution of these global DEMs is lower respect to their grid spacing and some level of upscaling may be required for their proper use. Finally, Fig. 14 Comparison of the global DEMs fine-scale morphology and short-range roughness with the reference DEM (site S2, Fig. 3). First row: shaded relief; second row: the residual DEMs; third row: rough-ness for the reference DEM and roughness ratios (ratio of global DEM roughness to reference roughness). The red box highlights an area with evident artifacts in ALOS the poor capability to reproduce fine-scale morphology, in specific morphological and land cover settings, represents a critical limitation in the adoption of these DEMs for geoenvironmental predictive/modelling tasks based on geomorphometric derivatives.
Author contributions S.T., T.N.S. and I.V.F. wrote the main manuscript text and reviewed the manuscript. S.T., and T.N.S. made the main statistical analyses related to the differences between global DEMs and the reference. S.T. made the analysis on residual DEMs and on Roughness. S.T. designed and created the figures presented in the paper with the review of T.N.S. and I.V.F..

Fig. 15
Shaded relief highlighting various artifacts in the ALOS and COP DEMs (Sites S7-S9, Fig. 3). Red arrows indicate areas with bumps or filled depressions; blue arrows indicate areas with artificial depressions; the red areas indicate the diffuse presence of artifacts Data availability There are no data to be provided. The data used for the paper are already available to public, because both the Global DEMs as well as the reference DEM can be dowloaded from the provided references.

Competing interests
The authors declare no competing interests.