Parameter-free delineation of slope units and terrain subdivision of Italy

Abstract Quantitative geomorphological and environmental analysis requires the adoption of well–defined spatial domains as basic mapping units. They provide local boundaries to aggregate environmental and morphometric variables and to perform calculations, thus they identify the spatial scale of the analysis. Grid cells, typically aligned with a digital elevation model, are the standard mapping unit choice. A wiser choice is represented by slope units, irregular terrain partitions delimited by drainage and divide lines that maximise geomorphological homogeneity within each unit and geomorphological heterogeneity between neighbouring units. Adoption of slope units has the advantage of enforcing a strong relation with the underlying topography, absent in grid cell–based analyses, but their objective delineation is still a challenge. A given study area admits different slope unit maps differing in number and size of units. Here, we devise an objective optimisation procedure for slope units, suitable for study areas of arbitrarily large size and with varying terrain heterogeneity. We applied the new approach to the whole of Italy, resulting in a map containing about 330,000 slope unit polygons of different sizes and shapes. The method is parameter–free due to objective optimisation using a morphometric segmentation function, and the map is readily available for general–purpose studies. A cluster analysis of slope units properties, compared with terrain elevation, slope, drainage density and lithology, confirmed that the terrain partition is geomorphologically sound. We suggest the use of the slope unit map for different terrain zonations, including landslide susceptibility modelling, hydrological and erosion modelling, geo–environmental, ecological, forestry, agriculture and land use/land cover studies requiring the identification of homogeneous terrain domains facing distinct directions.


Introduction
We investigate the automatic subdivision of a large and complex landscape into slope units (SUs) i.e., terrain units delimited by drainage and divide lines, using exclusively information obtained from a digital elevation model (DEM). Slope units are a valid alternative to grid cells as mapping units. They bear a strong relationship with topography, at variance with grid cells, and are potentially of interest in wide range of applications. They help dealing with heterogeneous data, solving problems stemming from positional accuracy and aggregating distributed data in a meaningful way.
In this work, to delineate SUs we adopt the first approach using the r.slopeunits software, a GRASS GIS module introduced by Alvioli et al. (2016) that adopts drainage and divide lines extracted from a DEM. The software implements an algorithm that considers landscape hydrology, the variability of terrain aspect, and geometric constraints for SUs (Alvioli et al., 2016). In an iterative process, the software defines SUs as "half-basins" using the hydrological model implemented in the GRASS GIS r.watershed module (Metz et al., 2011), which partitions a digital topography into drainage basins using the upslope contributing area, T. Using different T values at different epochs of the iterative process, the software adapts to the local morphology and produces SUs subdivisions with different spatial detail in different parts of a complex landscape. The software does not determine univocally the characteristic scale of the SUs partition. Instead, the level of detail of the SUs subdivision in different areas depends on a few software parameters, which need to be optimised in an objective and geomorphologically sound way. Working in the ≈2000 km 2 , Upper Tiber River Basin, Central Italy, Alvioli et al. (2016) optimised the parameters using objective functions. Here, we show that optimisation can be achieved independently from the geographical extent and the complexity of the study area.
The paper is organised as follows. After a brief summary of the background of the study (Section 2), we describe the methods used to prepare the SUs subdivision for Italy (Section 3). This includes the local optimisation procedure, and how the procedure is extended to areas of arbitrarily large size. Next, we present the results of the application of the method to the entire Italian territory (Section 4), and we discuss the results obtained (Section 5). We conclude summarising the lessons learnt (Section 6).

Background
Partitioning a complex landscape or a large territory into a mosaic of SUs is a form of terrain classification. A typical elementary unit for terrain classification is a "landform" (Minár and Evans, 2008;Evans, 2012), and a number of studies have defined landforms and investigated their topographic, geomorphological or thematic characteristics (e.g., Guzzetti and Reichenbach (1994); Drăguţ and Blaschke (2006); Drăguţ et al. (2011); Jasiewicz and Stepinski (2013); Dekavalla and Argialas (2017); Iwahashi et al. (2018); Luo and Liu (2018)). Landform classification is based typically on an a piori description of agenerally limitednumber of elementary morphological shapes, most commonly described in terms of discontinuities in terrain elevation, slope, aspect or curvature (Evans, 2006(Evans, , 2012. As an example, Jasiewicz and Stepinski (2013) defined "geomorphons", elementary morphological patterns singled out from a large set of theoretical patterns describing landform types of real natural landscapes.
A common problem in terrain classification is the selection of the characteristic scale of the terrain subdivision (Evans, 2003;Saito et al., 2011;Drăguţ et al., 2011;Drăguţ and Eisank, 2012;Ehsani and Quiel, 2014;Alvioli et al., 2016;Dekavalla and Argialas, 2017;Schlögel et al., 2018). Although examples exist of landforms that exhibit allometric scaling (Evans, 2012), partitioning of a complex landscape into SUs is not a scale invariant process, and it requires an appropriate, and not uniform, basin size for the drainage and divide networks.
Identifying the correct scale of a particular spatial analysis is a way out from what is known as modifiable areal unit problem (Openshaw (1984); Manley (2014)). Any study associated with the use of data aggregated within geographical areas is prone to the MAUP, and a objective link between mapping units and the underlying topography is highly desirable.
One way of addressing the characteristic scale issue is to use adaptive algorithms. For the selection of their "geomorphons", Jasiewicz and Stepinski (2013) did not decide the size of their landforms a priori; instead, they used a classification algorithm that adapts to the underlying topography. A limitation of the approach was that the software that implemented the algorithm relies on parameters that must be decided through expert supervision, with a trial and error approach. Although this is legitimate, the approach does not provide a clear definition for the characteristic scale of the landscape classification, due to the inherent subjectivity in the selection of the parameters.
Recently, Luo and Liu (2018) combined the "geographical detector method" (Luo et al., 2016) to the geomorphons method (Jasiewicz and Stepinski, 2013), and applied it to landslide susceptibility modelling . This is similar to the work of Alvioli et al. (2016) who used the outcome of multiple statistically-based landslide susceptibility models prepared exploiting SU-based terrain subdivisions with different spatial detail to decide an optimal scale for their terrain subdivision.
An alternative approach to decide the characteristic scale of a terrain classification relies on graph theory (Phillips, 2012;Phillips et al., 2015;Heckmann et al., 2015). Phillips (2012) showed that graph theory is able to describe scale hierarchies, explaining how geomorphic effects are transmitted across scales, despite the lack of evident links between geomorphological phenomena occurring at different scales. Heckmann et al. (2015) suggested that graph theory can be used to explore geomorphic and hydrologic connectivity, landscape evolution models, and various natural hazards.
In the literature, the geographical extent of terrain classifications span a wide range. Examples exist of terrain classifications for areas of a few square kilometres (Dekavalla and Argialas, 2017), of hundreds (Drăguţ and Blaschke, 2006;Zhu et al., 2018), thousands (Schaetzl et al., 2013;Melelli et al., 2017) hundreds of thousands (Guzzetti and Reichenbach, 1994) of square kilometres, for the conterminous United States of America (Luo et al., 2016), and even for the entire globe (Iwahashi and Pike, 2007;Iwahashi et al., 2018;Drăguţ and Eisank, 2012). The difficulty in classifying large and complex landscapes using consistent criteria does not stem from computational limitations. Modern computers can handle regional and even global data sets with spatial resolutions of a few tens of meters, and are able to execute complex algorithms on such large data sets (Drăguţ et al., 2011;Drăguţ and Eisank, 2012;Ehsani and Quiel, 2014;Marchesini et al., 2014;Alvioli et al., 2016). Instead, the main obstacle hampering the classification of complex landscapes or large territories is the selection of the parameters that control the segmentation process i.e., the characteristic scale and the related level of detailof the terrain classification. A related issue is the spatial resolution of the data used for terrain classification (e.g., Mashimbye et al. (2014), Kramm et al. (2017) and Schlögel et al. (2018)).

Methods
Delineation of SUs in a complex landscape adopting the hydrogeomorphological approach requires an a piori selection of the parameters used to extract from a DEM the networks of drainage and divide lines that bound the SUs, which determines the level of detail of the final SUs mosaic. Selection of the parameters is typically heuristic, hampering the implementation of automated procedures for SUs delineation (Erener and Düzgün, 2012;Sharma and Mehta, 2012;Zhao et al., 2012;Mashimbye et al., 2014). In an attempt to solve the problem, we propose to decide the characteristic scale of an SUs subdivision using solely information obtained from a DEM, (i.e., terrain aspect homogeneity constraints) and a minimum size for the SUs.
Working in the Upper Tiber River basin, Central Italy, Alvioli et al. (2016) tested the r.slopeunits software in association with an optimisation algorithm that considers the performances of (i) an SUs division based on terrain aspectusing an objective function adapted from Espindola et al. (2006), and of (ii) a statistically-based landslide susceptibility zonation )using the area under the receiver operating characteristic curve, AUC, as a performance metric.
In this work, we relax the performance criterion based on the landslide susceptibility assessment, and we keep the criterion that relies on terrain aspect segmentation. Application of the original algorithm of Alvioli et al. (2016) has conceptual limitations and computational constraints when applied over very large areas. We propose a new optimisation algorithm to overcome the conceptual and the operational issues, and to publish a national SU map. Optimisation criteria based only on morphometric quantities will expand the potential uses of the final SU map. Identification of SUs with terrain domains corresponding to a correct segmentation of the aspect map matches with the manual delineation of hillslopes, and using a hydrological subdivision at the correct level of detail reflects the geomorphological processes that shape natural landscapes. These requirements are embedded in the use of the two main parameters of the software, described in the following.

Main slope units delineation
The r.slopeunits software (Alvioli et al., 2016) implements an algorithm that defines a mosaic of SUS bounded by drainage and divide networks based on the output of the GRASS GIS r.watershed module (Metz et al., 2011). The module extracts the boundaries of hydrological "half-basins" from terrain information available from a DEM. To single out individual SUs from the half-basins, r.slopeunits assumes that homogeneity (heterogeneity) within (between) SUs is controlled by the variability of terrain aspect in each SU. As a result, parts of a half-basin with substantially different average terrain aspect are considered pertaining to distinct SUs, and the half-basin is split accordingly.
The software implements an iterative procedure ( Fig. 1) beginning with the delineation of a set of a few large half-basins characterised by a large upslope contributing area, T. The initial value of T should be large enough so that the corresponding basins are much larger than the final SUs, and the map does not keep memory of the value of T. The actual value depends on the DEM and cannot be established a priori. The half-basins that meet the internal homogeneity and size criteria (based on a few parameters) are recognised as SUs, and flagged accordingly. The remaining half-basins are sent to the next iteration, where a new set of half-basins is generated using a smaller contributing area, T. The new, and smaller, half-basins that meet the homogeneity and size criteria are recognised as SUs, and flagged accordingly. The procedure is repeated until no half-basins are left, and the entire study area is covered by a mosaic of SUs, of different sizes.
Half-basins are recognised as legitimate SUs based on two input parameters of the software, namely: (i) the SU minimum planimetric area, a [m 2 ], that represents the size below which, at each new delineation of half-basins during in the iterative process, a half-basin is considered a legitimate SU, and (ii) the terrain aspect circular variance, c, that represents the minimum terrain homogeneity for a half-basin to be considered an SU. By definition, 0 ≤ c ≤ 1, and the smaller the c the more uniform the half-basin, in terms of terrain aspect (Alvioli et al., 2016). Since the size of the contributing area T is reduced at each iteration, the size of SUs changes across the study area, adapting to the variable landscape morphology.
Locally, r.slopeunits may produce a few unrealistically small polygons. This is due to inconsistencies between the hydrographic network obtained from any DEM and the real network, as far as small details are concerned. The small areas just represent mistakes, and do not reveal any real topographic property. We can see appearance of small areas as random "noise" underlying the actual topographic "signal". We removed such areas adopting two strategies. A first strategy selects SUs smaller than 50,000 m 2 , and merges them to the adjacent SUs that share the longest boundary. This first strategy was used during the optimisation stage of the procedure. A second strategy selects SUs smaller than 100,000 m 2 , and merges them to the adjacent SU that has the most similar average terrain aspect, where the two SUs share a common boundary longer than half of the perimeter of the SUs to be merged. The second strategy is more accurate than the first one, and was used for the final stage of SUs delineation, after optimisation of the parameters was obtained. In both cases, we removed areas of slightly larger size than the minimum necessary to get rid of mistakes. This is justified by the need of producing useful SUs for practical applications. The number of removed polygons with respect to the total number of SUs is negligible, as well as their total area.

Local optimisation
The iterative procedure described in Section 3.1 can be run for multiple combinations of the (a, c) parameters, and an optimal combination can be selected by maximising an objective function, with respect to (a, c) ( Fig. 1). Alvioli et al. (2016) adapted a metric used to measure the quality of an image segmentation (Espindola et al., 2006) to the segmentation of a terrain aspect map. The segmentation metric imposes that individual slope units in the optimal SU set are as much as possible homogeneous with respect to terrain aspect. Optimization, thus, results in individual slope units to be facing a preferred direction, by selecting the SU sets containing polygons with appropriate shapes and sizes as to adapt to the topography.
The procedure considers circular variance of aspect, instead of terrain aspect itself, to avoid the complexity of the periodicity of angular measures. One obtains the aspect segmentation metric F(a, c) from the local aspect variance V(a, c) and the auto-correlation index I(a, c), as follows: V a;c ¼ ∑ n c n s n ∑ n s n ; ð1Þ where n = 1, …, N ac labels the SUs obtained setting the (a, c) parameters; c n is the circular variance of terrain aspect, and s n is the variance of the SU planimetric area; α n is the average circular variance for the n-th SU, and α is the average circular variance for the entire study area; ω nl is an indicator of spatial proximity, equal to unity where the (n, l) SUs are adjacent, zero otherwise. The quantity V(a, c) in Eq. (1) is the internal terrain aspect variance and I(a, c) in Eq. (2) is the external terrain aspect variance, with minima for SU sets exhibiting well-defined boundaries between adjacent SUs (Espindola et al., 2006). The overall optimisation function F(a, c) that measures the quality of the terrain aspect segmentation is given by the normalised sum: Optimisation of the SUs subdivision consists in maximising F(a, c) in Eq. (3), as a function of (a, c). One can devise a different objective function, for a different purpose. Domènech et al. (2019) for example used the difference between observed and modelled landslide size, where  Alvioli et al. (2016) for the local optimisation of the r.slopeunits (a, c) parameters for areas extending for ≲ a few thousands km 2 . the latter was dependent on the parameters a and c, in their particular model.
In previous works (Alvioli et al., 2016;Schlögel et al., 2018;Bornaetxea et al., 2018), optimisation was performed with F(a, c) in study areas of relatively small sizefrom ≈ 250 km 2 to ≈ 2000 km 2 determined by the extent of the available landslide inventory map. In this work, optimisation is attempted for a much larger area (i.e., the entire Italian territory, N300,000 km 2 ) encompassed by an administrative boundary.

Optimisation in an arbitrarily large area
Extension of the local optimisation algorithm of Section 3.2 to the entire Italian territory into a mosaic of SUs poses conceptual and computational challenges. First, the study area is delimited by an administrative boundary, and not by morphological or hydrological limits shaped by physical processes. Hence, there is no reason to infer that a single combination of the (a, c) parameters will be adequate for the entire study area. Second, the geographical extent of the study area is very large. Using the EU-DEM (https://www.eea.europa.eu) digital elevation model, at 25 m × 25 m resolution, it consists of about 5.30 × 10 6 valid grid cells, corresponding to ≈ 330,000 km 2 . Note that the area covered by the selected EU-DEM grid cells is larger than the geographical area of Italy, because hydrological calculations require cells outside the Italian administrative boundary.
Due to the vector operations required to find neighbouring SUs, calculation of Eq. (2) is very slow. A single run of r.slopeunits on one of the 439 basins require about half an hour, on average; evaluation of the F(a, c) metric for a single SU map requires up to one hour, in each basin. This must be multiplied for 49 different combinations of the parameters a and c, in first place. Secondly, we applied the optimization procedure not only on individual basins, but on a large number of spatial domains of increasing size (cf. Fig. 2). Slope unit sets, instead, were not computed again, but simply aggregated as needed. The overall running time of the procedure amounted to about three months of non-exclusive use of a machine equipped with 64 computing cores and 320 GB RAM memory. The vast majority of the required computer processes were executed automatically, by shell scripts, within GRASS GIS.
To make the problem computationally manageable, maintaining it geomorphologically meaningful, we devised the solution illustrated in Fig. 2. We first prepared a subdivision of Italy into 439 main hydrological basins (black polygons in Fig. 3), ranging in size between 50 km 2 and 4300 km 2 (mean = 741 km 2 , standard deviation = 570 km 2 ), a size comparable to the study area of Alvioli et al. (2016) (i.e., ≈2000 km 2 ), which we maintain is a rough threshold under which the optimisation procedure becomes unreliable. We established this in a test that showed that in much smaller areas the objective function F(a, c), Eq. (3), becomes erratic and does not necessarily have a single maximum.
To simplify the process, we removed the large plains, corresponding to 1.17 × 10 6 grid cells, 73,125 km 2 , shown in white in Fig. 3. We further excluded from the analysis all the islands, except for Sicily and Sardinia. The many other small islands would each represent an isolated, main hydrological basin for which separate parameter optimisation would be required. Since the size of the individual islands is small (b240 km 2 ), their optimisation would offer little insight on the performance of the method, or the results obtained.
Optimisation of the (a, c) parameters was performed separately in domains of increasing size, keeping the hydrological basins shown in Fig. 3 as the main analysis and optimisation domains. To obtain optimal values for the parameters in the j-th basin, we first maximised F(a, c) in the j-th hydrological basin itself. This is represented by the "Step 1" panel in Fig. 2. Next, we searched the maximum of F(a, c) within the areas represented by the geographic union of basins j and k, where k is any of the basins adjacent to j, for every possible k ("Step 2" in Fig. 2). Next, we considered all the possible triplets of adjacent basins including j ("Step 3" in Fig. 2). We repeated the procedure for up to three (K = 3) adjacent basins for all the 439 main hydrological basins shown in Fig. 3. In addition, for a subset of 30 basins located in the Central Apennines (and contained in a topographic unit, see Fig. 12), we extended the analysis to K = 5 adjacent basins, to check the effect of truncation at K = 3. Results revealed that extending the analysis beyond K = 3 adjacent basins did not result in significant differences in the optimised (a, c) parameters, and hence in the SUs terrain subdivision.
To consider that the level of detail of the SUs subdivision should vary depending on the local geomorphological setting, and that the geomorphological setting may change significantly across topographical boundaries, we further constrained the optimisation process within the boundaries of main topographic units. For the purpose, we used the topographic subdivision of Italy proposed by Guzzetti and Reichenbach (1994) who defined eight topographic provinces and 29 sections (red lines in Fig. 4), combining an unsupervised classification of morphometric variables with the visual interpretation of thematic information. This restriction is justified by the fact that even if the topographic provinces are large and contain several of the main hydrological basins shown in Fig. 3, they were defined to encompass regions with similar topographic, morphometric and geomorphological characteristics.
We computed a weighted average optimisation of the optimal parameters obtained at each step, using larger weights for the low iteration steps, and proportionally smaller weights for high iteration steps, as follows: where and K labels different steps, c i is the optimal c found in the i-th combination of n hydrological basins, and s i is the corresponding total planimetric area. For a given n, c n ð Þ j is the weighted average of all the optimal c i values calculated in the spatial domain defined by the i-combination of the j-th basin with the adjacent basins (combinations of two basins for K = 2, three basins for K = 3, …, five basins for K = 5).  Guzzetti and Reichenbach (1994), with red lines, superimposed to a shaded relief calculated from the EU-DEM used in this work. Topographic sections are further grouped in eight topographic sections, as suggested by the names of the 29 sections; names are shown in Table 2. Map is in EPSG:4326, datum in EPSG:3035.
Adopting this procedure for the j-th main hydrological basin, j = 1, …, 439, we obtained a sequence (c opt (1,j) , c opt (2,j) , …, c opt (5,j) ) of values optimised in domains of increasing size, centred in the j-th basin. The same procedure was applied to the parameter a, that controls the minimum size of an SU, to obtain the sequence (a opt (1,j) , a opt (2,j) , …, a opt (5,j) ). From Eq. (4), we observe that the value of c opt K, j for the j-th basin depends mainly on the value of c calculated for the basin j itself (c 1; j opt ¼ c j ), and progressively less from the values of c opt K, j calculated for K = 2, K = 3, …, K = 5.
Next, we investigated whether the two sequences, shown in Fig. 5, admitted asymptotes. Considering the circular variance sequence in the j-th basin, we hypothesised that it admitted an asymptote, c asy (j) .
We then calculated the positive differences between estimates of c opt at two subsequent steps of our sequence of optimal points in the j-th basin. For a decreasing series: while for an increasing series, to have positive differences we calculated instead: so that we can write: where the plus (minus) sign is valid for a strictly decreasing (increasing) series. Taking the limit of Eq. (8) for large K we have: In the last equality we have written d i as Ax i . Now, if |x | b 1, we have a geometric series whose sum is known analytically: Thus, fitting the first differences (up to, at most, i = 5 in our case) of Eqs. (6), (7) to Ax i as a function of i, we obtain the A and x parameters to replace in the following equation for the asymptote, where the plus (minus) sign holds for a strictly decreasing (increasing) series. The procedure presented in Eqs. (6)-(11) can be applied to find asymptotes for the circular variance of the aspect, c, and the minimum SU planimetric area, a. Existence of asymptotes for the two series, c asy (j) and a asy (j) , is not guaranteed for every hydrological basin considered in the analysis (Fig. 3), or for both parameters a and c. We accepted positive values for the asymptotes, providing the optimal value a opt or c opt for the corresponding main hydrological basin. Where an asymptote was not found, the search procedure failed, and we used the values of (a opt K=1,j , c opt K=1,j ) corresponding to the "local" minimum, in the j-th basin, of the F(a, c) function given by Eq. (3). In our experiment, this occurred for~150 basins,~34% of the 439 main hydrological basins shown in Fig. 3.

A national slope units subdivision of Italy
The obtained terrain subdivision consists of 331,926 SUs in Italy, ranging in size from 1874 m 2 to 15,493,140 m 2 , mean = 229,655 m 2 , standard deviation = 821,887 m 2 . Due to the complexity of the SUs subdivision, we cannot show the entire map for Italy. For illustrative purposes, in Fig. 6 (a) we show the SUs for the Nera River basin, Central Italy, and in Fig. 6 (b) the local search for a maximum of F(a, c), as a function of (a, c), for the same catchment, performed adopting two strategies: (i) a pre-defined grid of (a, c) pairs (purple surface), and (ii) a maximum-bracketing algorithm (blue line), adopted when the series of parameter values did not show asymptotic convergence. The purple and the blue dots show the maxima obtained by the two strategies. Fig. 7 shows examples of the SUs subdivision for four different and characteristics geomorphological settings in Italy, including basins (i) in the Northeastern Alps, near Bolzano, (ii) in central Italy, along the Adriatic Sea near Ancona, (iii) in Southern Italy, close to the Ionian Sea near Riace, and (iv) in Sicily, along the slopes of the Etna volcano. Stars in Fig. 3 correspond to the locations of the basins in Italy. Panels in the right column of Fig. 7 show the cumulative frequency-size distributions of the corresponding SUs areas. Basins close to the Etna volcano and the Bolzano basin are characterised by SUs larger than in the other  Table 2).
basins, but lithology and morphology are different. In the Bolzano basin, the relief is high and lithology is represented mainly by schists and carbonate rocks. The Mount Etna basin is partly dominated by volcanic morphology due to the presence of basalt and lava flows, and partly by mountains where turbidites and chaotic rocks prevail. In the Ancona basin, marls and overall turbidites outcrop in the mountains, where slope units are larger in size. The hills located closer to the Adriatic sea have smaller slope units characterised by unconsolidated sediments, typically sand, silt, and clay. The hilly part of the Riace basin is dominated by intrusive rocks, whereas, in the coastal part, consolidated and unconsolidated sediments crop out with turbidites, schists and chaotic rocks.

Performance of the slope units subdivision of Italy
Due to the complexity of the Italian landscape (Fig. 3), the large number of SUs (331,926) and their size variability (Fig. 7), assessing the c orrespondence of the SUs terrain subdivision with the geomorphological and hydrological reality i.e., the performance of the SUs subdivision, is not a trivial task. We attempt to perform the task in steps of increasing complexity, as follows: • First, we study the relationship between SUs size, terrain elevation and lithological types. • Next, we examine the size of the SUs within the topographic provinces and sections proposed by Guzzetti and Reichenbach (1994), and within the topographic classes proposed by Drăguţ and Eisank (2012). • Lastly, we compare this topographic subdivision of Italy with a new terrain subdivision obtained by unsupervised clustering of SUs morphometric variables.
The next three sections describe separately the results corresponding to the points above.  Fig. 3. Boxes on the right show the non-exceedance cumulative probability for landslide sizes. Maps are in EPSG:4326, datum in EPSG:3035.

Slope units, drainage density and terrain elevation
The size distribution of the SUs is strictly related to drainage density, a hydrological and geomorphological characteristic descriptive of many landscapes (Carlston, 1963). Drainage density of a given basin is defined as the sum of the lengths of all the streams in a basin, divided by the total area of the basin.
Existing approaches for the calculation of drainage density from a DEM rely on an arbitrary choice of a threshold value to initialise the calculation of the stream network. In this study, to give an abjective measure, we adopted the number of slope units per unit area (in the following, SU/km 2 ) as a proxy for drainage density. We maintain that the size of SUs produced by r.slopeunits changes across the study area, adapting to the local landscape morphology, and making SU/km 2 a robust objective approximation of the actual drainage density in different regions of the study area.
In practical terms, we used the geographical location of the centroids of SU polygons to calculate the values of SU/km 2 , providing an average drainage density for each area. The first assessment of drainage density distributions consisted in the comparison of the values of SU/km 2 in different classes of terrain elevation, obtained from the EU-DEM, all over Italy. The comparison revealed that the average drainage density decreases linearly with increasing terrain elevation, Fig. 8 (a). Results indicate that terrain elevation is linked to the average size of the SUs, which are (generally speaking) larger in the mountains than at low Fig. 8. The values of SU/km 2 , a proxy for drainage density, within ten percentiles of elevation, in (a), and within eight classes from the topographic classification based on elevation by Drăguţ and Eisank (2012), in (b). The cartoons in the insets show the geographical location of the different percentiles and of the classes, respectively. Topographic classes in (b) are as follows: 1, tablelands; 2, flat plains; 3, irregular plains; 4, smooth low hills; 5, rough low hills; 6, high hills; 7, low mountains; 8, high mountains. elevations e.g., near the large plains. For a length scale L equivalent to the average size A SU of an SU, L = ffiffiffiffiffiffiffi ffi A SU q = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 299; 655 m 2 p ≈ 696 m, at lower elevations relative relief is small, slopes are typically short and small, and drainage density is higher than at higher elevations where relative relief is high, slopes tend to be long and large, and drainage density is typically low.
The global topographic map of Drăguţ and Eisank (2012) is an independent classification which we can use to assess the performance of the SU map. They used Shuttle Radar Topography Mission elevation datathus, different from EU-DEM of our analysisand performed a terrain classification using elevation itself, and standard deviation of elevation. We analysed the relevant portion of the global map in a similar way as in the analysis of percentiles of elevation described in this section. Fig. 8 (b) shows the corresponding results, along with the list of topographic classes identified by Drăguţ and Eisank (2012); we disreagarded other terrain attributes in the map, for simplicity. Results are essentially consistent with the ones presented for percentiles of elevation, Fig. 8 (a).

Slope units and lithology
Drainage density is known to be controlled by lithology (Melton, 1957;Sangireddy et al., 2016). We calculated the values of SU/km 2 in different lithological domains. For the purpose, we used a lithological domain map of Italy, in 9 classes, modified from Bucci et al. (2018) and  at 1:100,000 scale, shown in Fig. 9. The analysis revealed that unconsolidated clastic rocks, even including moraine, glacial drift and landslide deposits, exhibit the largest drainage density, corresponding to a finer SUs granularity characterised by small and very small SUs, (≈ 2/km 2 Table 1), whereas foliated (e.g., schist, shale) and non-foliated (granular, e.g., gneiss, quartzite, marble) metamorphic rocks, have the lowest drainage density and a corresponding coarser SUs granularity characterised by many large SUs, (≈1.16/km 2 Table 1).
Pyroclastic rocks (mainly located in the volcanic provinces of Tuscany and Latium), consolidated clastic rocks (e.g. conglomerates) and chaotic rocks/ melange outcrop.

Slope units and topographic divisions
We checked whether the SUs partition, and particularly its spatial variation, was able to capture geomorphological differences outlined by the topographic subdivision proposed by Guzzetti and Reichenbach (1994). Fig. 10 shows the distribution of the SUs size in the eight topographic provinces and the 29 sections of Guzzetti and Reichenbach (1994), and Table 2 lists statistics for the SUs in the different areas.
Inspection of the figure and the table does not reveal large differences in the average size, or the range of the SUs area among the different topographic provinces and sections. This was expected, because a single metric (i.e., drainage density) cannot capture entirely the topographic and geomorphological complexity, and the geological variability, of the Italian landscape.
However, minor differences exist in the topographic provinces and sections, which are worth discussing. At the level of the eight topographic provinces of Guzzetti and Reichenbach (1994) (Fig. 10) we observe that Alpine-Apennines Transition Zone (3.), Adriatic Borderland (6.), Sicily (7.) and Sardinia (8.) show very similar values of mean Fig. 10. Box-plots show distribution of SUs area, [km 2 ], for the eight topographic provinces and 29 topographic sections of Guzzetti and Reichenbach (1994) (see Fig. 3). Top horizontal axis shows number of SUs in each topographic province, (a), and in each section, (b). Names and properties of the topographic provinces (a) and sections (b) are listed in Table 2. average SUs area and range. As expected, the drainage density (column E in Table 2) is also similar among the four topographic provinces and close to 1.6 SU/km 2 . Only the Alpine Mountain System (1.) and the Apennine Mountain System (4.) have higher values of the mean and the range of the SU sizes. The North Italian Plain (2.) and Tyrrhenian Borderland (5.) are characterised by smaller slope units. In the Alps topographic province (1.1, 1.2 and 1.3 in Fig. 10), SUs have the largest median area and the largest range in area, followed by the Apennines topographic province (4. in Fig. 10), confirming a dependence of the SUs size on terrain elevation and relative relief. Accordingly, drainage density is low (column E in Table 2). This is justified by the morphological and geological complexity, and the large extent of the two topographic provinces (Table 2). Conversely, the North Italian Plain (2., actually the hills close to the plains), the Tosco-Laziale section (5.2), the Campidano Plain (8.3) and the Gargano Upland (6.3 in Fig. 10) (Fig. 10). In the Apennines range, the Aspromonte Apennines section (4.7) has the lowest average and range of SUs area, indicating a strong homogeneity of the terrain.
The largest average SUs area is found in the Western Alps (1.1), that contains the Mont Blanc range (4808 m a.s.l.), the tallest mountain in the Alps, followed by the Molise (4.3) and the Molise-Lucanian Hills (4.4) Apennines sections. The three sections also exhibit the largest ranges in SUs area, a measure of the topographic heterogeneity of the sections.
Of the other topographic provinces, the Gennargentu Highland (8.2), the Sicilian Hills (7.1 & 7.2), the Central Apennines Slope (6.1) and Ligurian Upland (3.2) have medium to high values of the mean area and the range of the SUs area (Fig. 10). In Sicily, the Etna volcano section (7.4), with a low mean SUs area and a low range of variability of the SUs area, has a distinct signature from the Sicilian Apennines section (4.8), and the Sicilian Hills (7.1 & 7.2), with the former showing a smaller mean SUs area and larger drainage density than the latter. In Sardinia, the smallest average SUs are found along the foothills surrounding the Campidano Plain (8.3), and the largest average SUs are found in the Gennargentu Highland (8.2). The two topographic sections also exhibit the lowest (8.3) and the largest (8.2) variability of the SUs area, confirming that around large plains the SUs tend to be small and with high drainage density.

Slope units and terrain zonations
To further evaluate the performance of the SUs subdivision, we performed a final analysis comparing the topographic division of Italy proposed by Guzzetti and Reichenbach (1994) (Fig. 3) with a new terrain zonation of Italy prepared for this study through multivariate clustering based only on quantities calculated from slope units.
We devised a clustering procedure as follows. We calculated for each of the 439 basic hydrological basins (our starting point for nation-wide SU delineation, cfr. Fig. 3) the distributions of (i) the average aspect (actually described by circular variance) over the cells within each slope unit polygon contained in the basin, and (ii) the size distribution of SUs in the basin. From both distributions, we calculated the 10th, 20th, …, 90th percentiles of the distribution, obtaining a list of 9 + 9 = 18 values, for each of the 439 basins. The resulting 439 × 18 matrix was used as the only input for a clustering procedure. In order to let the automatic classification to tell something about how good was the SU delineation over the whole of Italy, we did not include as an input any other information, neither morphological nor thematic.
The clustering results were compared with the following quantities: (i) the frequency distribution of terrain elevation in the SUs, in five classes, corresponding to the 20th, 40th, 60th, 80th and 100th percentiles of the distribution of the elevation, (ii) the frequency distribution of terrain slope in the SUs, in five classes, corresponding to the 20th, 40th, 60th, 80th and 100th percentiles of the distribution of terrain slope, and (iii) lithological information obtained from a lithological domains map of Italy at 1:100,000 scale, in 9 classes, adapted from Bucci et al. (2018) and .
For the analysis, we adopted unsupervised K-means clustering (Hartigan and Wong, 1979), and specifically its implementation in the kmeans() package available for the R free software environment for statistical computing (http://www.r-project.org). A literature search revealed that no unique or preferred method exists to determine the number of clusters for K-means unsupervised clustering. Existing methods include the "elbow method" (Thorndike, 1953), the Akaike information criterion, AIC (Akaike, 1974), and the Bayesian information criterion, BIC (Schwarz, 1978;Ramsey et al., 2008), all of which evaluate an objective function to establish the optimal number of clusters. We experimented with different numbers of clusters, from 2 to 50, and with different methods for establishing the optimal number of clusters, obtaining the results shown in Fig. 11. Given that the three metrics gave three different results, we selected seven clusters, which is at the edge of the minimum shown by the BIC method. We maintain this is a good compromise between the complexity of the Italian landscape and the need to have a reduced number of clusters. Preliminary calculations with larger numbers of clusters showed that interpretation of the Table 2 Characteristics of SUs in the different topographic provinces and sections of Guzzetti and Reichenbach (1994). A: Topographic province (rows in bold) or section; B: Extent [km 2 ]; C: Percent; D: SU area [km 2 ]; E: SU/km 2 , a proxy for drainage density. See Fig. 3 for map of the geographical distribution of the topographic provinces and sections. results would be hardly understandable in simple, and useful, terms. We show in the following that already using seven clusters it is difficult to indentify striking differences between their characteristics. Here, we further note that our optimal number of clusters (seven) is similar to the number of topographic provinces (eight) in Italy proposed by Guzzetti and Reichenbach (1994).
To size the differences between the seven clusters we compared the generalised Euclidean distances between the cluster centroids, with smaller (larger) differences measuring more similar (different) clusters.
Inspection of Figs. 12, 13, 14 and 15, and of Table 3 and Table 4 reveals differences and similarities among the seven clusters. Clusters 1 and 7 match topography close to tabular. They have a similar topographic settingconfirming their similarity, Fig. 12 (b) with the proportion of low elevation terrain, and of flat and gentle terrain, significantly higher than the proportion of high and very high elevation terrain, and of steep and very steep terrain. Unconsolidated clastic rocks (L2) and (mostly massive) carbonate rocks (L5) are the main rock types in the two clusters, with carbonate rocks, locally carved by karst processes, particularly abundant in cluster 1 (Table 3). Cluster 1 is confined to the hills bounding large plains and to areas with a tabular morphology (e.g., in the Carso area, in the Apulia Lowland), and cluster 7 is practically absent in mountain areas. The lithological distribution of cluster 7 and of the entire Italian territory (last row in Fig. 12) is very similar, apart from the alluvial deposits, which, as a matter of fact, are underrepresented in the area covered by SUs. This was in some way expected since cluster 7 largely spreads all over the Italian territory. Clusters 2, 3 and 5 are typical of mountainous areas. They exhibit the opposite topographic setting, with the proportion of high and very high elevation terrain, and of steep and very steep terrain significantly higher than the proportion of low and very low elevation terrain, and of gentle and flat terrain. This is particularly evident in clusters 2 and 3 that have the highest abundance of high elevation and steep terrain of all the seven clusters. Metamorphic rocks (L9), followed by carbonate rocks (L5), unconsolidated clastic rocks (L2), and turbidites (L4), are the main rock types in clusters 2 and 3, whereas turbidites (L4) and subordinately, unconsolidated clastic rocks, metamorphic (L9) rocks and carbonate (L5) predominate in cluster 5. Clusters 4 and 6 exhibit an In the map, the 439 main hydrological basins used in the study (see Fig. 3) are coloured based on the results of unsupervised k-means clustering of SUs terrain variables and lithological information. Black lines show boundaries of eight topographic provinces of Guzzetti and Reichenbach (1994). A few unclassified basins are shown in white along with large plains, excluded from the cluster analysis. (b) normalised distance matrix selected automatically by the k-means algorithm. (c) extension of the different clusters, expressed as a percentage of the total area. Map is in EPSG:4326, datum in EPSG:3035. Fig. 11. Different metrics, denoted generically by f(N), for the optimisation of the number of clusters N, using the k-means algorithm. Red: within-cluster sum of squares, objective function for the elbow method; p ij and c ij represent the data points and the cluster centroids, respectively, while the i, j and k label the variables, the points in the cluster and the cluster itself, respectively. Green: the Akaike information criterion. Blue: the Bayesian information criterion.
intermediate topographic setting, with terrain (almost) equally abundant in all the considered elevation and slope classes, close to the percentile distribution of slope and elevation of the whole Italian territory. Sedimentary rocks are the most common rock type in the two clusters; mostly unconsolidated clastic rocks (L2) in cluster 4, and turbidites (L4) in cluster 6.
We further note that several of the hydrological basins draining the northeastern and eastern slopes of the Apennines range into the Adriatic Sea (Fig. 12) were classified alternatively as pertaining to clusters 3, 4, 5 or 6, despite the fact that the morphology of the catchments appears similar (Fig. 3). This is because these catchmentswhich are typically elongated in shapespan across two topographic provinces e.g., the Apennine Mountain System (4) and the Adriatic Borderland (6), and locally across two or more topographic sections characterised by different topographic settings (Fig. 7). As a result, the unsupervised clustering algorithm consistently attributed the single catchments to a cluster depending on the prevalent proportion of the clusters, and the corresponding topographic and lithological characteristics, in each catchment. This resulted in an apparently odd classification of the hydrological basins.
A few topographic sections contain significant proportions of only two or three clusters (excluding the large plains) indicating the overall topographic homogeneity of the sections. Sections with significant proportions (N5%) of two clusters include (i) the Gargano Upland (6.3) that Fig. 13. Histograms of elevation, slope and lithology classes in the seven clusters singled out by unsupervised k-means clustering (Fig. 12). See text for explanation. Blue and red bars show five terrain slope and elevation classes, respectively, both in order of increasing values; colour-filled bars show lithology classes, L1-L9, shown in Fig. 9. For each cluster, we specify the number of basins and the extent of the area covered by the basins. Fig. 15. Distributions of SUs areas, (a), and of SUs per km 2 (a proxy for drainage density), (b), within the clusters singled out by the k-means algorithm (cfr. Fig. 12). Labels on the top axis correspond to the number of SUs in each cluster, in (a), and to the average SU/km 2 , in (b), in each cluster. The total number of SUs shown in this Fig. (394,305) is larger than the total number of SUs quoted throughout the paper (331,926), due to clusters extending outside the Italian administrative borders (cfr. Figs. 3 and 12).  Guzzetti and Reichenbach (1994), the percentage of terrain assigned to each of the seven clusters singled out by k-means unsupervised clustering (Fig. 12). Code names used for the topographic sections are listed in Table 2. The "Plains" class (white) corresponds to plain areas in Figs. 3, 12, which were excluded from SUs delineation. See text for explanation.
contains a large proportion of cluster 1, and some cluster 4, (ii) the Etna section (7.4) with a large coverage of cluster 6, and some cluster 4, (iii) the Sardinian Hills (8.1) dominated by cluster 7 with some cluster 6, (iv) the foothills surrounding the Campidano Plain (8.3) that have equal parts of clusters 6 and 7. Sections with significant proportions (N5%) of thee clusters include (i) the Ligurian Upland (3.2), (ii) the Aspromonte (4.7), and (iii) the Iglesiente Hills (8.4). A few sections contain significant proportions (N5%) of six or seven different clusters (excluding the large plains) revealing the topographic complexity and variability of the sections. The Central-Eastern Alps section (1.2) contains all seven clusters, a result of the large extent and complexity of the province, whereas sections with significant proportions (N5%) of six clusters include (i) the Western Alps (1.1), (ii) parts of the Alpine foothills (2.3a), (iii) the Molise Apennines (4.3), (iv) the Lucanian Apennines (4.6), and (v) the Central Italian Hills (5.1).

General results
The purpose of our work was to devise and test a parameter-free procedure for the automatic subdivision of a large and complex landscape into a mosaic of SUs using solely terrain information obtained from a DEM. In previous works, for the optimisation of the parameters that control the automatic partitioning of a territory using the r. slopeunits software, Alvioli et al. (2016), Bornaetxea et al. (2018) and Schlögel et al. (2018) adopted an approach that optimised the parameters maximising jointly a terrain segmentation metric (Espindola et al., 2006) that considered the homogeneity of terrain aspect inside the SUs, and a measure of landslide susceptibility performance (AUC, Fawcett (2006)). Albeit successful, the approach had operational and conceptual limitations. The operational limitation was that it required the availability of landslide and thematic information used to prepare a sufficiently large set of landslide susceptibility models using different SUs terrain subdivisions. The conceptual limitation was that it provided an SUs subdivision optimised for landslide susceptibility assessment. However, the terrain subdivision may not have been adequate for other applications. Even for landslide susceptibility assessment, the SUs subdivision was optimised for a specific landslide type e.g., small shallow soil slides, deep-seated slides, very large deep seated slides, and multiple SU delineations would be necessary to account for all landslide types in an area.
The approach proposed and tested in this work proved capable of overcoming the operational and conceptual limitations. The procedure exploited solely terrain information extracted from the DEM, without any additional thematic information or expert judgement. The approach proved to perform well in different topographic and geomorphological settings, adapting to the local terrain conditions (Figs. 6, 7). We maintain this is a significant advantage over previous works (Alvioli et al., 2016;Bornaetxea et al., 2018;Luo and Liu, 2018;Schlögel et al., 2018).
We stress that the ability of the procedure to adapt to the local terrain conditions, and to produce a single terrain subdivision of Italy into a single set of SUs, has significant advantages for geoenvironmental modelling. Use of a single terrain subdivision facilitates the comparison and the exportability of the modelling results for a single hazard (e.g., for landslide susceptibility modelling), and for multihazards and multi-risk assessments (Gill and Malamud, 2014).

Specific results for Italy
Literature analysis reveals that the classification of the Italian landscape in terms of geomorphological and topographic characteristics is still poorly investigated. After Guzzetti and Reichenbach (1994), recently Smiraglia et al. (2013) published a "Land units map of Italy" with a subdivision of Italy in 149 land facets based on macrobioclimatic, lithological and geomorphological variables. To our knowledge, a detailed slope units subdivision of the entire Italian territory (more than 300,000 km 2 ) was not available before our work. Moreover, the delineation of 439 main hydrological basins and their classification in terms of topography is novel. The SUs subdivision is well correlated with terrain elevation (Fig. 8). This demonstrates that the adaptive algorithm used for the slope units delineation was able to capture the local terrain morphology (sec. 4.2.1). SUs tend to be larger when terrain elevation is high, and this obeys the general expectation that a mountain slope should be larger in size than a hill slope. Analysis of an independent terrain classification, available globally and performed with an objectoriented technique, confirmed the same tendency for the correlation of slope unit sizes with terrain types at different elevation.
SUs were also compared with a lithological classification of Italy. Overall, we observe that more hardened, foliated or fractured rocks are characterised by lower drainage densities (measured here by SU/ km 2 ) and larger SUs, whereas soft and poorly hardened rocks, including e.g., clayish and sandy, scarcely consolidated sediments, coarsely bedded and granular rocks, exhibit higher drainage density and smaller SUs (sec. 4.2.2). The findings suggest that the lithological characteristics affect topography, which in turn influence SUs delineation.
The 439 basins delineated in Italy were classified using 7 clusters defined from the SUs subdivisions (section 4.2.4). Clusters were also described in terms of the distribution of elevation and terrain slope, and based on the presence of different rock types. At least three different regions (mountainous, hilly, and nearly flat) can be singled out based on Table 3 Quantitative characteristics of SUs in the different clusters singled out by k-means clustering, shown in Fig. 12. See Table 3 Table 4 Descriptive characteristics of SUs in the different clusters singled out by k-means clustering, shown in Fig. 12. See Table 3 for additional properties; Table 2 and Fig. 4 for topographic sections; Fig. 9  elevation and terrain slope in the different clusters. These regions can be further distinguished based on the presence of different rock types in the clusters (Fig. 13). The map in Fig. 12 shows that the classification of the 439 basins exhibits a pattern where some of the morphological and lithological characteristics of Italy emerge. For example, the Alpine and Apennines mountain ranges, the carbonatic areas (e.g., Apulia and Gargano Upland, Iblei Plateau), the volcanic structures (Etna mountain), the hilly areas (Sicilian and Sardinian Hills) are classified as distinct clusters. This is further evidence of the effectiveness of the SUs subdivision, independent on any data beyond elevation. The subdivision is able not only to capture different morphology at the local scale (down to a single hillslope, i.e., a slope unit), but also at basin scale and intra-basin scale, as demonstrated by the clustering results.

Potential applications
We expect the obtained terrain subdivision of Italy in a dense mosaic of SUs to have several potential applications. For landslide studies, the primary application will be the production of statistically-based landslide susceptibility modelling , and related terrain zonations (Alvioli et al., 2016;Ba et al., 2018;Bornaetxea et al., 2018;Luo and Liu, 2018), at different geographical scales, from the local (catchment) to the synoptic (national) scales. Adoption of the same terrain subdivision in different areas will facilitate model comparison and the exportability of the model results. Additional applications for landslide studies include (i) improvement of the quality and completeness of event landslide inventory maps using satellite imagery , (ii) aggregation of pixel-based results of physically-based landslide simulations Domènech et al., 2019), (iii) enhanced evaluation of the quality and completeness of landslide maps and of the quality of landslide susceptibility models , (iv) extraction of bedding and structural information from stereoscopic aerial or satellite images, and (v) assessment of the influence of morpho-structural settings on landslide types, abundance and pattern Santangelo et al., 2015) and (vi) assessment and modelling of earthquake induced landslides (Tanyaş et al., 2019a(Tanyaş et al., , 2019b where directivity of slope dynamics response to seismic shaking is important (Del Gaudio and Wasowski, 2007).
The terrain subdivision of Italy into SUs can also be used e.g., for hydrological and erosion modelling, and for geo-environmental, ecological, forestry, agriculture and land use/land cover studies that require the identification of homogeneous terrain domainsor mapping units facing distinct directions (Smiraglia et al., 2013;Belyanin, 2017;Liu et al., 2018;Tracz et al., 2019;Hu et al., 2018). We further expect the subdivision to be useful for terrain visibility and illumination studies (Minelli et al., 2014;Notti et al., 2014). Of particular interest is the potential use of the SUs subdivision for enhanced landscape classification (Evans, 2003;Minár and Evans, 2008;Evans, 2012), combining the morphometric information captured by the hydrologically-based terrain subdivision with additional and complementary morphological, lithological and geological information, including e.g., the topographic subdivision of Italy proposed by Guzzetti and Reichenbach (1994).

Conclusions
We developed a method for the automatic delineation of slope units (SUs), geomorphological terrain units delimited by drainage and divide lines. The method relies on the r.slopeunits software, a GRASS GIS module introduced by Alvioli et al. (2016) to extract drainage and divide lines from a DEM in an adaptive fashion, and is capable of handling study areas of arbitrarily large sizes.
We applied the approach on the whole of Italy, covering more that 300,000 km 2 . The method is parameter-free, since the input parameters of the underlying software r.slopeunits used for local SU delineation can be optimised in an objective way. The method was designed to either reach for convergence towards optimal parameter values, in optimisation domains of increasing size, or to apply a local optimisation algorithm, wherever convergence did not emerge. We eventually produced an SU map for the whole of Italy, which we distribute in 20 (overlapping) vector layers, corresponding to the 20 Italian administrative boundaries. A sample of the results is shown in Figs. 6 and 7 for a single basin, for illustration purposes.
To investigate the robustness of the optimal SU map, and the performance of our parameter-free SU delineation algorithm, we have investigated the relation between seemingly unrelated quantities: the distributions of the SU sizes and the average aspect within them, on one side, and the slope, elevation and lithological distributions all over Italy, on the other side. We used a k-means clustering algorithm to perform such a comparison. Given that the variables used for classification (percentiles of distributions of SU sizes and aspect variability), were totally unrelated to the variables considered to validate the clustering results, we believe that classification returned an interesting output. It shows that SU have size, shape and aspect variability which is different in inherently different geographic areas of Italy. In our opinion, this is an additional indication that SU as delineated by the r.slopeunits software and whose input parameters were optimised with the algorithm presented in this work, provide a meaningful description of the terrain and have relevant information content.

Acknowledgments
This work was partially supported by RFI Rete Ferroviaria Italiana gruppo Ferrovie dello Stato Italiane.
Appendix A. Obtaining the software and map of the terrain subdivision for Italy Interested readers can download the SUs subdivision of Italy from the following URL: http://geomorphology.irpi.cnr.it/tools/slope-units. The terrain subdivision is provided in vector format for the whole of Italy (331,935 SUs covering 231,762 km 2 ) split in the 20 Italian administrative Regions. Note that SU polygons crossing the border of two Regions were included in both the corresponding vector layers, and the same goes for SU polygons falling partially outside the Italian national borders. Research scientists and practitioners interested in testing the proposed method in other geographical areas can obtain the software from the same URL as above, or can email the corresponding Author. The terrain subdivision and the r.slopeunits software are provided "as is", without warranty of any kind, either expressed or implied. Use of the software and/or the terrain subdivision does not imply the endorsement of the authors or their Institution, the Italian Consiglio Nazionale delle Ricerche.