2.1

Study Area

The present study focuses on the pilot (case-study) Natura 2000 site, “Dikti: Oropedio Lasithiou, Katharo, Selena, Krasi, Selakano, Chalasmeni Koryfi” (Figure 1), located at the central-eastern part of Crete, in Southern Greece. The site, codenamed GR4320002, is a mountainous area, covering 34,007.16 ha with altitudes ranging from 280 to 2148 m. It hosts diverse forests and distinct habitat types, supporting biodiversity and providing essential ES. Most of the area consists of steep mountains and canyons, including a cultivated plateau. Dominant vegetation types include phrygana, pine forests and shrubs, with a small portion of the site covered by vineyards. There is also riparian vegetation with common perennial and annual plants. The study area is of high ecological importance as the flora and fauna include a variety of common species as well as rare and vulnerable endemic species, some of which are found exclusively in the area or in the mountains of Crete¹⁶.

Figure 1.

Study area.

2.2

Methodology Outline

To address the aim of our study, we implemented a two-step process as depicted in Figure 2. The first step involved image classification, through object-based image analysis (OBIA) and a machine learning algorithm to generate the ecosystem type map of the study area for the year 2022. The second step involved ES modelling and production of the corresponding ES maps. Subsequently, trade-offs and synergies between paired ES were identified. The outcomes of this study were then integrated into the ppGIS/webGIS LIFE-IP 4 NATURA platform. The methodology employed in this study is illustrated in Figure 2 and further elaborated below.

Figure 2.

Study workflow.

2.3

Data and preprocessing

The study utilized the PlanetScope Analytic Ortho Scene Product, sourced from the Planet Labs PBC platform. These products provide orthorectified, calibrated multispectral imagery, facilitating automated extraction of quantitative information. The imagery is processed to remove terrain-induced distortions, ensuring precise spatial representation. The PlanetScope Analytic Ortho Scene is particularly well-suited for processing tasks such as ecosystem type and land cover classifications. Radiometric corrections are also applied to the imagery to address known sensor artifacts and convert the values to at-sensor radiance¹⁷. The multispectral images used in this study were acquired in August 2022 and cover four spectral bands (blue, green, red and near-infrared) with a spatial resolution of 3 m. The images were also mosaicked to ensure complete coverage of the study area.

Moreover, we used geospatial data including the European Digital Elevation Model (EU-DEM) developed for the Copernicus programme, with a spatial resolution of 25 m¹⁸ and the pan-European database of river networks and catchments developed by the JRC’s Catchment Characterisation and Modelling¹⁹. Furthermore, we employed CHELSA V2.1²⁰ products including bioclimatic data on precipitation and evapotranspiration, at a resolution of 1 km. Soil data were sourced from the data hub of ISRIC - World Soil Information²¹ at a resolution of 250 meters, as well as from the European Soil Data Centre- ESDAC²².

All ancillary datasets were resampled to standardize their analysis at 25 m and reprojected to the Lambert azimuthal equivalence area (LAEA) projection. This procedure ensured uniformity across all datasets, facilitating their integration and analysis as data uniformity is a prerequisite for their integration into InVEST toolset.

2.4

Ecosystem Types Classification

For the purposes of this study, we produced an ecosystem type map of the study area for the year 2022, utilizing the high resolution PlanetScope imagery. We adopted a seven-class classification scheme according to the “Correspondence between Corine Land Cover classes and ecosystem types” table compiled by MAES (2013)²³, and considering the regional characteristics of the study area (Table 1).

Table 1.

Ecosystem type classification used in this study, based on the “Correspondence between Corine Land Cover classes and Ecosystem types” table established according to MAES (2013)23.

Object-based image analysis (OBIA) was carried out in conjunction with the Random Forest (RF) algorithm to classify the images. In OBIA, objects are created through segmentation representing significant conceptual information. These objects are defined by their spectral attributes, shape, size, and the relationships they exhibit with their environment (adjacency). The combination of these factors is leveraged to ascertain the identity and characteristics of each object²⁴. The objects are then classified into distinct classes employing classification algorithms, including machine learning algorithms. Machine learning algorithms are highly effective classifiers of EO data, especially for ecosystem type and land cover classification tasks. Their popularity stems from their ability to achieve high accuracy while imposing minimal demands regarding the statistical properties of the data²⁵. Among these algorithms, Random Forest, introduced by Breiman in 2001²⁶, stands out as an ensemble classifier where multiple decision trees are constructed based on random subsets of training samples and variables²⁷.

The RF classifier was trained through training samples from visual image interpretation and ancillary geospatial datasets,. The spectral information taken into account included the four spectral bands of the original images. Furthermore, in order to improve the classification performance and highlight any discrete coverage within the study area, three spectral indices (NDVI, NDWI, NDSI) were calculated (Table 2). We also included information about the topography of the area derived from the EU-DEM. Elevation, slope, and the Topographic Solar-Radiation Index (TRASP)²⁹ were incorporated into the analysis process. The classification process was performed using the Trimble eCognition Suite, an advanced image analysis software available for geospatial applications³⁰.

Table 2.

Spectral indices and formulas used in the classification model.

2.5

Ecosystem Services Models

To quantify and map the ES of the study area, we developed six ES models using the InVEST toolset. InVEST is a geospatial modeling framework that assesses the impact of land use change on ES and is available at the Natural Capital Project¹⁴. The ES analysed in this study were (i) carbon storage, (ii) water yield, (iii) maximum potential water retention, (iv) soil protection from erosion (avoided erosion) and (v) nutrient delivery (nitrogen and phosphorous). The ES were selected based on data availability and their relevance to the study area.

(i) Carbon storage

The carbon storage module in InVEST employs a simplified carbon cycle to map and quantify the amount of carbon stored and sequestered³⁴. The model estimates the current carbon storage in a landscape or the amount of carbon sequestered over time using maps of ecosystem types combined with stocks in four carbon pools (soil, dead organic matter, belowground biomass, and aboveground biomass)³⁵. For each LULC type, the model requires an estimation of the carbon amount in at least one of the four fundamental pools mentioned above, given in metric tons per hectare (t/ha). In this study we used the ecosystem type map described in Section 2.4 and the carbon storage was estimated for each ecosystem type considering all the four fundamental terrestrial carbon pools applying literature values^34,36,37 and expert’s opinion.

(ii) Water yield

Water yield is the total amount of water in an area that contributes to maintaining the hydrological balance³⁸ and ecological equilibrium. From an ecosystem perspective, it is defined as the storage of water in rivers, lakes and aquifers and is perhaps one of the services most directly related to human well-being. The InVEST water yield model calculates the annual water yield from a basin, primarily targeted for supporting hydroelectric power generation from reservoirs³⁵. Although hydropower constitutes a relatively small portion of Greece’s energy sector, the total annual water yield has implications for various potential services, including agricultural irrigation, drinking water provision, hydropower generation, and industrial water extraction. The model is based on the assumption of the water balance equation of the Budyko hydrothermal coupling³⁹, taking into account the average annual precipitation and the actual evapotranspiration⁴⁰. Therefore, the water yield is calculated as follows:

Where Y_x is the water yield (mm), AET_X is the annual actual evapotranspiration and P_x is the average annual precipitation (mm) in each grid cell x.

The model inputs for this study were the ecosystem types map, the average annual precipitation (P), the Annual Actual Evapotranspiration (AET), the root restricting layer depth, the plant available water content and river basins boundaries. We also used as input the biophysical table containing biophysical parameters such as the plant root depth, crop coefficient (k_c) and information about vegetation for each ecosystem type. Finally, the seasonality factor (Z parameter) was defined setting the value 15 according to literature values³⁵ and expert’s opinion.

(iii) Maximum potential water retention

The maximum potential water retention (S) was calculated using the curve number (CN) approach as:

Where CN is calculated as the following empirical relationship:

Where iVEG represents vegetation class, iPERM represents permeability class, iSLOPE represents drainage capacity class. All indices range from 1 to 5⁴¹.

(iv) Soil protection from erosion (avoided erosion)

For this study the indicator avoided erosion was calculated as:

where AER_i is the amount of erosion avoided on pixel i, and the difference between RKLS_i and USLE_i represents the benefit of vegetation and good management practices. Sediment retention is computed by finding the difference between potential soil loss (USLE) of the landscape and the maximum potential soil loss (RKLS) which assumes the landscape is bare.

Where USLE_i is the potential average annual soil loss, R_i (MJ mm ha^-1 h^-1 yr^-1)is the erosivity factor, K_i (t ha h ha^-1 MJ^-1mm^-1) is the soil erodibility factor, LS is s the slope length and steepness factor, C_t is the land cover management factor and P_i is the supporting practice factor. The model’s inputs includes the ecosystem type map, the EU-DEM, Erosivity (R), Soil Erodibility (Κ)⁴³, and the biophysical table including the cover-management factor (C) and the support practice factor (P). The values of C factor were obtained from a literature search according to^44,45. The P factor was defined as 1 assuming that no erosion-reduction practices are being done³⁵.

(v) Nutrient delivery (nitrogen and phosphorous)

The model aims to map nutrient sources. It describes the transportation of a mass of nutrition through space using a straightforward mass balancing technique. A map of the types of ecosystems and the corresponding loading rates are used to identify the sources of nutrient burdens. Then, nutrient loads can be separated into parts that are dissolved and those that are linked to sediment. These parts will be carried by surface flow and subsurface flow, respectively, and will cease when they come to a stream. Subsequently, delivery factors are calculated for every pixel by utilizing the characteristics of pixels that are part of the same flow path.

The model’s inputs include the ecosystem type map, the EU-DEM, river basins, average annual precipitation (P) and the biophysical table including biophysical parameters (nutrient loading, maximum nutrient retention efficiency, critical length). Values were defined according to literature ^46,47. Threshold for flow accumulation was defined as 1000 after trial and error, and other parameters were set as default^14,35.

2.6

Tradeoffs and synergies

For the trade-offs and synergies identification between paired ES, we performed the Pearson’s correlation analysis using R programming language. First, we applied the min-max normalization method to standardize the ecosystem services. The strength of correlation between paired ES was determined by the absolute value of the correlation coefficient. A positive coefficient signifies a synergistic relationship, a negative coefficient indicates a trade-off relationship, and a coefficient of zero denotes no correlation ⁴⁸.

3.1

Ecosystem types

Ecosystem type map of GR4320002 is depicted in Figure 3. The natural ecosystem types (Grass, W-F, H-S, SVL, Wet) cover the largest part of the area (29904.88 ha), whereas anthropogenic types (Crop, Urb), including cropland and settlements, cover a smaller area (4983.75 ha) (Table 3).

Figure 3.

Ecosystem type map of GR4320002 (Urb: Urban, Crop: Cropland, Grass: Grassland, W-F: Woodland and Forest, H-S: Heathland and Shrub, SVL: Sparsely vegetated land, Wet: Wetland).

Table 3.

Area (ha) per ecosystem type of GR4320002.

More specifically, heathlands and shrub is the dominant ecosystem type covering the 53.74% of the study area and woodland and forest cover the 30.19%. Cropland covers the 13.84% and all the other types the 2.23% of the area (Figure 4). Regarding ecosystem types spatial distribution, forest vegetation is distributed throughout the whole area. Crops are distributed at the central and southeast part of the area. Wetlands (Wet) and urban (Urb) ecosystem types are also found in the core area of the site as well.

Figure 4.

Cover percentages of ecosystem types of GR4320002 (Urb: Urban, Crop: Cropland, Grass: Grassland, W-F: Woodland and forest, H-S: Heathland and shrub, SVL: Sparsely vegetated land, Wet: Wetland).

3.2

Ecosystem services

Table 4 illustrates the ES values for each ecosystem type. Carbon storage is 4443.905 10³t, soil protection from erosion is 7701.046 10³t, water yield and maximum potential water retention are 579.236 10⁶ mm and 158.296 10⁶ mm respectively and nutrient delivery is 25.624 10³ kgr for nitrogen and 3741.231 10³ kgr for phosphorus. The main sources of ES are the H-S, W-F and Crop ecosystem types. More specifically the H-S and W-F are the main pools of carbon storage and maximum potential water retention with their combined contribution being 92.56% and 92.35% respectively. Regarding the ecosystem service protection from soil erosion crop show the lowest values. Forests exhibit the highest values for these three ES. Conversely, cropland covered areas present higher values regarding nutrient delivery and water yield accounting for 60.09% for nitrogen, 75.13% for phosphorus and 26,03% for water yield (Figure 5).

Figure 5.

Ecosystem services provision (%) per ecosystem type.

Table 4.

Ecosystem services per ecosystem type of GR4320002.

In terms of their spatial distribution, carbon storage and maximum potential water retention exhibit similar patterns, with lower values distributed in the central and southeastern parts and higher values distributed throughout the area, especially in the south. In contrast, water yield and nutrient delivery show opposite patterns, with higher values concentrated in the central area. Soil protection from erosion display high pattern variation, with the lowest values observed at the central part of the area where the core of the cultivated land is located. (Figure 6).

Figure 6.

Spatial distribution of ES at GR4320002, (a) Carbon storage, (b) Water yield, (c) Maximum potential water retention, (d) Soil protection from erosion, (e) Nutrient delivery- nitrogen, (f) Nutrient delivery- phosphorus.

3.3

Trade-offs and synergies

The trade-offs and synergies of 15 paired ES are depicted in Figure 7. More specifically, 10 ES pairs display a trade-off relationship, whereas 5 pairs display a synergistic relationship. Regarding trade-offs, CS/WY (-0.878) show the highest degree of trade-off followed by CS/ND-N (-0.640) and CS/ND-P (-0.633). The paired ES MPWR/WY, MPWR/ND-N and MPWR/ND-P show high degree of trade-off relationship as well (r > |0.5|). All the other pairs indicate a low degree of trade-off (r < |0.1|). Concerning synergies, the highest degree is observed for ND-P/ND-N (0.978) followed by CS/MPWR (0.659) and WY/ND-N and WY/ND-N (r > |0.5|). Low degree of synergy shows the pair WY/SPER (r < |0.1|). Compared with other ESs, soil protection from erosion (SPER) showed the lowest degree of trade-offs and synergies with the other ES whereas CS showed the strongest.

Figure 7.

Pearson Correlation coefficients of paired ES. (CS: Carbon storage, WY: Water yield, MPWR: Maximum potential water retention, SPER: Soil protection from erosion, ND-N: Nutrient delivery- nitrogen, ND-P: Nutrient delivery- phosphorus).