Ecosystem Services Assessment for Their Integration in the Analysis of Landslide Risk

Abstract

Landslides are disasters that cause damage to anthropic activities, innumerable loss of human life, and affect the natural ecosystem and its services globally. The landslide risk evaluated by integrating susceptibility and vulnerability maps has recently become a manner of studying sites prone to landslide events and managing these regions well. Developing countries, where the impact of landslides is frequent, need risk assessment tools to address these disasters, starting with their prevention, with free spatial data and appropriate models. However, to correctly understand their interrelationships and social affection, studying the different ecosystem services that relate to them is necessary. This study is the first that has been attempted in which an integrated application methodology of ecosystem services is used to know in a systematic way if the information that ecosystem services provide is useful for landslide risk assessment. For the integration of ecosystem services into the landslide risk evaluation, (1) eight ecosystem services were chosen and mapped to improve understanding of the spatial relationships between these services in the Guerrero State (México), and (2) areas of synergies and trade-offs were identified through a principal component analysis, to understand their influence on risk analysis better. These are extracted from the models of the ARIES platform, artificial intelligence, and big data platform. Finally, (3) the similarity between the risk characteristics (susceptibility and vulnerability, already mapped by the authors) and the ecosystem services assessment was analysed. The results showed that the ecosystem services that most affect the synergy are organic carbon mass and the potential value of outdoor recreation; meanwhile, the possible removed soil mass was the most important trade-off. Furthermore, the lowest similarity value was found between landslide vulnerability and ecosystem services synergy, indicating the importance of including these ecosystem services as a source of valuable information in the risk analysis methodologies, especially with respect to risk vulnerability.

1. Introduction

Natural disasters have increased in recent years because of climate change. It is estimated that globally in recent decades, natural disasters have caused the death of approximately 150,000 people annually, with economic losses of more than 50,000 million dollars annually [1]. Furthermore, the different natural risks have not only direct effects on human well-being but also cause indirect damage to numerous services provided by ecosystems [2]. Specifically, in Mexico, 92% of the damage caused by natural disasters is produced by hydrometeorological phenomena [3].

Of the natural hazards associated with external geodynamics, slope movements are the second most dangerous after floods [4]. The downslope movement of a mass of rock, earth, or debris is known as a landslide [5]. The safety factor defines the stability of a slope, that is, by the balance between the driving forces that drag the ground materials down the slope by gravity and the resistance forces opposed to this movement, which is resistance in shearing (cutting) of the slope material. If this safety factor is less than 1, it means that the driving forces are greater than the resistance forces, causing the soil materials to be dragged down the slope [1]. Landslides usually occur after a period of intense rain [4].

These slope movements depend on both conditioning and triggering factors. Among the first, referring to internal factors that can influence the occurrence of landslides, are the structure, geomorphology, lithology, slope, climatic and hydrogeological conditions, and geomechanically properties of the materials. Triggering agents are external factors that trigger the event itself and can be erosive (waves, melting ice, rivers); dynamic loads (earthquakes); high-intensity precipitation; and volcanic, biological, or anthropogenic activity (excavations, mining, and overloads). These factors increase the vertical pressure or reduce the cohesion of the materials, accelerating the collapse of the rocky surface [4].

The main causes of these landslides are growing urbanisation, deforestation, and an increase in the magnitude of rainfall or the frequency of extreme rain events. The removal of vegetation cover increases surface runoff and soil degradation. The approach of population centres to areas of marked instability and greater precipitation, due to climate change, in areas prone to landslides causes a greater risk in sloped areas and high erosion. Therefore, these landslides occur mainly in mountainous ecosystems with high rainfall, where there is seismic activity; low drainage areas; or deforested, urbanised areas with inadequate infrastructure or overexploited by agricultural practices [1].

The most vulnerable countries to this phenomenon are developing, as poor spatial planning and deforestation for agricultural practices generate significant soil erosion, favouring floods and slope movements [2]. Among the damages caused by landslides are the effects on communication routes, buildings, forests, crops, and the economy and population of the area [1]. Therefore, natural disasters negatively affect ecosystem services [6].

Ecosystem services are benefits provided by systems to human beings and are classified according to the MEA (Millennium Ecosystem Assessment) in provisioning services (wood, water, food), regulation (available carbon, water, erosion), cultural (aesthetic value, source of inspiration, tourism) and support (photosynthesis, water, and nutrient cycles). Through them, it is possible to connect the human economy with the ecological dimension [7]. Natural capital refers to the different ecosystems from which these services originate, understanding them as a “stock” from which the service flows are produced. This natural capital is theoretically self-regenerative as it is independent of human actions [8].

It should be noted that most ecosystem services are quantified in biophysical and economic terms, but the lack of information or inadequate assessment of the ecosystem services provided limits the strategies that can be implemented in an area at risk of landslides [9]. The inclusion of the evaluation of ecosystem services in the risk assessment approach provides a basis for the analysis of those services that are lost or affected by the different disasters. This intervention in terms of planning and distribution of resources can delay ecosystem degradation and its prevention, thus increasing human well-being [10] and reducing the actual natural hazard risk [11]. However, few studies focus on how humans benefit from natural risk mitigation based on ecosystem service approaches [12]. Improved ecosystem services supply can be a key tool to optimise the human profit from nature at once protected infrastructures and livelihood through the mitigation of damage and the well plan of the economic activities in the territory [13]. Although, at present, the studies focused on ecosystem services and landslide risk assessment are based on the study of some ecosystem services, most of them related to soil erosion mitigation and changes in land cover to avoid land management, which increases the structural impact [11,14]. Otherwise, few studies focus on landslides, and ecosystem services focus on some ecosystems in special forests or watersheds [14,15,16].

Different studies in recent years on landslides have focused on the mitigation of damage and prevention of phenomena; however, risk assessment has been more relegated due to the difficulty of quantitative evaluation, both of the danger and vulnerability [17]. Nevertheless, these risk analyses allow the identification of regions with the greatest losses when a natural disaster occurs, possibly targeting mitigation measures and improving the effectiveness of risk management policies [18].

In order to mitigate landslide hazards effectively, new methodologies are required to develop a better knowledge of them and their evaluation, providing a systematic and rigorous practice to support infrastructure management and slope engineering [19]. Risk assessment consists of analysing hazards, vulnerability, and exposure. It is the product of a phenomenon’s probability and consequences, thus defining its potential damage [18]. The hazard represents the potential occurrence zones of the phenomenon based on the spatial distribution of past natural disasters [18], and vulnerability refers to the circumstances and characteristics that make an element susceptible to damage caused by a hazard [20]. Therefore, the maps representing the danger are approached by identifying the triggering and conditioning factors for the natural disaster occurrence [18]. As for the vulnerability maps, they are focused on the identification of environmental and socio-economic variables [21]. Correct risk assessment needs to use geographic information systems (GIS) since they allow for analysing and representing a wide extension of the territory [22]. This tool, widely used for landslide risks, enables the creation of maps on the spatial distribution of slope movements, the areas that could be potentially affected, and the probability of their occurrence in a spatial distribution [21].

The main objective of this article is to advance systematic and rigorous knowledge of landslide assessment, including an analysis of whether the information contained in ecosystem services is necessary to understand landslide susceptibility or vulnerability better using technologies based on data science and artificial intelligence. The specific objectives are mapping ecosystem services in the Guerrero area, using the artificial intelligence and big data tool ARIES. The treatment of ecosystem services through a principal component analysis to obtain areas of synergies and trade-offs. The identification of the ecosystem services is most related to the risk analysis from the perspective of the different components of the study (susceptibility and vulnerability). The evaluation of the importance of the various ecosystem services in the area of Guerrero and the assessment of the usefulness of ecosystem services to improve the analysis of natural risks to landslides through a similarity analysis.

2. Materials and Methods

2.1. Study Area

The study area corresponds to the State of Guerrero. This is a Mexican federal state in the southern part of the Republic of Mexico with an area of 63,794 km², which corresponds to 3.2% of the Mexican territory. It is bounded by the states of Morelos and Mexico to the north, with Puebla and Oaxaca to the east, Michoacán to the west, and the Pacific Ocean to the south [23]. The state of Guerrero Fin is located on the Guerrero–Morelos old marine platform, consisting of a series of extensive limestone outcrops [24] and the tectonostratigraphic complex of Xolapa and Guerrero, which are in tectonic contact [25]. Consequently, this state is frequently affected by hurricanes from the Pacific Ocean (the most common) and the Atlantic Ocean [26], which triggered massive landslides that have severely affected the population and infrastructures in recent years.

2.2. Materials Used for the Present Paper

For risk assessment, the components of hazard and vulnerability are necessary. In terms of landslide risk, the hazard is characterised by susceptibility maps. The susceptibility to landslides map of the Guerrero area was obtained from Bruzón et al. [27] based on the combination of different factors and landslide inventory through machine learning. In this map, the susceptibility values were established through 12 various conditioning factors (slope, orientation, distance to the drainage network, density of the drainage network, normal curvature, accumulated annual precipitation, lithology, distance to lineaments, density of lineaments, distance to road, density of roads, NDVI and land cover) and a trigger (precipitation). The landslide reference data were provided from an inventory made from three orthophotos taken on 12 August 2014, showing the consequences of the landslide period due to the Ingrid and Manuel meteorological phenomena.

The vulnerability to landslide map was provided by Arrogante-Funes et al. [28] based on the combination of ecological and socio-economic values. Ecological values were characterised by biodiversity, conservation status, habitat fragmentation, and regeneration delay (composed of erosivity of precipitation, soil erosivity, slope, and soil protection). Socio-economic values were estimated through marginalisation indices, population density, and building density.

Typically, the vulnerability assessment is based on ecological and socio-economic values. However, it is necessary to characterise the ecosystem services to have a holistic view of the values of the specific area [29,30]. Thus, these ecosystem services intend to provide a tool for unionising ecological and socio-economic values [31].

Artificial Intelligence for Ecosystem Services (ARIES) is the dynamic tool used to integrate and evaluate ecosystem services based on the software k.LAB through different values of fluxes, productivity, and societal uses [32]. Moreover, it provides models of supply and demand ecosystem services using empirical data and models to create flows in greater detail.

The different ecosystem services used are from provisioning (crop production, value of water from forests, value of non-wood forest products), regulation (organic carbon mass, the occurrence of pollinator insects, potential removed soil mass, and retained soil mass caused by vegetation), and cultural (potential value of outdoor recreation). Previously, a selection of ecosystem services affected by landslides was carried out (Table 1).

The value of water from forests is extracted through the data of the models of the ARIES platform (https://aries.integratedmodelling.org/, accessed on 15 October 2022) Regarding the regulation of water flows, the ecosystem service represents the ability of vegetation to intercept the rainfall, contributing to surface runoff. The regulation of these water flows on the earth’s surface depends on the components that store water in the surface layers and are influenced above all by the vegetation and physical properties of the soil [33]. In order to understand this service, ARIES simulates both surface and groundwater flows and groundwater extraction through wells, which connect human beneficiaries (agriculture, industry, domestic use) with upstream sources and sinks throughout the area thanks to high-resolution slope and elevation data [34].

The value of non-timber forest products was extracted from the ARIES platform (https://aries.integratedmodelling.org/, accessed on 15 October 2022). Regression models must be implemented to quantify the value of non-wood forest products (excluding firewood) for their monetary value. The value of these products is calculated as a function of the natural logarithm of the population density, the logarithm of gross domestic product (GDP) per capita, and a variable of the interest of the continent from the Global Human Settlement (GHS) (https://ghsl.jrc.ec.europa.eu/, accessed on 15 October 2022).

Crop production values provided by Crop Provisioning belonged to the ARIES Project (https://aries.integratedmodelling.org/, accessed on 15 October 2022). The contribution to agricultural productivity is estimated as the ratio of natural inputs to the ratio of the sum of natural and human inputs in energy terms (EcoConCrop), following the methodology of [35]. Thus, these models quantify the biophysical quantities of ecosystem services and are used by the different economic units. For the economic values, specific values are used for each type of crop based on the product prices published by the FAO (Food and Agriculture Organization of the United Nations) [36].

The organic carbon mass map provides the values of organic carbon stock from the ARIES platform (https://aries.integratedmodelling.org/, accessed on 15 October 2022). Globally, soils are estimated to hold up to 2400 GT of carbon [37], representing twice the atmospheric store of carbon and three times that of vegetation. The different anthropogenic impacts vary the capacity of the soil to retain carbon. The most important carbon storage relationship has been with climate change due to GHG emissions from soil destruction. Regarding land use, biomass burning and erosion are the main responsible for the release of CO₂ [38]. In ARIES, sources of ecosystem service are considered those capable of storing carbon: vegetation and soil. Areas where carbon can be released due to fires, land use change, deforestation, or other disturbances are called ecosystem service sinks. Models for organic carbon include global tables of plant carbon sequestration and global soil carbon sequestration with specific spatial data from ISRIC (International Soil Reference and Information Centre) (Aguilar, 2017). Total carbon sequestered is the sum of carbon accumulated above ground and below vegetation plus carbon accumulated in the top 200 cm of soil.

Potential removed soil mass map values were estimated from the sediment regulation data of the ARIES project (https://aries.integratedmodelling.org/, accessed on 15 October 2022). These data are calculated through the RUSLE universal soil loss equation (R-factor for the universal soil loss equation), which represents soil loss in tons of sediment per hectare per year [39]. Sediment flux modelling provides a spatial map of the connections between potential sediment fluxes, potential sediment sources, potential beneficiaries of sediment deposition, adsorbed sediment fluxes, loss of sediment value, and denied sediment sources. Thus, the model calculates the amount of benefit or harm for transported sediment.

ARIES project provided the retained soil mass caused by vegetation map through the vegetation values through sediment regulation data (https://aries.integratedmodelling.org/, accessed on 15 October 2022). As in the previous case, to estimate the contribution of the vegetation to the retention or removal of soil as an ecosystem service, the RUSLE [39] must be calculated twice, the first using the existing land cover and the second changing all the bare ground land use. Thus, a map shows the avoided erosion due to the vegetation in each cell.

The values of the occurrence of pollinator insects ecosystem service map were estimated through the data of pollination models based on the ARIES project (https://aries.integratedmodelling.org/, accessed on 15 October 2022). In this way, the concept of pollination in this field produces especially explicit data marked by their supply and demand by pollinating insects, providing services based on the distance to bodies of water, suitability of nesting, land cover, the percentage of land agriculture, weather patterns (including solar radiation and temperature) [40], and the occurrence of flowers that can be used as food for pollinating insects. This is represented through the maps of eco-floristic zones. The data on these are extracted from the FOOD and Agriculture Organization (FAO), also used by the Intergovernmental Panel on Climate Change (IPCC). Thus, the pollination occurrence map is obtained, which will provide a balance on the supply and demand of the region. The occurrence of pollinator insect map values is calculated as the product of the sum of the weight of the crops that depend on pollination and their production [41]. Then, the dependency on crop production is normalised to obtain the value of the deficit or the surplus quantity. Once all this has been obtained from the production model, the demand is extracted to produce the final values.

Finally, the potential value of outdoor recreation map values was estimated from the recreation models of the ARIES project (https://aries.integratedmodelling.org/, accessed on 15 October 2022). The models for the recreation ecosystem services are inspired by the nature-based outdoor recreation ESTIMAP models developed by [42]. In the simplest model, the theoretical potential of outdoor recreational attractiveness is considered a function that multiplies the factor of naturalness and accessibility at each point based on the distance to attractive areas of natural interest, established as a Euclidean distance to the protected areas, mountainous areas, and bodies of water. In order to estimate the ability of people to reach the study area, travel time data to the main cities extracted from the maps of the European Commission and The World Bank were used (https://www.worldbank.org/en/country/eu, accessed on 15 October 2022). Additionally, the population density data were provided by the Center of International Earth Science Information Network (CIESIN) (http://www.ciesin.org/, accessed on 15 October 2022). Thus, physical flows within ecosystem services can be quantified from the perspective of nature and the human benefits obtained by paying attention to potential flows.

Table 1. Summary table of the ecosystem services used downloaded from the ARIES tool. With a summary of the methodology used and the main unit of measurement.

Ecosystem Service	Method	Units
Organic Carbon mass	Sum of aboveground, belowground biomass carbon storage, and soil organic carbon.	Ton/ha
Crop production	Sum of production of maise, rice, wheat, barley, rapeseed, sugar beet, rye, soybean, sugarcane, potato, and sunflower.	Ton/ha
Occurrence of pollinator insects	Components of insect occurrence related to weather factors and landscape structure are combined to produce the pollinator occurrence map.	0–1
Potential value of outdoor recreation	Recreation potential values follow the ESTIMAP implementation for the recreation opportunity spectrum, which reclassifies the landscape by recreation theoretical supply and proximity to people.	0–1
Retained soil mass from vegetation	The potential value (supply) of the sediment regulation ecosystem service is computed by calculating RUSLE twice, first using the best land cover data available, then changing all land cover to bare soil and differentiating the results to estimate the avoided soil erosion attributable to vegetation.	Ton/ha
Potential removed soil mass	This implementation of RUSLE uses methods to calculate LS, based on contributing area, grid cell size, aspect, and slope length exponents, to calculate K, based on soil organic matter and clay, sand, and silt fractions, and global studies for C and P factors based on land cover type.	Ton/ha
Value of water from forests	Model-based on a regression function of the monetary value of water services.	USD 2013(PPP)
Value of non-wood forest products	Model-based on a regression function of the monetary value of non-wood forest products.	USD 2013(PPP)

Table 1. Summary table of the ecosystem services used downloaded from the ARIES tool. With a summary of the methodology used and the main unit of measurement.

Ecosystem Service	Method	Units
Organic Carbon mass	Sum of aboveground, belowground biomass carbon storage, and soil organic carbon.	Ton/ha
Crop production	Sum of production of maise, rice, wheat, barley, rapeseed, sugar beet, rye, soybean, sugarcane, potato, and sunflower.	Ton/ha
Occurrence of pollinator insects	Components of insect occurrence related to weather factors and landscape structure are combined to produce the pollinator occurrence map.	0–1
Potential value of outdoor recreation	Recreation potential values follow the ESTIMAP implementation for the recreation opportunity spectrum, which reclassifies the landscape by recreation theoretical supply and proximity to people.	0–1
Retained soil mass from vegetation	The potential value (supply) of the sediment regulation ecosystem service is computed by calculating RUSLE twice, first using the best land cover data available, then changing all land cover to bare soil and differentiating the results to estimate the avoided soil erosion attributable to vegetation.	Ton/ha
Potential removed soil mass	This implementation of RUSLE uses methods to calculate LS, based on contributing area, grid cell size, aspect, and slope length exponents, to calculate K, based on soil organic matter and clay, sand, and silt fractions, and global studies for C and P factors based on land cover type.	Ton/ha
Value of water from forests	Model-based on a regression function of the monetary value of water services.	USD 2013(PPP)
Value of non-wood forest products	Model-based on a regression function of the monetary value of non-wood forest products.	USD 2013(PPP)

2.3. Methods

Introducing ecosystem services would bring a new dimension to risk assessment (risk maps). This new approach would integrate the main drivers of land use change, allowing the evaluation of synergies and exchanges between different ecosystem services, proposing a simple method based on the sum of other correlations.

For the integration of these ecosystem services, an alternative methodology to the one typically used in risk studies is proposed here. First, we identified the main ES that could be related to the landslide risk, and we evaluated the ES in Guerrero State in Mexico. Then, we applied a principal component analysis (PCA) using the raster layers of each ES map individually. Our methodology introduces all the raster ES maps normalised into a principal component analysis (PCA) to comprehensively assess the territory’s synergies and trade-offs [43]. Finally, we applied analysis of similarities using the first (synergies) and second (trade-offs) PCA components and the maps from susceptibility and vulnerability achieved in Arrogante-Funes et al. [28].

With the study of similarities, we intend to know whether there is a relationship between the information they provide us and the susceptibility and vulnerability maps [28] or not. The point is that if there is already a clear relationship between the indicators of the ecosystem services summarised from the PCA and the susceptibility and vulnerability maps, the inclusion of the ecosystem services assessment maps in our risk methodological framework would only repeat information when it would be unnecessary to enter it. This study aims to understand better the existing relationship between ecosystem services and the maps already made of vulnerability and susceptibility to include them in a future version of risk.

An overview of the process followed in the present study it is shown in Figure 1.

2.3.1. Exploratory Data Analysis

Exploratory data analysis detects the structure of the data [44]. This analysis is divided into parts: the study of the data structure, cleaning and filtering the data, and finally, a graphic study of the elements of interest [45].

First, we explored the data structure of each ES map representing everyone with the help of GIS software and analysed their spatial coherency.

Then, we normalised and standardised some ES (variables). In order to carry out the PCA analysis and the analysis of similarities, the standardisation of the different variables must be given. The data of the variables to be treated are placed at the origin to rule out bias; that is, they are normalised. The rest of the data are scaled to eliminate the differences in the measurement systems. The z-score standardisation minimises the deviation of the variables whose numerical contribution is greater when class segregation occurs. Thus, in this way, the mean and standard deviation of the data set are used to rescale the data so that the resulting variables have a mean of zero and a variance of one. This is used for normalisation and subsequent comparison of the different PCA analysis service maps and similarity analysis [46].

Regarding the case of the potential removed soil mass and the retained soil mass caused by vegetation, in addition to the z-score standardisation, a logarithmic transformation was applied due to the asymmetric distribution of the variable, producing highly biased data. This logarithmic transformation is used when the data have very large values and a skewed positive distribution and allows the variance of the data to be more constant. Thus, it reconstructs the plot of the data towards an exponential curve, whose matrix that describes the curve represents the relationship between the values. Furthermore, this logarithmic transformation allows you to create contour maps, with intervals based on a logarithmic scale, by filling the maps with a gradient with the same scale as the axes. Thus, it is possible to create accurate maps of the data set when these are extended in different orders of magnitude [47,48].

Then, we performed a univariate analysis of each ES, including a graphical analysis of the distributions. On the other hand, a multivariate analysis of Pearson’s correlation coefficient was conducted.

We applied a min–max scaler algorithm to perform the graphical analysis of the distributions. This algorithm lets us transform the data by scaling each variable to a linear pre-established range, and it is performed to allow the best interpretation of the distributions since negative values are eliminated by introducing an acceptable minimum of 0 and a maximum of 1, unlike with the z score normalisation method, which goes from −1 to 1 [49].

The z-score, the min–max scaler, and the logarithmic transformation, as well as the principal component test and the analysis of similarities tests, were performed with the help of the python Scikit-learn tool [50].

2.3.2. Principal Component Analysis

Principal component analysis methods are widely used when very large data sets with many variables are available because they allow the dimensionality to be reduced. Thus, by having many related variables (being able to account for redundant information), this method reduces their number by creating other transformed variables known as principal components. These principal components collect most of the capture of the original variability, thus reducing their dimensionality but preserving as much information as possible [51]. These principal components are linear combinations of the original variables. The first principal component is the direction along which the samples show the largest variation. The second principal component is the direction uncorrelated to the first component, along which the samples show the largest variation [52]. In order to analyse the reliability of the variables in the principal component analysis, the statistics of Bartlett’s test of sphericity, Kaise–Meyer–Olkin (KMO), and omega hierarchical asymptotic were calculated.

Bartlett’s test of sphericity is a statistical measure to check the redundancy between the variables associating the observed correlation matrix with the identity matrix, and the alternative hypothesis is that variables are correlated enough to perform a data reduction technique [53]. Bartlett’s test of sphericity showed in Equation (1), where det(R) is the determinant of the correlation matrix, N is the number of observations, and p is the number of variables.

$$\mathrm{Barlett}=-\log \left(\det \left(\mathrm{R}\right)\right)\times \left(\mathrm{N}-1-\frac{2\mathrm{p}+5}{6}\right)$$

(1)

The KMO test is a statistical measure to establish how suited data are for factor analysis, values less than 0.5 indicate the sampling is not adequate [54]. KMO equation showed in Equation (2), where r_jk is the correlation between the variable in question and another, and p_jk is the partial correlation.

$$\mathrm{KMO}=\frac{\sum {\sum }_{\mathrm{j}\ne \mathrm{k}}{\mathrm{r}}_{\mathrm{jk}}^2}{\sum {\sum }_{\mathrm{j}\ne \mathrm{k}}{\mathrm{r}}_{\mathrm{jk}}^2+\sum {\sum }_{\mathrm{j}\ne \mathrm{k}}{\mathrm{p}}_{\mathrm{jk}}^2}$$

(2)

Omega hierarchical asymptotic is a measure of internal consistency reliability of data that is the overall consistency of a measure. We selected the omega hierarchical asymptotic test as a reliability measure because, for this type, the data are more appropriate than others [55]. The omega hierarchical asymptotic test is a modification of the omega hierarchical test that uses an infinite-length test with a structure like the observed test. Omega hierarchical equation showed in Equation (3), where is the ratio between the true score variance divided by the sum of variances and covariances of the data.

$$\mathsf{\omega }=\frac{{(\sum {\mathsf{\lambda }}_{\mathrm{j}})}^2}{{(\sum {\mathsf{\lambda }}_{\mathrm{j}})}^2+\sum {\mathsf{\sigma }}_{\mathsf{ϵ}\mathrm{j}}^2}$$

(3)

The first principal component marks the maximum direction of variation on the projections, representing the most synergistic areas. The second principal component is a new vector that passes through the origin, perpendicular to the primary component, thus forming a plane with a larger and smaller than three-dimensional space. Therefore, the second main component is where there are trade-offs among the services [43].

2.3.3. Similarity Test

The similarity between the risk characteristics (susceptibility and vulnerability, already mapped by the authors of this study [27,28] and the ecosystem services maps were analysed. The application of the similarity test was important to see how integrating these maps into the risk analysis could add valuable information.

Cosine similarity measures the similarity between two vectors of an inner product space. It is measured by the cosine of the angle between two vectors and determines whether two vectors point in nearly the same direction. In the present study, we used the similarity cosine test. This test measures the similarity using the cosine of the angle between two vectors in a multidimensional space [56]. It is given by:

Consider two vectors of features, x and y; it is given by Equation (4):

$$\mathrm{similarity} \left(\mathrm{x}, \mathrm{y}\right)=\cos \left(\mathsf{\theta }\right)=\frac{\mathrm{x} · \mathrm{y} }{\left|\mathrm{x}\right| · \left|\mathrm{y}\right|}$$

(4)

This metric determines the numerical similarity between two vectors, calculating the angle between them in a space with dimensionality equal to the number of layers we included. The value similarity is maximum when the vectors are in the same direction (angle equal to 0°) and becomes minimum when perpendicular (angle equal to 90°). This similarity metric can recognise patterns. Therefore, if the similarity metric between vulnerability or susceptibility and synergies ES and trade-offs ES maps were low, it was valuable to include these in the risk assessment work.

3. Results

3.1. Exploratory Data Analysis: Ecosystem Services Maps of Guerrero State

We observed a certain spatial coherence from the spatial analysis of ecosystem services. Figure 2 shows the maps from the eight ES cartographies downloaded by the ARIES platform (https://aries.integratedmodelling.org/, accessed on 15 October 2022).

From the value of water from the forest ES map (see Figure 2a), we could extract how a higher value occurs in mountainous areas since these have a lot of vegetation, the highest rainfall, and plutonic soils with little porosity. As for the Crop production ES map (Figure 2b), we could see how the values were higher in areas with little slope and soil rich in carbon. Regarding the map of the importance of non-wood forest products ES (see Figure 2c), we observed that it occurs above all in arboreal jungle areas. Regarding the occurrence of pollinator insects ES map (see Figure 2d), we observed how the areas with less value of this ecosystem service were found in agricultural areas and more in the jungle areas with greater vegetation biodiversity. From the map of ES called Retained soil mass caused by vegetation (Figure 2e), we saw how the values were greater in cloud forest areas where the roots retain the soil, and the tree canopy stops the dragging by the rains. In terms of organic carbon mass ES cartography (Figure 2f), we saw how it was higher values in areas with less deforestation. The potential removed soil mass ES map (Figure 2g) showed higher values found in areas with the highest slopes. Finally, regarding the potential value of outdoor recreation (Figure 2h), we could see how in areas with greater accessibility and regions where natural values predominate, it is where it is greater, see, for example, coastal areas and mountain areas.

The Pearson correlation test values graph (Figure 3) between the different ES maps shows how their interrelationships are relatively low, with almost all values lesser than 0.25. We can find how there is a negative correlation between the value of water from the forest and the value of non-wood forest products. It makes sense because the smaller the fraction of forest, the lesser amount of water it will be able to retain. We also see a positive correlation between organic carbon mass and the potential value of outdoor recreation.

3.2. Principal Component Analysis

Figure 4 shows the synergies and trade-offs of ecosystem services for the study area. The largest amount of independent information corresponds to ecosystem services (PCA1) synergies followed by (PCA2) trade-offs.

Figure 4a shows that the most synergistic zones are in the northeast of the study area. These areas represent points where high values of synergies of the different ecosystem services converge and, therefore, areas of high conservation and management interest. Conversely, the least synergistic zones are in the centre and coast of the study area. Although low values of synergies show the result of concentrating few ecosystem services, they are scattered throughout the region.

Figure 4b shows the areas of greatest trade-offs of ecosystem services in the north and northeast of the study area. These trade-offs imply that the increase in one of the services causes the detriment of another. On the other hand, it is observed that the regions with lower trade-off values are those related to the coast.

As we can appreciate in

Table 2. Results of PCA based on ecosystem services.

Number of PCA	1	2
Explicability	0.41	0.31
Eigenvectors	194.46	146.88
Input layer	Synergies	Trade-offs
Organic carbon mass	0.35	−0.38
Crop production	−0.36	0.34
Occurrence of pollinator insects	0.13	0.47
Value of non-wood forest products	−0.47	−0.24
Potential value of outdoor recreation	0.39	−0.35
Potential removed soil mass	−0.25	0.50
Retained soil mass caused by vegetation	0.29	0.28
Value of water from the forest	0.46	0.34

, the two principal components related to synergies and trade-offs reach together 0.72 of explicability. The results of statistical analyses are for KMO 0.58, Bartlett analysis p value 0.0, and omega hierarchical asymptotic 0.41.

As we can see, the value of water from the forests, organic carbon mass, and potential value of outdoor recreation reach the highest values of synergies (Table 2). This may be because the higher vegetation density produces higher organic carbon mass and a greater potential value for outdoor recreation activity. Hence, it makes greater soil retention by vegetation that prevents erosion and a higher value of water from forests.

On the other hand, the results show that the occurrence of pollinator insects and potential soil loss obtain the highest trade-off values (Table 2). These trade-offs may be associated with areas with less vegetation, such as cultivated areas or pastures, where we can find a greater occurrence of pollinator insects linked with more potential removed soil mass. This results in less organic carbon mass, reduced non-timber products, decreased potential outdoor recreation value, and less soil retention by vegetation.

3.3. Analysis of Similarities

The values in the landslide vulnerability (Figure 5d) distribution graph are displaced to the right, implying that most of the State of Guerrero area presents a medium–high vulnerability. However, in the susceptibility distribution graph (Figure 5c), the opposite occurs, i.e., the values are distributed mostly to the left, meaning that low–medium susceptibilities predominate in the area. On the other hand, we see how the synergies (Figure 5b) and trade-offs (Figure 5a) are distributed in a similar way to the normal distribution. Finding a greater number of trade-offs than synergies.

Table 3. Value from the cosine similarity test.

First Characteristic	Second Characteristic	Cosine Similarity Metric
Synergies	Synergies	1
Trade-offs	Synergies	−0.0003
Landslide Susceptibility	Synergies	0.3587
Landslide Vulnerability	Synergies	0.0333
Landslide Susceptibility	Trade-offs	−0.0108
Landslide Vulnerability	Trade-offs	−0.1955
Landslide Vulnerability	Landslide Susceptibility	0.0786

present the value from the Cosine similarity test. When the value of this test is close to 1, it implies that the similarity between the two curves is maximum, but when the value is close to 0, it implies that the similarity is minimal. In the case that concerns us, a low similarity would be telling us that the information provided by the study of synergies and trade-offs of ecosystem services is not being taken into account in the analysis of susceptibility and vulnerability. Therefore, from the similarity test (Table 3), we observed how the highest similarity was found between the landslide susceptibility and the ES synergies achieving a value of approximately 0.36. On the other hand, we also observed how the similarity between the synergies and the landslide vulnerability achieved a low value (0.033). We also observed how the landslide susceptibility and the trade-offs were not very similar since the cosine similarity test reported a value of −0.01. The value achieved by the landslide vulnerability and the trade-offs was slightly higher than the value risen by the landslide susceptibility and the trade-offs pair (−0.19), but it was inverse.

4. Discussion

In this study, an analysis of the importance of ecosystem services in the integration of landslide risk was developed through a similarity test and a principal component analysis (PCA). In recent times, ecosystem services are beginning to be included in evaluating and managing natural risks [11,57,58]. The inclusion of ecosystem services in these analyses allows the integration of human well-being in the evaluation of natural risks, improves the characterisation of ecological and socio-economic values, and facilitates the development and implementation of policies, regulations, and programs [59]. However, the methodology to assess, map, and evaluate ecosystem services does not rise to address the biophysical and socio-economic dynamics at different multi-temporal scales [32,60].

Accurate assessment of ecosystem services is demanding in terms of data. That is why modelling platforms based on artificial intelligence that allow identifying, personalising, and connecting data have a special interest in addressing this challenge [32]. Moreover, the characteristics of ARIES as a semantic and automatic reasoning network allow us to connect distributed spatial data from multiple sources (quantitative, semi-quantitative, expert opinions) in the modelling of ecosystem services, providing the opportunity to advance in a framework for calculating and evaluating ecosystem services, multi-spatial, multi-temporal, customisable, scalable, and replicable [61].

Regarding the correlation of ecosystem services, those most related will be the value of water from forest and forest production [11], as is obtained in the study when presenting the highest correlation value, the value of water retained with the production of non-timber products. Similarly, it occurs with organic carbon mass and the potential value of recreation because people prefer natural environments to carry out their recreational activities, which in our case focus on coastal areas and forest areas, which have a large amount of organic carbon.

In this study, the synergies occur for services that increase and improve the value of several others [43]. The ecosystem service that has the most synergies is the value of water from the forests since this is affected by vegetation (influencing factor in crop production, pollination, and non-timber production) and the soil type for the infiltration (affecting organic carbon mass, soil removal, and retention). Consequently, the value of water from forests involves synergies with many ecosystem services and functions, which implies greater importance in the model.

Trade-offs increase disadvantages between ecosystem services [43]. In this way, the ecosystem services that produce high trade-offs are removing potential soil and the occurrence of pollinating insects. Removing the soil produces losses in organic carbon mass, crop production, and recreation values. Moreover, the loss of soil is generally considered significant and irreversible. A reduction in the occurrence of pollination insects is related to lower crop production and the demand for more agricultural areas to the detriment of other types of ecosystems.

Concerning the evaluation of ecosystem services as a tool in the assessment of the risk of landslides, we observed that it is a subject of growing interest in recent years [11,62,63]. Our results found that synergies and susceptibility are the two most similar variables. This means high susceptibility in areas with high synergy and vice versa. However, there is no close relationship between the vulnerability and the synergies or trade-offs, being more influenced when the vulnerability is of medium range and producing greater influence into trade-offs than synergies.

In this way, risk analysis can be used to identify regions that need protection where service values are high, but landslide occurrence is also high [62]. Lower risk impacts were found for risk in terms of crop production and pollination services and higher for soil retention and avoided potential erosion [62].

The value of water from forests and organic carbon mass is typically presented among the services most involved in risk analysis [64]. Within the risk analysis, the relative position of ecosystem services among them is highly important, which justifies the need to establish synergistic and trade-off points [65].

Concerning landslide risks, it is noticed that these severely affect forests, leading to the total elimination of the ecosystem services their produce [66]. These forests involve high organic carbon mass, non-timber production, the value of water, and soil retention. Thus, it was observed, as in the study by Alqhadi et al. [63], that the vegetation effectively produces a relationship with the different ecosystem services and how these are related to the areas of greatest susceptibility.

Regarding the weaknesses of our study, it was observed that the explainability of the first components (synergies) and second components (trade-offs) of the PCA analysis reaches only 72% of the variability of the sample. This might be because some services may not be identified, and multiple connections, feedback, and outcomes in various service patterns have not obtained much consideration [67]. On the other hand, the little development of big data and artificial intelligence in applications related to the environment and, in particular, ecosystem services make it difficult to export the results with the necessary efficiency to face the current challenges of the synthesis of information [68]. Modern open science infrastructure, public data, and model repositories are useful starting points, but without shared semantics and common standards for machine actionable data and models, our collective ability to build, grow and share a collective knowledge base will remain limited [69].

5. Conclusions

In this study, we presented an assessment of how ecosystem services can have a place in the analysis of natural risks, especially in the analysis of landslides. At present, ecosystem services are starting to be included as part of risk analyses in hazard and vulnerability measurement. However, there is no study on whether ecosystem services incorporate new information in risk analysis. Or on the contrary, there is a reiterated source of information with current measures in risk analysis.

Our results show how the inclusion of ecosystem services treated jointly in the form of synergies and trade-offs have a relatively low similarity compared to typical characteristics of risk analysis, in this case, susceptibility and vulnerability. This leads us to conclude that including ecosystem services in risk analysis is important as it gives us information not captured by the metrics typically included in risk analysis. On the other hand, it is observed that the greatest similarities are seen between susceptibility and synergy, which tells us that it is not covered by the metrics ecosystem services as a link between ecological vulnerability and socio-economic vulnerability is the path that provides the most information to the system. Several risk studies currently include ecosystem services in their vulnerability analysis, which is key to capturing the flow of services between nature and human well-being.

Finally, in this study, we included some insights into one of the problems we have to face in environmental studies and in the analysis of natural risks, which is the integration of new big data technologies and artificial intelligence, which slows down the synthesis of information and compromises efficient decision-making.