Mapping Geodiversity at a National Scale: the Case Study of Italy

In order to assess the geodiversity of the Italian Peninsula, which covers approximately 300.000 km2, a semi-quantitative method based on the use of grids recording several indicators and indices was developed. The variety of geological, geomorphological, and pedological elements, characterizing the Italian territory, has been assessed with a two-step procedure. Firstly, the variety algorithm has been applied using grid cells with variable size, related to the spatial resolution of the input data, then the resulting variety values were averaged with a fixed cell size functional to the extent of the study area and the output scale of the geodiversity map. This procedure made it possible to preserve the spatial resolution of the input data (Digital Terrain Model, lithological and soil maps) providing as output a geodiversity map that faithfully reproduces the features of the Italian territory. In case of discrete data (rivers, lakes, glaciers, etc.), a procedure that assigns to each cell the maximum area or length values out of all its elements has been implemented. It made possible to preserve the hydrological elements that shape the landscape (e.g., the longest rivers, largest lakes, etc.) and represent important freshwater resource. An overview of the geographical distribution of geodiversity classes over the whole Italian territory has been elaborated. The resulting geodiversity map is a valuable tool for environmental planning, in particular for the identification of areas to be preserved, for the proper management of geo-resources and natural services.


Introduction
Geodiversity is commonly defined as the variety of geological (rock, mineral, fossil), geomorphological (landform, processes), soil and hydrological (springs, swamps, lakes, rivers) features (Gray 2004(Gray , 2008;;Serrano and Ruiz-Flaño 2007) typifying a territory.It represents the basis upon which all creatures, from plants to humans, live and interact (Gordon and Barron 2013;Matthews 2014;Lawler et al. 2015) and, as the expression of the variability of the Earth's physical environment, together with climate, has a crucial influence on biodiversity (Hjort et al. 2015;Crisp et al. 2021;Comer et al. 2015;Zarnetske et al. 2019).
Geodiversity provides support for geo-conservation actions (Dias et al. 2021), allows the identification of geo-diverse areas that need priority protection (Brilha et al. 2018) and helps in understanding the history of the Earth and the evolution of life.Furthermore, geodiversity actively contributes to several ecosystem services and abiotic ecosystem services, also named geosystem services (Gray 2012;Fox et al. 2020;Van Ree and Van Beukering 2016), strictly related to the Earth's geodiversity (Prosser 2013).These latter are able to satisfy human needs with the provision of rare-earth, geothermal resources, sites of cultural interest and other natural resources (Gray 2019).
Geodiversity also contributes to public health moulding unique landscapes in which to spend time outdoors (Crisp et al. 2021;Gordon and Barron 2013;Gray et al. 2013;Kaskela and Kotilainen 2017;Najwer et al. 2022).
In this context, the importance of studying the geodiversity of Italy emerged.The complex long-term geological evolution of Italian Peninsula coupled with the unique history of human civilization has resulted in a great variety of rocks and landscapes (Marchetti et al. 2017).The human presence has profoundly contributed to the shaping of Italian landscapes, leaving a clear sign in the land use as for example cultivated terraces, exceptional modifications of mountains and hills, widespread within the Mediterranean area (Cicinelli et al. 2021).
In spite of this, geodiversity in Italy has not yet been studied at national scale, since the attention was mainly focused on the analysis of geosites and geological heritage (ISPRA 2005).In this framework, our aim was to evaluate and map the geodiversity of the whole Italian territory (Fig. 1).It is extremely variable in terms of lithology and morphology as testified by the presence of Alps and Apennines chains, active and inactive volcanoes, lagoons and wetlands along 7.500 km of coastline and of climate variability due to its 1300-km elongated shape.
According to Sharples (2002), we considered the geology and topography the basic systems to identify the abiotic components of the environment useful for geodiversity assessment.The lithological units (expression of geological system) have a primary influence in defining landforms and soil forming.The topography, determined by past geomorphological processes (gravitative, fluvial, glacial, and periglacial, aeolian, coastal/marine and karstic ones), in turn controlled by climate (i.e., quantity and physical state of water that determines weathering, erosion and depositional processes), and tectonics (i.e., uplift, volcanoes, earthquakes), is suitable to map the current landform (expression of the morphological system).Furthermore, the elements of hydrological system (rivers, lakes, glaciers, wetlands) are considered as source of freshwater.
A semi-quantitative procedure based on the use of indicators and indices, which takes advantage by the map algebra algorithms to spatialize data, was defined.These algorithms are able to define and illustrate in proper maps the spatial variability of input data (lithological, morphological, pedological and hydrological) through matrices (grids recording indicators and indices compiled from input data) from whose sum the geodiversity map was derived.The methodological novelty is the use of a two steps procedure which made possible to preserve the spatial resolution of the lithological, morphological, pedological data providing a geodiversity map that faithfully reproduces the territorial features.For the hydrological data a procedure that made possible to preserve the natural elements that characterize the landscape (e.g., the longest rivers, largest lakes, etc.) was also implemented.
The procedure presented here, is also a valid approach for studies at regional and local scale.

Study Area
The Italian peninsula has an elongated shape extending into the Mediterranean Sea in a NW-SE direction.It is characterized by many islands; the largest ones are Sardinia and Sicily.The wide latitudinal extent (from 35°29' to 47°05' North), the altitudinal range (from sea level to over 4800 m) and the coastline length (ca.7500 km) result in very different climatic and environmental conditions in the country.The mountainous and hilly areas of the Alpine and Apennine chains dominate the landscape with their continuity and significant average altitude.Along the coasts are present several alluvial plains, the biggest is the Po Plain (constituting 70% of total plain areas) resulted from the accumulation of fluvial sediments deposited by the Po River system during the last one million years (Federici 2000;Fredi and Lupia Palmieri 2017;Marchetti et al. 2017).
The Alpine Mountain range is the result of the convergence and collision of the European and African (Adria) continental margins during Middle Cretaceous to Late Eocene.The Alps are a thrust belt constituted by two different mountain chains which developed in opposite directions ("double vergence") that are separated by a series of major faults (i.e., the Insubric Line).The "Alpine chain" is formed by a pile of crustal nappes of continental margin and metamorphic ophiolites with European vergence, overlapped northward since the mid-Cretaceous.The "Southern Alps" are a younger tectonic system consisting of Mesozoic sedimentary rocks similar to the Apennines, which since the Miocene has developed a southern vergence, i.e., toward the Po Plain.This large plain is an alluvial flat region resulted by an earlier marine and more recent fluvial sedimentation, fed by the Po River and its tributaries.Large part of the plain has a substratum of "buried mountains" which are the front of both northern Fig. 2 Schematic map of the peri-Mediterranean orogens and foreland areas (Cosentino et al. 2010, modified) Apennine and Southern Alps thrust belts (Fantoni and Franciosi 2010).
The Apennine thrust belt is the result of the collision of the western continental margin of Adria (i.e., the so-called African Promontory) with the Sardinia-Corsica block, mainly during Miocene-Pliocene time (Castellarin and Cantelli 2010).The structural edifice of the Apennines consists of a series of east-verging nappes, interposed between the back-arc Tyrrhenian basin to the West and the undeformed Apulian-Adriatic foreland to the East (Bonardi et al. 2009).The Ligurides, which include ophiolites and oceanic sedimentary rocks like radiolarites, originated from the so-called Ligure-Piemontese Ocean that has disappeared during oceanic lithosphere subduction and Alpine orogeny.
The Apennines are a segment of the peri-Mediterranean orogenic belt and can be subdivided into three principal sectors (Northern, Central and -Southern Apennines) bounded by regional tectonic lineaments (i.e., Sestri-Voltaggio Line, Ancona-Anzio Line and Sangineto Line).
The Northern sector is characterized by the abundance of flysch formations ranging in age from Cretaceous to Miocene, while the Central Apennine sector is characterized by the presence of large shallow-water carbonate platforms of Late Triassic-Cretaceous age, which constitute the highest mountain of the Abruzzi region (Gran Sasso, Maiella), followed by mainly pelagic cherty limestones and marls sedimentation (Parotto and Praturlon 1975;Castellarin et al. 1984).
The Southern Apennine sector is characterized by both large carbonate platforms and cherty-calcareous-pelitic pelagic basin of Triassic-Cretaceous age and Miocene flysch deposits (Bonardi et al. 2009).The Apulian region is formed by a subhorizontal Jurassic-Cretaceous shallow-water carbonate platform sequence, with the transition to the Adriatic deep-water basin exposed in the Gargano promontory; this region represents the undeformed foreland of the Apennine thrust belt.
The Calabrian-Peloritan Arc extends from the Sibari Plain to the northeast corner of Sicily, trough the Messina Strait.It is a segment of the Alpine chain that before the opening of the Tyrrhenian basin was located close to Sardinia.The arc consists of a pile of east-verging nappes including metamorphic basement and granites of Paleozoic age.At the eastern margin, relatively undeformed Miocene-Pliocene terrains cover the front of the various nappes.The island of Sicily (except the Peloritani Mountains in the NE sector) is the easternmost segment of the Maghrebian chain of north Africa and belongs to the northern continental margin of this continent.It consists of several south-verging nappes and a foreland area in the southeast corner (Iblei Mts.) (Abate et al. 1978).
The island of Sardinia is a fragment of the European continent, together with the Corsica Island (France).It separated from Catalunia (Spain) about 30 Ma ago and reached the present position about 18 Ma ago.Large part of the eastern side of the island consists of granites and Paleozoic rocks formed during the Hercynian orogenesis (Carmignani et al. 1992).The eastern side of the island is characterized by the presence of widespread volcanics, mainly of Tertiary age.

Materials and Method
The geodiversity has been assessed using lithological, topographic, pedological and hydrological data as illustrated in the workflow of Fig. 3.These input data have been transformed into grid (matrices with cells or pixels organised in rows and columns) useful to support the spatial analysis and map drawing.
These grids record the indicator and index values, which measure the spatial variability of the input data.An indicator gives a measure, generally quantitative, useful to simply illustrate and communicate complex phenomena, including trends and progress over time (EEA 2005).An index is an aggregation of indicators (OECD 1993).

Input Data
The source and the main features of input data for the geodiversity assessment, are listed below.
a) Lithological data.The lithology has been computed from the geo-lithological Map of Italy at 1:100.000 scale (Bucci et al. 2022).The authors derived this map by the union of 277 geological map at 1:100.000 scale of the Geological Service of Italy (ISPRA, Servizio Geologico d'Italia, http:// sgi.ispra mbien te.it/ geolo gia10 0k/).It reports nineteen lithological classes identified by adopting the dominant lithotype concept, in accordance with the Global Lithological databases (Hartmann and Moosdorf 2012), and considering also the geotechnical characteristics of the rocks.Taking into account the purpose of present work, we have modified the map by merging some classes (i.e., the "carbonate rocks" with the "marlstone" and the "non-schistose metamorphic rocks" with the "schistose metamorphic rocks") (Fig. 4a).Furthermore, the class of "mass wasting material", which includes landslides, was not taken into account and the small polygons of this class have been assigned to the lithology of area in which they are included.A total of sixteen classes were considered to calculate the spatial variety of lithology (Appendix, tab.1).
b) Morphological data.The morphological data was derived from the Digital Elevation Model (DEM) of Italy with cell of 20 m side (Fig. 4b) available at Ministry of Environment website (http:// www.pcn.minam biente.it/ viewer/ index.php? servi ces= dtm_ 20m).This model was obtained by using the elevation data reported in the topographic maps at 1:25.000 scale of the Italian Geographical Military Institute (IGMI).For the present study, the DEM c) Pedological data.The pedological map of Italy at 1:1.000.000scale (Costantini 2012) represented the source of soil data (Fig. 4c).The authors drew this map by processing the data recorded in the soil geodatabase of Research Centre for Agrobiology and Pedology.It is an update of database compiled by Mancini (1966).The soil map, showing the spatial distribution of main soils in Italy, is supported by a detailed legend which explain the soils features (Appendix, table 3).d) Hydrological data.According to Serrano and Ruiz-Flaño ( 2007), seas, rivers, lakes, springs and wetlands are the main elements of the hydrosphere system (Fig. 5a, f).The ability of these elements to shape the landscape has already been considered by the morphological analysis, so in this context we want to consider them as a natural resource.In this perspective, we have also included the glaciers, because they represent the main reservoir of nonliquid fresh water on the Alps, and the aquifers that retain one of most precious freshwater resources.
i) Seas -The Italian Peninsula is bordered by Ligurian, Tyrrhenian, Ionian, and Adriatic seas, (ISPRA, 2010) with a coastline about 7.500 km long.The coastline was considered an appropriate geographical element to identify the zones that can benefit from the sea resource.It was derived from the official administrative boundary published by Istat (https:// www.istat.it/ it/ archi vio/ 222527) (Fig. 5a).The spatial resolution cannot be univocally certified by Istat, because the data sources used to draw the coastline, mainly aerial photographs and cartographies, have different scales, usually equal or more detailed of 1:25.000,i.e., the scale used for morphological data.
iii) Wetlands -Wetlands consist of marsh, peat, natural and artificial ponds, permanent or temporary, with still or flowing water, fresh salty or brackish, including marine waters whose depth during low tide does not exceed six metres (Ramsar, 1971) (Fig. 5c).These Italian zones can be accessed through the Web Feature Service offered by ISPRA (http:// sgi1.ispra mbien te.it/ zoneu mide/ viewer/ indic ator.html) and are also available on RAMSAR web site (https:// rsis.ramsar.org/).In Italy, these zones, extending for 73.982 hectares, are 57 and encompass 15 regions.iv) Lakes -Italian lakes were acquired from dataset "lakes and other inland waters" available on the Ministry of the Environment website (http:// www.pcn.minam biente.it/ viewer/ http:// www.pcn.minam biente.it/ mattm/ visua lizza tori/) (Fig. 5d).It consists of hydrographic elements reported on the Regional Technical Maps, at 1:10.000 scale, integrated, for data lacking, with those reported on the IGM (Istituto Geografico Militare) topographic maps at 1:25.000 scale.v) Glaciers -The glaciers were recovered from the glacier inventory published by Paul et al. (2019) (available at https:// doi.org/ 10. 1594/ PANGA EA. 909133) (Fig. 5e).These authors mapped the Alpine glaciers, extended for about 1806 ± 60 km 2 , from Sentinel-2 images (Paul et al. 2020).Only the glaciers of the Apennine chain, Calderone Superiore and Calderone Inferiore, were digitized, checking their limits on Google Earth images, by the interactive map of Italian glaciers (https:// repo2.igg.cnr.it/ ghiac ciaiC GI/ ghiac ciai_ new.html), an output of the Italian Glaciers Database (Smiraglia and Diolaiuti 2015).vi) Aquifers -The aquifer map (Fig. 5f) describes the capacity of different lithologies to favour the recharge of aquifers and it depends on how easily rainfall can flow through rocks or unconsolidated sediments.It was derived from the lithological Map of Italy at 1:100.000 scale drawn by Bucci et al. (2022).The permeability classes were reported in Appendix, table 4.

Methodology
In the present work, a quantitative method based on the use of grids to discretize the study area (Pereira et al. 2013;Santos et al. 2017;Bétard and Peulvast 2019;Dias et al. 2021;Carrión-Mero et al. 2022) was adopted (Fig. 3a).
In geodiversity map processing, the effect of the choice of the grid size and the grid resolution on the accuracy of spatial modelling have been discussed by several authors (e.g., Hengl 2006 and references therein).Small size grids (derived by large-scale maps) are not appropriate for national or regional studies, as the crowding of shapes coupled with a too detailed legend would make the maps difficult to manage and read, on the counterpart, large size grids (derived by small scale maps) does not highlight the details needed to characterise the study areas (Pereira et al. 2013 and reference therein).Nevertheless, the choice of grid resolution is rarely based on the characteristic of spatial variability of the input data (Vieux and Needham 1993;Bishop et al. 2001).Indeed, in most GIS projects, the grid resolution is chosen subjectively, based on the sensitivity and experience of the researchers.
In our study, the suitable grid resolution (p) was estimated by using the concept of Minimum Legible Delineation (MLD) of the input data and the Scale Number (SN) of output maps.The value of MLD is considered of 0.25 cm 2 by Vink (1975): (1) According to (1) and setting the scale of the geodiversity map to 1:1.000.000(SN = 1.000.000), the pixel size (p) that allows to depict MLD is equal to 2.5 km.
A rule of thumb, that considers the value of p to be 10% of the grid cell size was used.Considering that p for the geodiversity map at 1:1.000.000scale is 2.5 km, a grid cell of 25 km side was adopted and named CS F , that stands for fixed cell size.In order to verify the reliability of CS F , we analysed the relationships between the cell sizes and extent of study areas, derived from scientific literature on geodiversity assessment (Appendix, fig.1).A linear correlation with a regression coefficient (R) of 0.92 and a standard deviation of ±1.5 km emerged between the cell size and extent of study areas, resulting in a cell with 30 km side for our study area characterized by an extension of about 300.000 km 2 , which is very close to the value of 25 km calculated with the rule of thumb.
The same procedure was also applied to the input data involved in the geodiversity assessment, whose scale varies from 1:10.000 to 1:1.000.000(Table 1) to calculate for each one the more appropriate cell size (CSv, variable cell size).For the morphological data at 1: 20.000 scale, the cell dimension derived from the role of thumb (0.2 km side) was not adopted because of computer processing power and time-consuming.Therefore, the elevation grid was firstly resampled into cells with 100 m side than the slopes were calculated and reclassified according to classes proposed by Steinke et al. (2016) All grids record indicator or index values (Fig. 3a), calculated from the input data, and show their spatial variability in thematic maps.They have been ranged into five classes, with the Natural break method (Jenks and Caspall 1971), varying from 1 (very low variety) to 5 (very high variety) except for the hydrological data for which "1" indicates very low and "5" very high contribution to freshwater resource.This last step is essential to make all data comparable and synthesizable (Lirer et al. 2010) in the geodiversity map.
The single indicator or index (hereafter, I) has been assessed by applying different algorithms according to both the geometric features (polygons, lines, or points) and attributes type (qualitative and quantitative) describing the input data.

Polygon Features
For polygon elements (i.e., geology, topography, soil) spanning the entire study area, a two-step procedure (hereafter, VM-A), based on both the "variety" (Ferrando et al. 2021) and "mean" zonal statistics algorithms, was applied (Fig. 3a).The "variety" has been performed by using a "CS V " grid to highlight as finest as possible the spatial variability of the input data.This function counts the number of geographic features with unique values within each cell (e.g., the number of different lithology classes).The more classes fall into a cell, the greater is its spatial variability and thus the contribution to the geodiversity assessment (Fig. 3b).The "mean" has been applied to the output of "variety" by using a "CS F " grid of 25 km side.This latter is an essential step to make all grids comparable and summarized into the Geodiversity map of Italy (Fig. 3a).This double step made it possible to preserve the greater spatial variability of data avoiding the loss of information (Fig. 3b).This result was evidenced by the comparison of the grids derived from the VM-A method with those calculated by using only the "variety" algorithm (hereafter, V-A) and grid with "CS F " (Fig. 3a, c).The VM-A procedure was not applied for the soil data (Fig. 3a) because they are in the same scale of the geodiversity map, in this case the "CS V" and "CS F " coincide.
For the hydrological data, the goal is to quantify the water resource in each cell.Regarding the aquifers, the weighted average (2) was applied to the rock permeability map to calculate the Permeability indicator (PM): where Pd is the degree of permeability and E p is the extent of the single polygon (Fig. 3a).
For lakes, wetlands and glaciers, two different procedures, resulting in as many grids, were followed (Fig. 3a): i) NV (Normalized value): the extent (E p ) of all elements (e.g., lakes) falling inside a single "CS F " were summarized and normalized according to the total area of geographic elements distributed on the whole study area (A t ) (i.e., sum of extension of lakes falling inside a cell divided by the total extent of all lakes present in the Italian territory) (Fig. 3a). (2) ii) MV (Maximum value): the maximum extent of polygon (E pmax ) between those falling inside the single "CS F " cell was considered as indicator value and assigned to it with the spatial join algorithm.This procedure preserves the lateral continuity of the geographical elements that shape the landscape (Fig. 3a).
The comparison of "NV" and "MV" procedures was illustrated in Fig. 6d.

Linear Features
The rivers were used as indicator of freshwater resource.The "NV" procedure (relation 3) was applied as follow: where L r is the length of single river falling inside a "CS F " and L t the total length of all rivers.
Similarly, the MV procedure was applied according to the relation (4): where L max is the length of the longest river falling inside a "CS F ".
The length of Italian coastline has been used as indicator of sea resource for the zones facing the sea.The intersection of Italian coastline with the "CS F " grid resulted in a single segment of coastline for each cell (L c ).The ratio (NV procedure) between the length of each segment and the total length (L t ) of the Italian coastline was calculated (5) and assigned to each cell of "CS F " grid (Fig. 3f).

Results
Four procedures ('VM-A', 'V-A', 'MV' and 'NV' -Fig.3a) have been implemented in our approach aiming at preserving the maximum detail of the input data and the elements shaping the territory as well as possible in the indicator/index maps.
The "VM-A" procedure ensured a better preservation of the shape of the input data in the indicator/ index maps.The difference between "VM-A" and "V-A" procedure has been illustrated for the lithology (Fig. 6a, d) and morphology input data (Fig. 6e, h). (3) The lithological variety derived from 'VM-A' procedure clearly preserved the shape of Alps and Apennines Mountain areas, characterized by the presence of a lot of lithotypes, in both output maps of 'VM-A' procedure (see black rectangles in Fig. 6b,c).Details are also well evident in Sardinia and Sicily islands, where it is possible to distinguish the areas dominated by few lithologies (see green rectangles in Fig. 6c) to those with several lithologies (black rectangles in Fig. 6b, c).On the counterpart, the map of 'V-A' procedure (Fig. 6d) shows the merging of the very high and high variability classes and the loss of classes from medium to very low (black rectangle in Fig. 6d).
Similar comparison was also made for the morphological data (Fig. 6e-h).The spatial variability of input data (Fig. 6e) is preserved in both maps of 'VM-A' procedure: the variety map with "CS V " (Fig. 6f) and the mean map with "CS F " (Fig. 6g).In both maps, along the Apennines chain it is possible to identify the relief with variable slope (orange and red colours), less complex morphologies (yellow colour) and lowland areas (light blue and blue colors).
Instead, the map drawn with 'V-A' procedure (Fig. 6h) displays the Italian peninsula dominated mainly by the high and very high spatial variability classes (orange and red colours), losing the morphological shapes reported in the input data (Fig. 6e).

Lithological Indicator (Lt) Map
The analysis of lithological variety map (Fig. 6c) shows a lithological variability (e.g., different rock types in the same cell) ranging from high to very high in the most part of the Alps, in the northern part of the Apennines (between Liguria and Tuscany regions), along the coastal area between Tuscany and Latium regions, in the southern Apennines between Campania and Calabria regions, in the area of Etna volcano, between Licata and Agrigento in Sicily, between Portoscuso-Villasimius and Fertilia-Isola Rossa in Sardinia.Medium variability characterises a large part of the central Apennines and the western sector of Sicily, while low or very low variability is mainly showed by plains (i.e., Po plain and Tavoliere delle Puglie).For localities, see Fig. 1.

Morphological Indicator (M) Map
The very high and high geomorphological variability classes mainly encompass mountainous areas of the Italian territory (Fig. 6g).These zones are topographically heterogeneous landscapes typified by the presence of mountain slopes, escarpments, terraces, structural reliefs, ridges, cliffs, canyons or very narrow valleys.On the counterpart, the low and very low geomorphological variability classes encompass the lowlands and the zones lying between mountain area and plains.In these zones, there is a greater spatial continuity of dominant landforms (floodplains, terraces, surface erosion, valley bottoms, and tabular surfaces).

Pedological Indicator (P) Map
The very high, high and medium pedological indicator classes encompass the hilly lands and they are due to the local variability of climate and lithology, the two main soil-forming factors (Fig. 7).
In Italy, soil regions located on hills are generally smaller than those located in mountains or plains zones (Costantini and Dazzi 2013).Medium variability classes characterize several sectors of Alpine foothills, only their eastern sector is typified by high variability class.The greatest part of Sardinia Island is characterized by the medium classes while the other zones pertain to low and very low pedological variability.

Hydrological Index (HI) Map
The Hydrological index map considers the contributes of both surface waters index (SWI) and terrain permeability indicator (PM) (Fig. 3a).For SWI, the comparison of "MV" (Fig. 8a-d) and "NV" (Fig. 8e-h) procedures (see paragraph 3 for detail) highlighted the advantages of using "MV" to preserve the spatial continuity of the elements shaping the territory (i.e., rivers, lakes, wetland, etc.).
In this way we have preserved the shape of Garda Lake (Fig. 8c), Po River and other similar features that would have been lost with the 'NV' procedure (Fig. 8e, h).For the latter procedure, a cell containing a little pond or a sector of the Garda Lake with the same extension pertain to the same variability class (Fig. 8g).
The SWI map summarises the contribution of the following five indicators drawn in as many maps (Fig. 9a-e): -Rivers (R), (Fig. 9a) -The most important resources are the Po River, followed by several Po tributaries, Tevere and Arno rivers.All other tributaries of the Po, Liri, Garigliano, Volturno, Sele, Ofanto, and Bradano rivers are resources of medium importance.Resources of lower importance are the Apennine rivers that flow into the Adriatic Sea, the "Fiumara" and "Rio" types of torrents that originate a short distance from the coastline.-Seas (S), (Fig. 9b) -The five classes of this resource are ranked according to the length of the coastline.Areas Fig. 7 Variety map of pedological data with a coastline length higher than 75 km fall into the very high resource class while those with a length of less than 14 km into the very low resource class.-Lakes (L), (Fig. 9c) -The most important resource is the Garda Lake, followed by all the lakes in northern Italy and the Trasimeno Lake in central Italy.Among the several lakes and lagoons characterising the coastal areas, the most important resource is the Venice Lagoon, while several minor lakes are along the coasts of Campania and Latium Regions.Both Sicily and Sardinia have few minor lakes, and both regions belong to the low and very low resource classes.
-Wetlands (W), (Fig. 9d) -Very high and high resource classes characterise the Venice Lagoon, the lake and marsh of Massaciuccoli in the Tuscany, the Saline di Margherita di Savoia in Apulia Region, the southern and western sectors of Sardinia near the city of Cagliari.
Smaller wetlands, belonging to the medium to very low resource classes (i.e.Volturno River mouth, Sele-Serre Persano WWF Oasis, Biviera di Gela, etc.) are widespread throughout Italy.-Glaciers (G), (Fig. 9e) -All glaciers are concentrated in the Alps, only two smaller glaciers pertain to the "very low resource" class are located on Gran Sasso Mountain (central Italy).
In the SWI map, the very high and high contribution (red and orange colours) to the water resource comes from the Alpine rivers, large glaciers in the Alps and wetlands close to the Po River mouth.In the central Italy, the high resource classes (orange colour) are due to the presence of Tevere and Arno rivers.A medium resource class (yellow colour) characterize the southern Apennine and the inner sector of main islands whose coastal belts are characterized by very high and high resource classes (Fig. 9f).The sum of PM and SWI resulted in the Hydrological index map (Fig. 9h).The very high and high resource classes characterize: i) the north-eastern sector of Italy, where rivers (Po and left tributaries), glaciers, lakes, and wetlands are present; ii) the area of Po River mouth and all the territory of Venice Lagoon; iii) the Arno and Tevere Rivers together with the neighbouring wetlands and zones with high permeability of rocks.The high variability of southern Apennines is mainly related to the high and very high permeability of rocks and to the presence of several rivers of medium resource class.Several sectors of Sicily and Sardinia islands fall within high and very high resource classes mainly due to the high permeability of aquifers and to the presence of several small ponds (Sicily) and wetlands (Sardinia).

Geodiversity Index Map
The geodiversity map of Italian territory resulted from the sum of the four indicators/indices maps (Lt, M, P, HI) (Fig. 3a), described in the previous sections.Two partial summary of maps were made to show the effect of individual matrix (lithology, morphology, soils and hydrology) on the assessment and mapping of Italy's geodiversity (Fig. 10a, c).
The first summary map displays the contribution of lithology and morphology maps (Fig. 10a).For the Alps, the very high and high variability classes of lithology and morphological indicators coincides (see red and orange colours in Fig. 10a).The Po Plain retain a very low variability according to the classes displayed in the two summarized matrices.Along the Apennine chain, the complementarity of medium variability class of lithology matrix and of very high and high classes of morphology matrix produced a wide area with high variability class that continues in both the islands (orange color in Fig. 10a).The southern Italy and the islands are Adding to this first output the pedological indicator map (Fig. 10b), it is possible to observe a reduction of very high variability on Alps, in the northern part of Apennine chain and in Sicily Island.As last step, the sum with the hydrological index map (Fig. 10c) has promoted the occurrence of medium class on the Alps and medium-low classes along the Tyrrhenian and Adriatic coastal zones.In the islands, the high class becomes middle class and the middle class turns low.

Discussion
In this work, a semi-quantitative approach suitable for the study of geodiversity at national scale was proposed and tested to the Italian Territory that is characterised by a complex long-term geological evolution and unique history of civilization.
Geological, morphological, pedological and hydrological input data were processed with diverse procedures conditioned by the type (continuous or discrete) of input data.The The two-stage procedure, "VM-A," applied to continuous data (i.e., geology, morphology, and soil), has made it possible to preserve the highest graphical resolution of input data (spatial detail of geographic elements) into the single indicator/ index maps (Fig. 6c, g) considering both their scale and the extent of the Italian territory.This result is clearly pointed out by the comparison with the indicator/index derived with the one-step procedure "V-A" that has calculated the variety of input data directly with cells of 25 km side (Fig. 6d, h).The two-step procedure has not been applied to the soils because the source data map has the same scale of the output geodiversity map and therefore the cell with a side of 25 km was already the optimal dimension to represent the data in detail (Fig. 7).
The procedure "MV", used for discrete data (i.e., hydrological data), has allowed the retention of shape of environmental elements modelling the territory.It has made it possible to preserve the main rivers (Po, Adige, Tevere, Arno, etc.), lakes (Garda, Como, etc.) and wetlands (Fig. 8a,  d), which would have been overlooked by the normalisation procedure (NV) (Fig. 8e, h).
These procedures ("VM-A" and "MV") allowed us to preserve the diversity richness of the parameters chosen for the geodiversity study while maintaining a link to the input data.Many authors have generated large area (i.e., regional/ national scale) geodiversity maps by proposing different methodologies and dealing with similar issues (Benito-Calvo et al. 2009;Carrión-Mero et al. 2022;Chrobak et al. 2021;de Paula et al. 2021;Pereira et al. 2013).One of the main limitations of the grid analysis is related to the scales of the input data, since larger scales show greater resolution than smaller ones and they could affect the results of spatial analysis (Marceau 1999).Finding a single grid size able to extract all the essential information from the data and turn them into comparable matrices may not be possible when the scales of such data are too different from each other, as in our case study (from 1:25.000 to 1:1.000.000).In this context, the methodology we proposed has proved to be a valid approach.
The comparison of indicators/index maps (Figs. 7, 6c, 6g, and 9h) with the partial summaries' maps (Fig. 10a, b) made it easy to identify the different contribution of single input data to geodiversity assessment (Fig. 10c) and to distinguish between the most geo-diverse areas for the hydrological aspect from those characterized by a high degree of variety for lithological, morphological or pedological features.
The Geodiversity Map points out the distribution of high and very high geodiversity classes along the Alpine and the Apennine Mountain chains and of very low and low classes in the coastal and alluvial plains (Fig. 10c).Morphological variability plays a key role in geodiversity classes distribution, the highest values characterize Alps and Apennines chain and the low -very low values the alluvial plain, as can be noticed for the Po and Tavoliere delle Puglie Plains, for which a low lithological variety, due to the prevalence of alluvial deposits, dominates (Fig. 10c).The soils system increases the variability of several site along the Apennine chain while the hydrological index increased the geodiversity of Alps and of Tyrrhenian and Adriatic coastal zones.

Conclusion
The geodiversity assessment is an arduous task in the light of the many components of the abiotic system and its operationalisation can be complicated by the need to integrate data of different types and resolution.
However, geodiversity is a valuable tool for environmental management exposed to extensive changes and increasingly intensive exploitation of natural resources.In this context, this research has proposed a GIS methodology, based on the use of indicators and indices to assess the geodiversity, at national scale and applied to the Italian Peninsula.
This methodology has processed four gridded indicator/ index maps that preserve both the maximum spatial resolution enabled by the scale of geology, geomorphology and soil input data and the elements structuring the physical environment for hydrological data.
The visualisation of the geodiversity indicators is userfriendly, the GIS procedure also allows to interact with the input data and to view step by step the indicators and indices calculated during the workflow.Furthermore, the indicator/ indices maps, all classified in the same number of classes, are easy to read and comparable to each other.
The main benefit of a geodiversity map at national scale is the possibility to have an overview of the geographical distribution of geodiversity classes.The availability of such maps can promote collaborative spirit between stakeholders and the overcoming of sectoral approaches in environmental management.It also provides the tools to identify a priority list of areas to be preserved and properly managed to ensure future generations the geosystem services we enjoy today.
The proposed methodology can be exported in other national territories.It can be also adopted for geodiversity study at regional and local scale even if more detailed input data are required.This type of documents will be crucial in the next feature due to the increasing importance of geodiversity in environmental conservation and management under the climate change effects.Furthermore, this map, combined with geoheritage and geosite maps, can be used for environmental education, because it provides stakeholders the basic tools for the knowledge of the relevant geological and geomorphological aspects of an area.It can provide an overview of abiotic nature elements that can be attractors for tourism and useful for the assessment of natural resources, serving as a basis for identifying areas of exceptional natural value (Chrobak et al. 2021) and for setting geopark, as already experienced in Brasil (de Paula et al. 2021).
The GIS framework will be updated by future research focused on the mapping of other geological and natural resources (e.g., mineral deposits, building stones, geothermal resources, oil and gas, fossils…) and on the landscape changes influenced by ongoing climate change.

Fig. 1
Fig. 1 Geographical location of the study area (a) and general framework at European scale (b) are displayed.The places cited along the manuscript and name of Italian Regions and listed on the left side of figure and reported in the map of Italy

Fig. 3
Fig. 3 Flow chart illustrating the methodology proposed in the present work (a).The white boxes indicate the input data, the green boxes are the pre-processed data, the cyan boxes are indicators and the blue boxes are indices.The grey boxes display the algorithms applied for indicator/index assessment.In the lower part of image are

Fig. 4
Fig. 4 The maps of lithological (a), morphological (b) and pedological (c) input data are shown.The map legends are in the Appendix (Tabs.1, 2 and 3, respectively)

Fig. 6
Fig. 6 Comparison between lithological (a) and morphological (e) input data and the variety maps calculated with two steps (VM-A) (b, c and f, g) and one-step (V-A) (d and h) procedure, respectively.

Fig. 8
Fig. 8 The maps display the indicator maps for wetlands (a), glaciers (b), lakes (c) and rivers (d) by applying MV procedure and for wetlands (e), glaciers (f), lakes (g) and rivers (h) by applying NV procedure

Fig. 9
Fig. 9 The grids of rivers (a), seas (b), wetlands (c), lakes (d), glaciers (e) are reported in the upper part of images.They are summarized in the Surface water index map (f).The Permeability indicator and hydrological index maps are reported in figures g and d, respectively

Fig. 10
Fig. 10 The figure illustrates the contribution of each index/indicator to the definition of the geodiversity map.In (a) the sum of the lithological (Lt) and morphological (M) data is showed, in (b) the pedological (P) data has been added and finally in (c) the hydrological (HI) ones

Table 1
Summary of input data, scale and source are indicated