**1. Introduction**

According to the European Environment Agency, "Land take is the process in which urban areas and sealed surfaces occupy agricultural, forest or other semi-natural and natural areas" [1] (p.117). "The most intense form of land take is soil sealing, which is an essentially irreversible process that leads to the destruction or covering of soils by buildings and other construction, and layers of completely or partly impermeable artificial material (asphalt, concrete, etc.). Soil sealing accompanies land take, but areas subject to land take are usually not entirely sealed" [2]. At the end of 2021, the European Commission approved the new European Union (EU) Soil Strategy for 2030, which highlights how healthy soils are essential for achieving the objectives of the European Green Deal concerning climate and biodiversity. The strategy defines a general framework and concrete measures to protect and restore soils to ensure their sustainable use [3]. One of the long-term objectives to be achieved by 2050 is to reach "no net land take" [3] (p. 3). The general framework of the strategy termed "land take hierarchy" concerns four consequential actions: avoid, reuse, minimize, and compensate. The first action, i.e., "avoid", aims at preventing further land take as much as possible. If land take cannot be prevented, then the second action, i.e., "reuse", should be implemented, with a view to reusing land that has already been urbanized or sealed; for example, through soil remediation or densification. If land take cannot be prevented and land cannot be reused, the third action should be looked at, to minimize the effects of land take by impermeabilizing land that is already in unfavorable conditions. If all the precedent actions cannot be taken, the fourth action provides for applying compensation and mitigation measures.

In relation to the EU-28 (which means the 27 EU Member States plus the United Kingdom), although land take decreased from 2000 to 2018, in the period 2012–2018 it

**Citation:** Isola, F.; Lai, S.; Leone, F.; Zoppi, C. Land Take and Landslide Hazard: Spatial Assessment and Policy Implications from a Study Concerning Sardinia. *Land* **2023**, *12*, 359. https://doi.org/10.3390/ land12020359

Academic Editors: Matej Vojtek, Andrea Petroselli and Raffaele Pelorosso

Received: 30 December 2022 Revised: 18 January 2023 Accepted: 24 January 2023 Published: 28 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reached the amount of 539 km<sup>2</sup> per year; 78% of land take takes place in agricultural areas, such as arable lands and permanent crops (50.5%), pasture and mosaic farmland (27.2%), forests and transitional woodlands (14.3%), and grasslands (6%) [4]. The main causes of land take are to be attributed to expansion of industrial and commercial areas, and to enlargement of residential zones and construction sites [2]. In Italy, in the 2020–2021 period, land take took a value of 69.1 km2, corresponding to an average value of 19 hectares per day. Only a tiny fraction of this growth in artificial surfaces was compensated by the restoration of natural areas, which equaled 5.8 km2, due to change from consumed soil to unconsumed soil through to the recovery of building sites, areas, and surfaces in cases where "reversible" land take [5] took place.

Moreover, land-taking processes entail several problems, such as the loss of multifunctional and fertile soils, biodiversity degradation and loss of ecosystem services (ESs) [6–8]. ESs are defined as benefits provided by ecosystems to human beings [9]. Since 1970, when the term "ecosystem service" was coined [10], this theme has increasingly been studied and analyzed in conceptual terms [11,12], in terms of classification [9,13], and in relation to their assessment and mapping [14,15]. Moreover, ESs provide protection against hydrogeological hazards [16]. According to Notaro and Paletto [17], this protection can be direct or indirect, where the former concerns the defense against natural phenomena such as floods and landslides.

A landslide is defined as the movement of materials such as rock, soil, debris, and artificial fill downward and outward along a slope [18] when the forces of gravity exceed the slope resistance [19]. Moreover, although landslides occur primarily in mountainous areas, this phenomenon may happen in low-relief zones [20]. The United States Geological Survey (USGS) [20] identifies three types of landslide causes: geological, morphological, and human. Indeed, the drivers of landslides may be either natural, caused by the intrinsic properties of rocks and soils or by physical processes, such as heavy rains and seismic activity, or artificial when human activities bring about changes in slope stability, as happens through deforestation or excessive soil sealing [21].

From this perspective, landslides are strongly connected with land use and land cover dynamics [22] and, particularly, with human-driven processes [23]. Human-induced activities, such as land use/land cover changes, may alter vegetation structure and modify soil characteristics and hydrogeological processes [24,25]. Indeed, despite the slowness that characterized geological and geomorphological changes, land use/land cover changes can occur in a short period of time due to their high dynamicity [26,27]. Land use/land cover changes can influence landside events in terms of frequency and spatial configuration due to their potential negative impacts on hydrological and mechanical processes involving soils [28,29].

The relation between land use/land cover changes and landslides has been studied by various authors [23,30,31]. Hao et al. [23] investigated the extent to which the landslide disaster occurred in 2018 in Kerala, India, was influenced by land use/land cover changes through a comparison between land use/land cover changes before (2010) and after (2018) the disaster. Pisano et al. [30] investigated how land cover changes influenced landslide susceptibility in the past and how they might influence future events, by carrying out a landslide susceptibility analysis implemented through a spatial multi-criteria evaluation in relation to three past land cover maps (1954, 1981, and 2007) and three future scenarios (one in 2030 and two in 2050) in the Rivo Basin, Italy. Muñoz-Torrero Manchado et al. [31] studied the influence of deforestation and related agricultural activities on landslide susceptibility using remote-sensing techniques and free satellite data in Nepal.

However, although several authors have studied the relation between land use/land cover changes and landslide hazard, the relation between land-taking processes and landslide hazard is still under-researched. Therefore, this study aims at analyzing the relations between land-taking processes and landslide hazard in order to understand to what extent land-taking processes increase landslide hazard through a regression model that relates the level of landslide hazard to a set of land cover variables that includes artificialized

land; that is, land taken up through urbanization processes, and a set of covariates that represent the land cover types associated with the LEAC (land and ecosystem accounting) classification. The methodological approach is implemented in the Sardinia Region, Italy.

The study is structured into six sections as follows. The second section describes the study area, the methodological approach used, and the input data for the regression model (landslide hazard, LEAC land cover groups, geological characteristics and elevation). The results are presented in the third section and discussed in the fourth section. The fifth section provides recommendations and implications for spatial planning policies stemming from the results, while the sixth section provides concluding remarks and future directions of the research.

#### **2. Materials and Methods**

This section is structured into three subsections. The study area is described in the first subsection. In the second subsection, the discrete-choice Logit model is used to estimate the relations between land-taking processes and landslide hazard. The third subsection describes input data for the regression model.

#### *2.1. Study Area*

In Italy, a sectoral planning tool termed PAI, an acronym of "piano di assetto idrogeologico" (whose word-for-word translation would be "hydrogeological setting plan"), must identify areas prone to natural hazards; i.e., both landslide and flood hazard, and define measures to reduce their magnitude and to prevent or mitigate their impacts. The PAI can be regarded as part of the broader river basin managemen<sup>t</sup> plan envisaged by the Water Framework Directive (WFD) [32], and its responsibility lies with ad hoc established competent authorities. Accordingly, Italy has been divided into eight river basin districts [33,34], each having its own competent authority.

One of such seven districts coincides with Sardinia, an Italian island around 24,000 km<sup>2</sup> in size and the second-largest island in the Mediterranean Sea Basin. Sardinia is further divided into seven subdistricts [35], one of which, the so-called "Coghinas-Mannu-Temo" subdistrict (hereafter, CMT), is the area chosen for this study (Figure 1).

The reasons for choosing the Sardinian CMT subdistrict for this study are twofold. First, Sardinia is included in the CORINE Land Cover (CLC) inventory coordinated by the European Environment Agency under the Copernicus program of the EU; this makes it possible to retrieve a regularly updated series of land cover maps, of which the most recent one refers to the year 2018 [36]. Second, a comprehensive and detailed spatial assessment of landslide hazard and risk concerning the whole CMT was carried out and officially validated in 2014 [37] and it is publicly available from the regional geoportal.

Located to the north-west of Sardinia, with an area of 5575.5 km2, CMT stretches over more than one-fifth of the island and it comprises around forty watersheds, of which the largest and most important ones are the four ones from which its name originates, i.e., Coghinas River, Mannu River, Mannu River in Porto Torres, and Temo River (Figure 1, panel C). The prevailing morphology is hilly, heavily marked in the southern border by the Marghine-Goceano mountain chain and by the Mount Limbara rocky granitic massif to the east, with the exception of the Nurra coastal plain to the north-west and the smaller plains around the mouths of the rivers Coghinas and Temo (Figure 2, panel A). As in all of Sardinia, in CMT the climate is typically Mediterranean, mostly lower meso-Mediterranean, but with coastal areas included within the upper thermo-Mediterranean zone and mountain chains in the upper meso-Mediterranean [38] (Figure 2, panel B). As for vegetation, according to the study by Bacchetta et al. [39], more than 41% of the CMT host species belong to the Sardinian thermo-meso-Mediterranean cork tree series, while the Sardinian oak tree series and holm oak tree series occupy around 11% of the CMT each and the other vegetation series take lower percentages (Figure 2, panel C). Geological instability is diffuse in the study area, where landslide events have been recorded for decades: the Italian landslide inventory (IFFI [40,41]) has documented 398 landslide events occurring up to 2007 in the

CMT. The most prominent category is that of diffuse falls or topples (228 events), followed by simple falls or topples (95 events); third comes diffuse superficial instability (42 events), followed by rotational or translational slides (16 events). Very small numbers concern the other categories; due to the geological and geo-lithological characteristics of the study area, no flow events have been reported in the study area (Figure 2, panel D).

**Figure 1.** The eight WFD river basin districts in Italy (panel **A**), Sardinian watersheds (panel **B**), and the Coghinas-Mannu-Temo subdistrict with its watersheds (panel **C**).

#### *2.2. Regression Model*

The relation between landslide hazard (LH) and the size of land taken up (L\_TAKE) is assessed through a linear regression model that uses the LEAC land cover groups as explanatory variables, whose detailed definitions are given in Section 2.3. The covariate representing land take is one of the LEAC groups, namely the variable associated with the artificialized land LEAC group. Dependent and explanatory variables refer to the elements of a 300 m square grid that overlays the study area, and are measured as their percentage share of a grid cell. The model operationalizes as follows:

$$\text{LHL} = \alpha\_0 + \alpha\_1 \text{L\\_TAKE} + \alpha\_2 \text{ARA} + \alpha\_3 \text{PMF} + \alpha\_4 \text{FOR} + \alpha\_5 \text{GRSH} + \alpha\_6 \text{DEPOQ} + \alpha\_7 \text{VOLSE} + \alpha\_8 \text{ELEV} + \alpha\_9 \text{DIN} + \alpha\_{10} \text{FOR} + \alpha\_{11} \text{FOR} + \alpha\_{12} \text{FOR} + \alpha\_{13} \text{FOR} + \alpha\_{14} \text{FOR} + \alpha\_{15} \text{DFOL} + \alpha\_{16} \text{DFOL} + \alpha\_{17} \text{DFOL} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{10} \text{FOR} + \alpha\_{11} \text{FOR} + \alpha\_{12} \text{FOR} + \alpha\_{13} \text{FOR} + \alpha\_{14} \text{FOR} + \alpha\_{15} \text{FOR} + \alpha\_{16} \text{FOR} + \alpha\_{17} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{18} \text{FOR} + \alpha\_{19} \text{FOR} + \alpha\_{19}$$

where the dependent variable and covariates are identified as shown below:


**Figure 2.** Some features of the Coghinas-Mannu-Temo subdistrict: elevation (panel **A**); climate zones (panel **B**, based on [38]); simplified vegetation series (panel **C**, based on [39]); documented landslides events (panel **D**, based on data from the IFFI project [40]).

The estimates of the coefficients of the multiple linear regression show the correlations between landslide hazard and the land covers of the LEAC taxonomy and, in particular, the interdependence of LH and the size of land take.

The use of a multiple linear regression is motivated by the fact that prior assumptions are not available as regards the functional form of the relations between dependent and explanatory variables, which is consistent with several studies aimed at identifying the interdependence between spatial variables [42–45]. From this point of view, a spatial phenomenon, related to *n* variables, represented by a surface in an n-dimensional space whose equation is unknown, can be approximately detected, in each of its points, by the tangent hyperplane. The linear equation estimated through the regression model, which relates dependent and explanatory variables, identifies the tangent hyperplane in a small neighborhood of a point of the surface, and in such neighborhood it represents an approximation of the unknown equation of the surface [46,47]. As a consequence, model (1) represents the trace of a hyperplane on a surface in a ten-dimensional space, which reports the correlations between LH and the nine covariates listed above.

The covariates DEPOQ and VOLSE control for the ground substrate, by considering if, and to what extent, landslide hazard is influenced by the conditions of the substrate, which in the study area is mainly featured by cohesive and compact rocks such as volcanic sedimentary successions (VOLSE) and, secondly, by deposits from the quaternary era (DEPOQ), i.e., loose incoherent materials. ELEV controls for the altitude impact on landslide hazard. If the estimates of their coefficients in (1) are significant, this will entail that substrate and elevation are likely to influence LH, at least to some extent. The magnitude of the coefficients will show the size of the impacts, in terms of increase or decrease in the landslide hazard measure.

The sign of ELEV is expected to be negative, since in the study area, on average, landslide hazard conditions are more frequent in lowlands rather than in mountainous locations, as further discussed in Section 3.3, whereas the expected signs of DEPOQ and VOLSE are positive and negative, respectively, since it is intuitively likely that LH will increase as long as the substrate incoherence and looseness increases, and the other way around.

The variable HGLAGGED represents the spatially lagged values of LH, and controls for spatial autocorrelation of the dependent variable in model (1). The HGLAGGED definition is based on the methodology implemented by Zoppi and Lai [48], which builds on Anselin's studies [49,50].

Moreover, a *p*-value test is used to check the level of significance of the estimates of the coefficients of model (1).

#### *2.3. Input Data for the Regression Model*

The dependent variable and covariates needed to feed into the regression model (1) were calculated with reference to a 300 m vector square grid that covers all of the CMT subdistrict (Figure 3) and comprises a total of 62,231 cells, using three main input spatial datasets listed in Table 1.

**Figure 3.** The 300 m vector square grid used in this study: extent of the grid with reference to the Coghinas-Mannu-Temo subdistrict (panel **A**, no. of cells: 62,231), and detail (panel **B**).


**Table 1.** Input datasets used to compute dependent and explanatory variables.

#### 2.3.1. Landslide Hazard

In compliance with national law no. 183/1989, in Italy each competent authority for a WFD river basin district must approve, as part of the comprehensive river basin managemen<sup>t</sup> plan, its PAI, which only focuses on landslide and flood risks in the district. The PAI has a dual character: on the one hand, it is a knowledge-oriented and regularly updated tool, which provides for the spatially explicit assessment of flood and landslide risks and hazards, as well as of exposures, hence vulnerable infrastructure, buildings, and land. On the other hand, it is a legally binding plan, which contains provisions that restrict land uses and land transformations in areas prone to landslide or flood hazard: the higher the hazard level, the stricter the restrictions.

The Sardinian Basin Competent Authority approved a first version of its PAI in 2004; since then, the assessment of the hazard level has continuously been updated to integrate new studies in previously non-analyzed areas, or to revise locally the hazard level when a new infrastructure that mitigates natural risks is realized. Accordingly, the landslide hazard map has been revised 42 times so far, and the flood hazard map 59 times.

Landslide hazard levels in the study area were assessed in a study commissioned by the Sardinian Basin Competent Authority in 2011, whose documents are publicly available on the institutional website [37]. The outcomes of the study were approved in 2014 and, as far as the spatially explicit assessment is concerned, integrated within the 36th updated revision of regional PAI spatial dataset available from the regional geoportal.

As per the methodology used in the PAI, landslide hazard (HL) classes range in the 0–4 interval, as follows: no hazard: HL = 0; moderate hazard: HL = 1; medium hazard: HL = 2; high hazard: HL = 3; very high hazard: HL = 4. For each cell in the 300 m square grid shown in Figure 3, the independent variable LH in model (1) was calculated as the percentage of the cell's area having non-null landslide hazard (HL -= 0) in the vector data retrieved from the geoportal.

Moreover, LH's spatially lagged variable (HGLAGGED), included in model (1) as a covariate, was calculated using GeoDA (version 1.20) [51], developed by Dr Luc Anselin and his team, based initially at the University of Illinois at Urbana-Champaign and currently at the Center for Spatial Data Science, University of Chicago, United States of America.

#### 2.3.2. LEAC Land Cover Groups

The CORINE (acronym for "Coordination of Information on the Environment") land cover is one of the several spatial datasets made available by the EU through the Copernicus Land Monitoring Services, covering a total of 39 countries, i.e., both members of the European Environment Agency and cooperating countries, and regularly updated every six years using a standardized nomenclature, hence allowing for consistent classifications and measures across time and space, and enabling time-series analyses.

In this study the 2018 CORINE Land Cover vector map (CLC2018) was used. The map provides information on land covers, i.e., on the biophysical characteristics of the Earth's surface, through a hierarchical classification that comprises 44 classes at the third (and lower) level, 17 at the second level, five at the first level, with a minimum mapping unit equaling 25 hectares.

The CLC2018 was next reclassified so as to group the third-level land cover classes following the taxonomy used by the European Environment Agency for land cover accounts [52] and comprising eight groups. Information on how the reclassification was performed is provided in Table 2, whose last column lists the CLC classes that were assembled within a single LEAC group. For the purpose of this study, only five groups out of the eight listed in Table 2 were mapped because three (open space with little or no vegetation; transitional woodland and shrub; wetlands, water bodies and marine waters) are not relevant within CMT. Furthermore, the latter group was not relevant with respect to the aim of this study: indeed, the absence of any relationships between marine or inland waters and landslide hazards is quite straightforward. CLC classes listed in Table 2 that only contain one digit (for instance, "1.") refer to first-level land covers and comprise all of the second- and third-level land covers that detail the first-level one (for instance, 1.1.1, 1.2.1, and so on); likewise, classes containing two digits (for instance, "2.2.") refer to second-level land covers and comprise all of the third-level land covers that detail the second-level one (for instance, 2.2.1, 2.2.2, and so on). For instance, the "standing forests" group includes all of the sub-levels of the 3.1 class, which, in the study area, comprise three third-level land cover classes as follows: 3.1.1 (broad-leaved forests), 3.1.2 (coniferous forests), and 3.1.3 (mixed forests). More specific information on wood types and management can be found in another, and older, land use/land cover map produced in 2008 by the regional administration of Sardinia [53], which further details the CLC taxonomy up to the fifth level. According to this dataset, approximately 28% of the surface covered by the LEAC "standing forest" group in the study area was managed in 2009. Managed forests were almost completely made up of cork oak woods (27%), while negligible percentages concerned other types of managed woods, either broad-leaved (for instance, eucalyptus woods) or coniferous (for instance, pine woods, especially in coastal areas). While many cork oak woods are still managed for production purposes, especially in North-Eastern Sardinia [54], eucalypti and pine trees (both non-native species in the island) were planted mainly for swamp reclamation, slope stability, and erosion control in coastal dunes in the XX century; as of today, they are often unmanaged, to the extent that some have undergone a renaturalization process and have evolved into mixed forests, as a result of successional processes [55] and native species' regaining their spaces.

**Table 2.** LEAC groups and corresponding CLC classes.


Once a vector map of the LEAC groups was retrieved, for each cell in the 300 m square grid shown in Figure 3, the explanatory variables L\_TAKE, ARA, PMF, FOR, and GRSH, in model (1) were calculated as the percentage of the cell occupied by LEAC groups listed in Table 2, respectively, as nos. 1, 2, 3, 4, and 6.

#### 2.3.3. Geological Characteristics and Elevation (Control Data)

A 1:25,000 regional geological map of Sardinia was produced at the beginning of the year 2000 building upon geological data collected by the former regional agency for mines and quarries. The spatial dataset, available from the regional geoportal [56], identifies geological characters in compliance with the "CARG" national mapping program initiated in the 1980s by the Italian Geological Society. The taxonomy of the Sardinian geological map is hierarchically structured into five main classes and five levels ([57], pp. 49–108), and a simple reclassification was carried out in this study, whereby i., first-level classes only were considered and, ii., three main groups were retrieved by merging first-level classes. The three groups are as follows: i., quaternary deposits (also comprising lakes); ii., volcanic sedimentary successions; iii., intrusive complexes and metamorphic basements. Finally, for each cell in the 300 m square grid shown in Figure 3, the explanatory variables DEPOQ and VOLSE were calculated as the shares of the cell occupied, respectively, by quaternary deposits and by volcanic sedimentary successions. For any terrestrial cell in the grid, the share occupied by the third group is, fairly obviously, the difference between 100 and the sum of DEPOQ and VOLSE.

Elevation was retrieved from the 10 m resolution digital terrain model (DTM) available "off the shelf" from the regional geoportal [58]. The Sardinian DTM was produced in the early 2010s, based on elevation points and contour lines contained in the 1:10,000 regional technical map (CTR, acronym for the Italian "Carta Tecnica Regionale"). Because the production process was implemented in compliance with the national guidelines issued in 2009 [59], horizontal and vertical accuracy, though not explicitly stated in the DTM metadata, are as follows: horizontal tolerance: 2 m; vertical tolerance in open fields: 2 m; vertical tolerance in densely wooded areas (i.e., in areas where tree canopy cover is over 70% of the surface): 1 2 of the mean height of the trees. Next, for each cell in the 300 m square grid shown in Figure 3, the explanatory variable ELEV in model (1) was calculated as the average elevation in the cell.
