Remote Sensing and Phytoecological Methods for Mapping and Assessing Potential Ecosystem Services of the Ouled Hannèche Forest in the Hodna Mountains, Algeria

Abstract

Regardless of their biogeographic origins or degree of artificialization, the world’s forests are a source of a wide range of ecosystem services (ES). However, the quality and quantity of these services depend on the type of forest studied and its phytogeographic context. Our objective is to transpose the concept of ES, in particular, the assessment of forest ES, to the specific Mediterranean context of the North African mountains, where this issue is still in its infancy and where access to the data needed for assessment remains difficult. Our work presents an introductory approach, allowing us to set up methodological and scientific milestones based on open-access remote sensing data and already tested geospatial processing associated with phytoecological surveys to assess the ES provided by forests in an Algerian study area. Specifically, several indicators used to assess (both qualitatively and quantitatively) the potential ES of the Ouled Hannèche forest, a forest located in the Hodna Mountains, are derived from LANDSAT 8 OLI images from 2017 and an ALOS AW3D30 DSM. The qualitative ES typology is jointly based on an SVM classification of topographically corrected LANDSAT images and a geomorphic-type classification using the geomorphon method. NDVI is a quantitative estimator of many plant ecosystem functions related to ES. It highlights the variations in the provision of ES according to the types of vegetation formations present. It serves as a support for estimating spectral heterogeneity through Rao’s quadratic entropy, which is considered a relative indicator of biodiversity at the landscape scale. The two previous variables (the multitemporal NDVI and Rao’s Q), completed by the Shannon entropy method applied to the geomorphon classes as a proxy for topo-morphological heterogeneity, constitute the input variables of a quantitative map of the potential supply of ES in the forest determined by Spatial Multicriteria Analysis (SMCA). Ultimately, our results serve as a useful basis for land-use planning and biodiversity conservation.

1. Introduction

Natural ecosystems are essential support for the well-being of human populations and all living things through the multiple services they provide [1,2]. For example, mountains, as the “water towers” of the world [3], are ecosystems that play a crucial role in the well-being of human populations, more than 50% of which depend on freshwater that is captured, stored, and purified in mountain regions [4]; this is especially true in arid and semi-arid areas, such as those characterized by a Mediterranean climate, where mountains act as a reservoir for the dry season and their importance depends on the durability and volume of their snow cover [3]. In addition to fresh water, mountain ecosystems provide various ecosystem services. These include providing raw materials (wood, fiber, and mineral resources), climate and natural hazard regulation, pollination [5,6], and cultural services [7]. The demand for ES is distributed differently worldwide, with ES for food, timber, fodder, and water in high demand in developing countries, whereas in rich countries, the highest demand is for recreational ES [6].

Whatever the topographic configuration in which they are located, forest ecosystems are also a source of the most fundamental ES for humans [8,9]. More than a billion people, including the world’s poorest, directly depend on forests for their livelihood [10,11]. They are home to over 80% of the world’s terrestrial biodiversity [10] and contribute to primary production [12]. Forests are the largest carbon reservoir of all terrestrial ecosystems and store more than 80% of all terrestrial airborne carbon (C) [13], thus contributing to climate change mitigation [10,11,14].

Mediterranean forests, considered a unique global natural heritage in terms of biological diversity [15], particularly contribute to reducing poverty in societies and maintaining local economies. These forests provide a wide range of ES [16]. In terms of provisioning services, they provide wood products (charcoal) and non-wood products (food, medicines, resins, game, etc.). The latter is extensively collected in many Mediterranean forests [17,18]. In terms of regulating services, these forests play a fundamental role in soil and watershed protection and contribute to microclimate improvements. Finally, these ecosystems have a great capacity to provide cultural and tourism ES [19].

Today, forests and mountains are under pressure and their capacity to provide all these goods and services is increasingly affected by environmental and social changes. The demand for ES continues to grow, whereas the supply of mountain ES decreased between 2000 and 2010 worldwide [6]. Even if mountain areas covered by vegetation remain better protected, especially against the risk of erosion, Mediterranean forests have to face various threats of anthropogenic and/or natural origin, such as urban development, fires, and climate change, which cause a rapid change in the water cycle (increased evaporation and decreased precipitation), drought, accelerated desertification, the displacement of biodiversity, and floods. This has an impact on the dynamics of land-use changes [17].

Thus, maintaining the availability and sustainability of ES under current global change scenarios appears to be a challenge and requires integrating ES into management strategies [20,21]. Hence, the aim of this work is to study the potential ES that mountain forests of the Mediterranean basin—one of the richest in endemism [22]-could provide. We particularly target the southern part of the basin in North Africa to contribute to research in this under-documented area. Indeed, the literature has shown that only 9% of the research articles on the multiple ES provided by Mediterranean forests come from the southern and southeastern countries of the Mediterranean basin [19].

In Algeria, forests are rare, their overall area representing nearly 2% of the total national territory and 11% of the country’s north, where there is a climate favorable to the presence of tree cover [23]. In this country, the concept of ES remains little studied, and only a few forest ecosystems have been the subject of studies on ES. We cite, among others, the forests of Chréa and Senalba (2011–2016) [24] and the national park of El Kala studied by Khaznadar M. (2016) [25].

To assess potential ES, i.e., the general capacity of our study area to provide ES, we used botanical data and remote sensing for the delineation of “ES Provisioning Units” (SPU) [26].

Geospatial Data from Remote Sensing for Ecosystem Services Characterization

Earth remote sensing in general and the resulting geospatial data offer multiple potentialities for assessing ecosystem services [27,28]. Indeed, the spatially exhaustive character of all these techniques helps assess a territory continuously and at different scales of analysis. The temporal scale is also apprehended by the repetitiveness of the satellite’s passage, providing time series capable of evaluating the seasonal and/or inter-annual variations in the remotely sensed information [29].

Schematically, there are two main approaches to mapping the supply of ES using remotely sensed data. The oldest and simplest approach uses land cover (LC) categories produced by satellite image classification as indirect support [27,30]. Each thematic land cover class is thus assumed to qualitatively represent a biophysical or other value (economic, social) of one or more ecosystem services [28,31]. Although criticized because of the standardization of ecosystem qualities or conditions simplified by a single cartographic representation [27,31], this approach nevertheless retains the major interest in integrating the concept of ES at the geographical level of the landscape [32]. A landscape assessment of ES is not necessarily limited to land-use/land-cover information (LULC) but allows for the integration of other environmental descriptors, such as soil, climate, surface geomorphology, substrate nature, or hydrological characteristics [32,33], which can be modeled by various geospatial data [28,34]. In this sense, the geographic dimensions of the landscape constitute a Service Provision Area (SPA) [33] or Service Provision Unit (SPU) [26] that has ecological significance at the scale considered. The concept of service providers was conceived as a continuum including different organizational ecological levels from populations of a single species to biomes, passing through communities, habitat types, or landscapes, among others. This multiscalar notion implies questioning the level of the aggregation of the cartographic unit used to spatially assess the supply of ES [33,35]. Does the map unit correspond to the minimum area for estimating a homogeneous and constant supply of ES over the dimensions of this area, or should all these minimum units be aggregated to evaluate the supply of ES over a spatial unit of greater dimension [36]? The answer to this question logically depends on the type of ES to be mapped and the purpose of the assessment. For some ES, the pixel size will be considered relevant, whereas for others, it will be necessary to opt for a more integral mapping unit such as a watershed or an ecoregion [33].

This constraint on the spatial suitability of mapping units is not only restricted to the indirect qualitative approach to estimating ES provision, but it also applies to the quantitative measurement of ES by remote sensing, which is the most credible and currently recommended methodological alternative [27,31]. The variety of sensor types used in satellite imagery, each with its own acquisition characteristics, implies a disparity of spatial resolutions that must necessarily be adapted according to the ES to be quantified. Nevertheless, under this condition of spatial adequacy between the satellite data and the mapping unit of the ES, there is currently a wide range of quantitative remote sensing products that can be directly or indirectly related to ecological and ecosystem processes and, more particularly, to their associated functions [27,31,37]. In particular, Pettorelli et al. (2018) [37] emphasize the key element of ecosystem functions in linking biodiversity, biophysical properties that can be measured by remote sensing, and ecosystem services. From fieldwork in ecology, it has been shown that there is a strong relationship between biodiversity in its taxonomic, structural, and especially functional aspects through the expression of associated traits and the functioning of ecosystems [38,39]. Among the set of data and methods of Earth observation applicable to ES, the transposition of functional and structural diversity in the form of a set of biotic or abiotic traits contained in the electromagnetic signal measured by a remote sensing sensor, called a spectral trait by Lausch et al. (2016) [40], appears as a most relevant concept to quantify some ES. However, the successful application of spectral traits most often requires access to either very high spatial resolution (VHR), imaging spectroscopy, or LIDAR data for which the acquisition costs are not very affordable, and/or experience with the appropriate inversion methods to convert the signal to biophysical/biochemical value [27,30]. These constraints lead ecologists or geographers who are relatively unfamiliar with these complex techniques to prefer free data and simpler methods [41] in line with their objectives and the particularities of their study areas. In the context of our study, given our objective to outline the mapping of ES in the context of Algerian forests, this is the methodological choice we have made (Figure 1).

2. Materials and Methods

2.1. Study Area

The forest of Ouled Hannèche is located in the central part of the Hodna Mountains. The latter, draining the Chott El Hodna (area classified RAMSAR in 2001), form a junction point between the Tellian Atlas in the north and the Saharan Atlas in the south [42]. Administratively, the forest is located south of the wilayas of Bordj Bou Arréridj and Setif (Figure 2). It is located over six municipalities covering an area of 6232 ha [43].

The study area presents a mountainous landscape that is rather rough with an average altitude of about 1500 m. The most important relief is the crest of the Echelendj Jebel (highest point: 1875 m).

The forest is mainly composed of holm oak (Quercus rotundifolia) in a semi-arid bioclimatic environment [44]. It also contains the Atlas Cedar (Cedrus atlantica), an endemic species of North Africa, considered endangered by the International Union for Conservation of Nature (IUCN), and many other species protected in Algeria such as Pistacia atlantica, Juniperus oxycedrus, and Juniperus phoenicea [45].

2.2. Phytoecological Data

The vegetation surveys were carried out by stratified sampling in the Ouled Hannèche forest. Several main criteria were taken into account to select the stations to be sampled: altitude, exposure, topography, and physiognomic homogeneity of the plant cover. Floristic surveys were conducted during different seasons for two years (2018 and 2019). The forest was surveyed most often in the spring and at least once during the other seasons. A total of 59 phytoecological surveys with a minimum area of approximately 400 m² were conducted and geolocated. Each plot surveyed is characterized by a floristic list of species in order of importance of their cover and the dominant species. The dominant species is the most abundant species in the plot. It is the species that contributes the most to the appearance of the vegetation plot by dominating the other species with its physiognomy [46].

In total, 68 surveys were considered in this study (59 carried out by us and 9 carried out by the Bureau National d’Études pour le Développement Rural (BNEDER)).

2.3. LANDSAT Images and Topographic Correction

Three LANDSAT 8 Operational Land Imager (OLI) images were acquired from the USGS database. Their acquisition dates correspond to three distinct seasons in 2017: an image from 17 February for winter, 8 May for spring, and 25 June for summer. The year 2017 was chosen rather than the concurrent 2018–2019 survey dates because it was the only year that guaranteed cloud-free images in all survey seasons. These images are provided pre-processed in surface reflectance at ground level using the LEDAPS procedure [47] and orthorectified. The Ouled Hannèche forest is located in a mountainous region with steep relief; it was deemed necessary to perform a topographic correction of the images’ reflectances. The objective of the topographic correction or illumination correction is to compensate for differences in solar illumination (irradiance) between locations with different slopes and exposures to obtain the values of luminance or reflectance that the sensor would have acquired if the observed surface was perfectly flat [48]. This radiometric normalization procedure is unavoidable when the images are acquired during winter when the low zenithal solar angle can cause serious shadow effects [49]. Among the different categories of topographic correction methods, those that present the best compromise between the complexity of implementation and the quality of the results obtained are the semi-empirical methods [50]. They are based first on the calculation of the cosine of the solar incidence angle representing the illumination conditions for each pixel and calculated from the acquisition of the solar and topographic angles (slope and exposure) estimated from an ALOS AW3D30 Digital Surface Model (DSM) [48]. In our study, we applied the semi-empirical SCS + C method [51]. It combines the modeling of the geometry of the “sun–canopy–sensor” (SCS) complex with a corrective factor specific to each spectral band (c-correction) relating the reflectance to the illumination conditions. In addition, to improve the correction, we chose to vary parameter C according to stratification by slope classes [52]. The topographic correction formula is as follows:

${\rho }_{SCSC,\lambda }={\rho }_{\lambda }\frac{\cos \beta \cos {\theta }_s+c_{\lambda }}{\cos {\gamma }_i+c_{\lambda }}$

where ${\rho }_{SCSC,\lambda }$is the ground reflectance for the λ-band corrected by the SCS + C method,${\rho }_{\lambda }$ is the uncorrected ground reflectance for the $\lambda$ band, $\beta$ is the slope computed from the DSM, ${\theta }_s$ is the solar zenith angle at the date of image acquisition, and $c_{\lambda }$ is the empirical coefficient of the $\lambda$ band, calculated as the ratio of the intercept and slope of the illumination cosine regression ($\cos {\gamma }_i$) by the ground reflectance of each band. The coefficient $c_{\lambda }$is determined individually for each regular 5° slope interval [50,52].

2.4. Supervised Support Vector Machine Classification

The main land cover (LC) types present in the study area were obtained by a supervised classification employing the Support Vector Machine (SVM) algorithm using a Gaussian kernel [53] on the LANDSAT 8 spectral bands calibrated in surface reflectance normalized for the effects of topography. The classification was applied for each date, taking into account a land cover typology considering the physical nature of the material covering the ground. This typology gathers a set of seven classes: GV, which corresponds to photosynthetically active “green” vegetation; NPV, which is assimilated here to senescent or non-photosynthetically active vegetation; SOIL for bare soil; SNOW for snow cover; ROCK for outcropping rock; URBAN for built-up areas; and NA for unclassified artifacts. The training data used to calibrate the model, 237 samples, come either from a photo-interpretation procedure or from field surveys. The classifications for each period are then combined to obtain a final classification that provides information on the succession of physical states in which the land cover is found at different times of the year. We consider that this type of classification is of interest for the mapping of ES because it provides bio-geophysical-based qualitative information on the land cover, which an LC classification defined according to a more traditional nomenclature does not directly indicate. It differs, however, from a spectral unmixing method for which a sub-pixel proportion of the material can be estimated [54].

K-fold cross-validation with 10 folds repeated 10 times was performed to estimate the accuracy of the model (Table 1).

2.5. Classification of Landforms Using the Geomorphon Method

Jasiewicz and Stepinski (2013) [55] developed a method for the classification, mapping, and analysis of landform features from a Digital Terrain Model (DTM) called geomorphons. It makes an important contribution to the knowledge of the physical environment and provides an overview of the landscape structure at different scales. This method is based on the principle of pattern recognition. It relies first on texture analysis of the DTM elevation according to the concept of the local ternary pattern applied to the eight neighboring pixels of a central cell to define the visibility neighbors. These are determined, according to a search distance (L) from the line-of-sight principle along the eight main directions to better fit the local topography [55].

The line of sight is used to estimate the elevation angle. For each direction, a set of elevation angles is calculated. These elevation angles are then used to estimate the zenith angle and the nadir angle.

The zenith angle $\left(D\varphi L\right)$ is defined as the difference between the vertical and the maximum elevation angle in any of the eight directions. The angle at nadir $\left(D\psi L\right)$ is defined as the difference between the vertical and the minimum angle of elevation in a neighborhood direction. In a viewing direction and along a search distance, the value of the difference between each pair of zenith and nadir angles with respect to a flatness threshold (t) defines an element ($D\mathsf{\Delta }L$) of the local ternary pattern. It provides information about the direction of the variation of the local relief within a certain scan distance and can be summarized as follows:

$D\mathsf{\Delta }L=\left\{\begin{array}{c}1 if D\psi L-D\varphi L>t Raised side towards the central pixel.\\ 0 if\left|D\psi L-D\varphi L\right|<t Plane relief towards to the central pixel.\\ -1 if D\psi L-D\varphi L<-t Lowered side towards the central pixel. \end{array}\right.$

The computation of each feature in the eight viewing directions results in a local ternary pattern named “geomorphon”, which is representative of a landform morphology. For simplification, Jasiewicz and Stepinski (2013) [55] summarized the 10 most common geomorphons of a landscape on the Earth’s surface.

The geomorphon classification was performed on the ALOS AW3D30 DSM [56] at a 30 m spatial resolution from the QGIS software using the GRASS extension “r.geomorphon”. The search distance (L) was set to 48 m and the flatness threshold to 1°.

2.6. Normalized Difference Vegetation Index (NDVI)

The NDVI is a vegetation index reflecting the properties of photosynthetic vegetation [57,58]. Indeed, generally, the vegetation absorbs solar radiation in the blue and red wavelengths due to the presence of photosynthetic pigments in the leaves (chlorophyll, xanthophyll, carotene) and it reflects it in the near-infrared (NIR) wavelength due to the structure of the mesophyll of the leaf [59].

Its formula is:

$\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$

where ${\rho }_{Red}$ and ${\rho }_{NIR}$are, respectively, the reflectances in the red and near-infrared wavelengths. In our study, topography-corrected ground reflectances were employed.

It is a simple vegetation index to calculate and interpret, ranging from −1 to +1 and gradually indicating the degree of “greenness” provided by plant activity [60]. Negative values generally indicate the presence of water, whereas those close to 0 indicate the absence of vegetation. It has the advantage in that it can be computed from a wide range of passive optical sensors, including older ones, which gives it the ability to provide information on past ecosystem conditions [60].

Despite some shortcomings due to saturation effects or the influence of soil or other noisy sources on the signal, NDVI is positively correlated with biophysical variables such as the net primary production (NPP), leaf area index (LAI), and fraction of photosynthetically active radiation (FAPAR) [60]. These general capabilities make it an integrative indicator of plant ecosystem functions and justify its frequent use in mapping different categories of ES, as demonstrated, for example, by Krishnaswamy et al. (2009) [61] with a multidate distance of NDVI in relation to biodiversity and carbon storage services, or Paruelo et al. (2016) [62] and Jullian et al. (2021) [63] with the Ecosystem Services Provision Index (ESPI) synthesizing several ES at different spatial scales.

2.7. Rao’s Quadratic Entropy

In addition to biodiversity measurements made in the field, the estimation of landscape biodiversity from satellite images is a credible and relevant approach [64,65,66]. It is based on the Spectral Variation Hypothesis [67], which states that the higher the environmental heterogeneity of a territory, the higher the species diversity. On the scale of satellite images, this environmental heterogeneity, indicative of the various ecological niches, is revealed by the spectral variability of reflectance values and/or indices. The most frequently used measure to quantify this spectral diversity is the Shannon entropy [68]. However, Rocchini et al. (2017) [69] demonstrate that the Shannon entropy has serious flaws in that it only considers the relative abundance of reflectances, neglecting the magnitude of the differences in values between these same reflectances. They propose to use Rao’s Q index instead [70,71], which explicitly considers the distance between the numerical values of the pixels. This method is implemented in R from the package rasterdiv and also offers the possibility to compute the Rao’s Q index of multidimensional data such as the several spectral bands or vegetation indices acquired at different dates. In this work, we applied this index to the NDVI images of 8 May and 25 June 2017, previously masked to keep only the forest surfaces. The NDVI values are calibrated in integer values to avoid the error where the diversity measure is only an artifact due to the length of the significant numerical cipher of the original NDVI kept in the floating type [72]. The distance chosen is the Euclidean distance, and the convolution window chosen is of dimension 3 × 3 pixels.

Its formula is the following:

$Q_{rs}={\displaystyle \sum }_{i=1}^{F-1}{\displaystyle \sum }_{j=i+1}^F{d}_{ij}\times p_i\times p_j$

where $Q_{rs}$ is Rao’s quadratic entropy applied to remote sensing images. $p$ is the relative abundance of the pixel value in the focal window (F). $d_{ij}$ is the spectral distance between the ith and jth pixel values ($d_{ij}=d_{ji}$ and $d_{ii}=0$). i denotes pixel i, j denotes pixel j.

2.8. Topographic Landscape Diversity: Shannon Entropy of Geomorphon Classes

From the morphological units “geomorphons” previously created (Section 2.5), we estimated the topographic landscape diversity on the territory of the Ouled Hannèche forest. It was chosen to apply the Shannon entropy [69] to estimate the topographic environmental heterogeneity of landscape elements [73]. Doxa and Prastacos, (2020) [74] emphasized the advantage provided by Rao entropy over Shannon entropy to characterize the landscape heterogeneity from continuous variables, but due to the categorical nature of the geomorphon classes, the Shannon index remains more relevant in our case. The higher the value of this index, the more it reflects the diversity of morphological forms within a landscape of variable geographical extent. In this work, the geographic dimension of the landscape is reduced to a 9 × 9 pixel convolution window, where one pixel is at a resolution of 30 m on a side. Within this window, the Shannon index is calculated according to the following formula:

$H=-{\displaystyle \sum }_i^n{p}_i\times \ln \left(p_i\right)$

where $H$ denotes the Shannon entropy, $n$ is the total number of individuals present in the window under consideration, and $p$ is the proportion of individuals belonging to the geomorphon category $i$ to the total number of individuals.

2.9. Potential Ecosystem Services Favorability Map

The provision of potential ecosystem services is a synthetic decision combining different factors whose relative importance is prioritized according to a value judgment involving the evaluators. From a geospatial perspective, one of the approaches to perform this synthesis is to use a multi-criteria analysis method incorporated in a GIS forming a system called Spatial Multicriteria Analysis (SMCA) or Spatial Multi-Criteria Decision Support System (MC-SDSS) [75]. It was therefore decided to adopt such a system to map the potential ES of the Ouled Hannèche forest by combining the three quantitative indicators described above. Specifically, these three factors are the average NDVI between 8 May 2017 and 25 June 2017, the multidimensional Rao entropy between these same dates, and the Shannon entropy of the geomorphon classes. All three are masked to keep only the pixels located on the land cover classes corresponding to the GV-GV-GV, GV-NPV-NPV, SNOW-GV-GV, and NPV-GV-GV classes.

The MC-SDSS procedure implies that the indicators are standardized to be established in a common evaluation scale [75]. A linear normalization according to the rank method (min–max) was applied to the three rasters using the spatialEco package [76]. This transformation is applied to data without negative values. The formula is the following:

$v\left(a_{ik}\right)=\left(\left(a_{ik}-mi{n}_i\left\{a_{ik}\right\}\right)/\left(ma{x}_i\left\{a_{ik}\right\}-mi{n}_i\left\{a_{ik}\right\}\right)\right)$

where $v\left(a_{ik}\right)$ is the value for pixel $i$ of the $k\mathrm{th}$ criterion,$a$, normalized to a score from 0 to 1. $mi{n}_i\left\{a_{ik}\right\}$ and $ma{x}_i\left\{a_{ik}\right\}$ are the minimum and maximum values of the $k\mathrm{th}$ indicator, respectively.

The decision rule combining the normalized raster criteria is then performed according to the linear weighted combination (WLC), assigning a relative importance weight according to each indicator [75,77]. The evaluation of the importance of the criteria weights is a primordial step for which we used the analytic hierarchy process (AHP) of Saaty (1980, 2004) [78,79]. It should be noted at this point that we decided to apply a single weight per criterion, which corresponds to the global approach of weighting, implying a certain spatial homogeneity of the criterion. Otherwise, a local weighting adjusted spatially within the criterion can be applied. Saaty’s pairwise comparison uses a scale of 1 to 9 to assess the preferences of one indicator over another per pair of criteria. The relative importance of each factor in the assessment of ES is thus transcribed into a numerical value scale and organized in matrix form (Table 2). From this matrix, a vector containing the criteria weights is obtained by calculating its largest eigenvalue. The weights are then normalized so their sum is unitary (Table 3). The coherence of the judgement deciding the attribution of the weights is estimated thanks to a ratio whose value must not be higher than 0.1 (Table 4). These steps were performed using the ahpsurvey package of R [80].

Finally, the WLC decision rule is translated as an algebraic combination of maps by the relation $V\left(A_i\right)={{\displaystyle \sum }}_k^1{w}_{k}v\left(a_{ik}\right)$. $V\left(A_i\right)$ denotes the overall evaluated value associated with each pixel i; $w_k$ is the weight relative to the kth attribute; and $v\left(a_{ik}\right)$ is the value at the location of the same kth criterion.

Our decision synthesis classifies each indicator according to the degree of importance in the global estimation of the ES supply. This relative degree of importance, summarized in Table 2, is justified as follows. We consider that the average NDVI over the dates in May and June is the most important factor in characterizing the integrity of the forest ecosystem, which gives it an essential role in the delivery of provisioning, regulating, and supporting services. Rao’s quadratic entropy applied to multidate NDVI is ranked second because of the information on biotic diversity that it can reveal, synonymous with the variability in the vegetation-related functional traits and thus multifunctionality in terms of ES. Finally, we ranked the Shannon entropy of geomorphon classes as the least important factor because of its indirect character in the provision of ES, although it is characteristic of landscape structures favoring ES of a distinct nature.

3. Results

3.1. Land Cover

The land cover map represents an appropriate reference for the assessment of ES. Each class is considered to be a unit of ES provision [32]. The land cover map of the Ouled Hannèche Mountains is presented in Figure 3.

Dense woody formations, such as forests and dense or tall matorrals, are represented by the combinations containing the GV class, at least for the dates in May and June. The GV-NPV-NPV combination preferentially represents open formations composed of dry-grass-based vegetation or sparse matorrals with low foliage density. The combinations indicating the presence of snow (SNOW-ROCK-ROCK, SNOW-NPV-NPV, and SNOW-GV-GV) show locations favorable to infiltration that could recharge the groundwater table. In this particular case, this snow class characteristic is indicative of an ability to provide water cycle regulation ES. In general, the characteristics of each class allow for the extrapolation of ES in relation to the physical nature of each combination.

Distribution of field surveys describing plant formation types:

ERM (Erme), ACDR_F (Cedrus atlantica forest), ALP_F (Pinus halepensis forest), STP (steppe formation), ACR_MT (matorral of Acer monspessulanum), LW_MT (low matorral), CST_MT (Cistus albidus matorral), SP_OAK_MT (Quercus rotundifolia open matorral), DS_OAK_MT (Quercus rotundifolia dense matorral), HG_OAK_MT (Quercus rotundifolia high matorral), PHL_MT (Phillyrea angustifolia matorral), JP_MT (Juniperus phoenicea matorral), PIS_MT (Pistacia lentiscus Matorral), GRS (grassland).

From our field surveys, 106 plant species belonging to 37 families are identified. The main species, ranked according to an R-ratio corresponding to the number of occurrences of each species in each survey relative to the total number of sites observed, are summarized in Table 5.

3.2. Geomorphons

The comparison of the vegetation surveys to data from the r.geomorphon approach (Figure 4) allowed us to distinguish the following pattern: the summits are mainly devoid of vegetation; the ridges (in the northern part) and the upper slopes are populated by cedar forests (Cedrus atlantica) with an undergrowth of holm oak (Quercus rotundifolia), juniper (Juniperus oxycedrus), and a rich and diverse herbaceous layer. The spurs are often occupied by steppe formations based on Stipa tenacissima and Ampelodesmos mauritanicus with the presence of juniper (Juniperus phoenicea). The slopes are the best represented; they are often populated by dense formations, such as matorrals of holm oak (Quercus rotundifolia) and maple of Montpellier (Acer monspessulanum) with the presence of Crataegus monogyna and Rosa canina. In the southern part of the massif, the slopes are covered with pine forests (Pinus halepensis) and matorrals of white rockrose (Cistus albidus), pistachio tree (Pistacia lentiscus), Aleppo pine, Phillyrea angustifolia, and Olea europea, with the presence of Jasminum fruticans. The hollows on the slopes are often populated by grasslands. In the lower slopes, we find low matorrals of Calycotome spinosa and matorrals based on Rosmarinus officinalis and Globularia alypum with the presence of Juniperus phoenicea and Ampelodesmos mauritanicus. Finally, the valleys are dedicated to agriculture (cereal growing, arboriculture, market gardening).

3.3. Normalized Difference Vegetation Index (NDVI)

In general, we find moderate NDVI values (between 0.2 and 0.5) throughout the Ouled Hannèche forest for all three dates (Figure 5). Excluding places covered with snow during winter, the map of the seasonal variation of NDVI indicates that some territories in the forest have a high value of NDVI of above 0.6 for the three periods considered. These are the dense forest formations. This underlines the capacity of these territories to maintain ES related to active vegetation and the amount of aboveground biomass. These services are, for example, the climate-regulation service resulting from carbon storage and sequestration and the erosion-regulation service, since the NDVI is a useful indicator for predicting the risk of erosion and quantifying the state of soil fertility [30].

The NDVI values are distributed according to the type of plant formation and its floristic composition. They were extracted from the different existing plant formations based on our field surveys. The results are shown in graphical form in Figure 6.

The formations with a high NDVI are the Cedrus atlantica forests, with a mean value for the surveys corresponding to this species close to 0.6 for the periods in May and June. Acer monspessulanum matorrals have a relatively constant average NDVI value in each period, close to 0.55. The NDVI of Pinus halepensis forests and Quercus rotundifolia matorrals has an average value of 0.5. Open and more or less clear formations growing on rocky ground, such as steppe formations and Juniperus phoenicea matorral, have among the lowest NDVI values (between 0.2 and 0.3 on average). The Erme is made up of geophytes, i.e., plants whose vegetative buds spend the unfavorable period (winter here) in dormancy buried in the ground and whose aerial parts only appear at the beginning of the favorable period in spring [81]. This explains their negative NDVI values for the date of 17 February 2017.

3.4. Rao’s Quadratic Entropy

Our results presented in Figure 7 indicate the level of spectral diversity within the Ouled Hannèche forest. In general, the two monodate indices, i.e., the one from 8 May 2017 and the one from 25 June 2017, show the same geographical distributions of biodiversity. We note low values of the Rao index in the southeastern part of the forest, which represents a certain spectral homogeneity that reveals a low plant diversity. This is consistent with the species present being dominated by the Oleaceae and Anacardiaceae families, which are characterized by structural uniformity and low canopy density. The locations with the highest Rao entropy values are preferentially distributed in the north and west of the forest area. This is consistent with the vegetation organization at these same sites where we found alternating cedar (Cedrus atlantica) forests with grasslands and matorrals of Quercus rotundifolia, Acer monspessulanum, and Juniperus oxycedrus.

The multidate index calculated from the two dates simultaneously indicates a global heterogeneity in the studied territory, with a dominance of high values reflecting a temporal diversity linked to the consideration of phenology in the calculation of the quadratic entropy of Rao. However, we can see that the geographical distribution of the spectral heterogeneity remains roughly similar to that of the monodate indices, with the most diverse sites being kept at the same locations.

3.5. Shannon Entropy of Geomorphon Classes

The Shannon entropy map of the geomorphon classes indicates a high global variability of the topo-morphological forms in the studied territory (Figure 8). In the southeastern part, we note that the landscape geodiversity, rather high in this location, does not necessarily coincide with the biotic diversity expressed by the Rao entropy. Geomorphology is probably not the determining factor in species distribution in this location, for which the aridity of the climate is probably preponderant. However, in the northern part, the Shannon entropy is likely indicative of the rugged terrain in this area, where a chevron-shaped geomorphic structure predominates. This structure is at the origin of the differences in exposure and solar illumination that determines species diversity according to the level of adaptation to these changing ecological conditions.

3.6. Ecosystem Services Supply Potentiality Map

Logically, the map of the potentiality score for the provision of ES shows a spatial variation representative of the favorability of each indicator and their combination (Figure 9). This score reveals a difference in the general capacity to provide ES according to the ecosystems present in the Ouled Hannèche forest. Its geographical distribution is related to the spatial distribution of the main plant formations identified by our phytoecological surveys and classified according to the dominant species in the survey area.

The highest score is concentrated in the formations based on cedar (Cedrus atlantica). These are true forest formations with trees higher than 7 m and a dense and diverse herbaceous layer. The other sites with high scores correspond to the high matorral of Quercus rotundifolia and the matorral of Acer monspessulanum. The medium scores correspond to less dense formations, such as the matorral with holes of Quercus rotundifolia and the forest with holes of Pinus halepensis, distributed uniformly over large areas of the territory. The sites with low scores are generally formed by low-sparse matorral of Cistus albidus, which usually develops after fires, or clear matorral of Pistacia lentiscus, Phillyrea angustifolia, or steppe formations.

4. Discussion and Conclusions

In a particular biogeographic context, subject to a Mediterranean bioclimate under the intermittent but increasingly pronounced Saharan influence [44], our work highlights the potential provision of ES in an Algerian mountain forest. The limited availability of georeferenced environmental maps in this region of the world, especially at the spatial scale of the Ouled Hannèche forest, led us to create qualitative and quantitative geospatial indicators that can be assimilated to the spatial units of ES provision [28]. Given the multiplicity of ES that can be estimated by remote sensing, it appears that a generic approach, as proposed in our study, is often necessary upstream to better identify the overall capacity of a territory in terms of productivity and diversity. This methodological approach is relatively simple but remains entirely relevant in a logic of complementarity and dependence between the datasets used (see Figure 1).

Each indicator thus provides biotic or abiotic information related to the ES, which can be broken down according to different spatial scales and different types of SPU [32,33]. From a qualitative point of view, the land cover classes derived from the SVM classification provide information at the smallest scale, followed by the geomorphon classes describing the landscape in a more detailed way. Locally, the vegetation formations indicate the ES supply in the geographical space where they are present.

This same logic of spatial nesting is found in the quantitative indicators. The size of the mobile window defines the geographical scale of assessment of the supply of ES and conditions the local estimate. The abiotic environmental heterogeneity [82] represented by the Shannon entropy of geomorphon classes according to a 9 × 9 window (72,900 m² for LANDSAT) reveals the largest spatial landscape scale. The Rao entropy of the 3 × 3 focal window (8100 m² for LANDSAT) is at the intermediate scale, whereas the NDVI incorporates information on the minimum area that is the ground resolution of the sensor (900 m²). The neighborhood relationship of varying sizes reflects variations in the abiotic and biotic components of ecosystems along gradients. For Shannon heterogeneity, it synthesizes the topographic complexity of a mountain range, which is a determining and explanatory factor of biogeographic speciation at the origin of taxonomic and functional diversity [82,83]. Similarly, the measure of dissimilarity between NDVI values described by Rao’s quadratic entropy quantifies the variations in the biotic signal according to different portions of the considered analysis space, which is synonymous with the estimation of biodiversity β [64,84,85].

This biodiversity β reflects the fluctuation of the response of functional traits according to the adaptation of ecosystems to environmental factors and exerts significant control on different ecosystem services through the modulation of traits [39]. Figure 6 illustrates this principle by the variability of NDVI values according to the type of plant formations. Indeed, the spatial adequacy between the survey area (400 m²) and the pixel size means that the local NDVI signal reflects the expression of the functional traits of the dominant species in the survey area. A difference in the NDVI value per formation can therefore be equated to a difference in the integrated signature of functional traits [35,40,86] and thus to a disparity in the level of ES provision. In summary, the variation in NDVI due to plant formation is a surrogate for expressing functional diversity in the Ouled Hannèche forest. Despite using a simple index restricted to the red and near-infrared spectral bands and tested on a relatively limited sample of surveys, these results are interesting for this type of Mediterranean forest.

Indicators other than NDVI could have been used to more effectively estimate ecosystem functions. Variables such as gross primary productivity (GPP), net primary productivity (NPP), leaf area index (LAI), fraction of photosynthetically active radiation (FAPAR), surface temperature (LST), albedo, and evapotranspiration are more directly indicative of major biogeochemical cycles, such as that of carbon and water, which underlie many ES [27,37]. They are classically derived in the form of TERRA/AQUA MODIS time series, offering the possibility of continuous seasonal monitoring of the phenology of ecosystem functioning attributes and even detecting reversible or irreversible disturbances that can affect them [86,87]. However, the relatively coarse spatial resolution of these sensors (between 500 m and 1 km) is not well suited for monitoring the ecosystems of the Ouled Hannèche forest. At the local scale of our study area, our experience has indeed shown us that the spatial resolution of LANDSAT probably corresponds to the upper limit of pixel size for which the measured spectral features could match the dimensions of our field surveys. Moreover, another advantage offered by LANDSAT is its long historical coverage [88]. At the desired spatial resolution (i.e., 30 m on a side), the archive of this Earth observation program is accessible back to the 1980s [88], allowing monitoring and comparisons of forest conditions by transferring our method to earlier dates, which is consistent with the concept of spectral feature variation [40,86].

In perspective, methodological improvements are to be considered, especially first by applying them to our existing input dataset. The leaf area index (LAI), which provides information on canopy structure and is related to many functional processes, can be estimated from empirical relationships with reflectance or vegetation index values from LANDSAT images [89]. However, the empirical parameters need to be calibrated to the specificities of the ecosystem under study, and these estimates sometimes lack precision, especially when the vegetation cover is sparse [89]. An alternative to these approaches consists of inverting a Radiative Transfer Model (RTM), which can better return some of the biophysical and/or biochemical properties of forest ecosystems [90]. Among the set of RTMs available, one of the most popular and technically affordable is the PROSAIL model [91,92], which combines the PROSPECT leaf optical properties model with the SAIL bi-directional canopy reflectance model. It allows us to describe, in the solar spectrum (400–2500 nm) and at the canopy level, the spectral features related to leaf biochemistry, mainly chlorophyll, water, and dry matter content, as well as those related to canopy architecture, such as the leaf area index (LAI), leaf angular distribution, and relative leaf size [91]. In addition, it is also possible to return the FAPAR, a key variable for primary production, and the albedo, an essential element for estimating the surface energy balance used to calculate evapotranspiration. Generally, the inversion of this RTM provides the most satisfactory results when applied to hyperspectral data but it can also be adapted to multispectral sensors [92]. More than for LANDSAT, PROSAIL is frequently inverted from SENTINEL-2 images; an inversion procedure using a regression algorithm derived from machine learning is specifically implemented in the SNAP software [93]. Given the spatial resolution of SENTINEL-2, 10 m, this approach should be seriously considered to study the properties of plant communities in the Ouled Hannèche forest.

Ideally, most of the indicators described above could be retained, either as a replacement for or in addition to the variables created during this work, to carry out a decisional synthesis of the potential of the forest’s provision of ES according to a method similar to that indicated in Section 2.9. Advances in the field of spatial multicriteria analysis (SMCA) [94] allow for a better representation of the expert opinions and trade-offs necessary to evaluate the multifunctionality of the forest and the services it provides [95]. On the other hand, Bayesian networks, because of their ability to integrate both qualitative and quantitative data and their explicit consideration of uncertainty, are a powerful tool for spatial modeling of ES [96]. In particular, they can be used to map human demand for ES, defining in combination with the spatialization of supply, the global and complete approach to the socio-ecological concept of ES [28]. However, in some cases, the synthetic assessment of the potential of supply and/or demand of ES may be difficult because of the different physical meanings of the indicators, and it is then necessary to focus on one or a few easily identifiable and interpretable ES. In this case, some biophysical indicators that can be estimated by remote sensing are quite appropriate.

In any case, as emphasized by Andrew et al. (2014) [27], the approach to the assessment and mapping of ES is multidisciplinary in nature and involves collaboration between researchers from different disciplines such as ecology, remote sensing, GIS, geography, and social sciences. In this configuration, we have initiated a mapping of the ES provision of a mountain forest in Algeria, a country for which these kinds of works are rare and for which any attempt to provide useful data for forestry management of the territory is valuable. We hope that this work will encourage further assessment of forest ES in Algeria, allowing us to implement some of the methodological perspectives advocated.

Author Contributions

Conceptualization, A.L. and F.M.; methodology, A.L. and F.M.; software, F.M. and A.L.; formal analysis, F.M.; investigation, A.L.; writing—original draft preparation, A.L. and F.M.; writing—review and editing, F.M. and A.L.; project administration, Y.D. and R.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

The authors thank the foresters of the wilaya of Bordj Bou Arréridj, in particular those of Bordj Ghedir and Ras El Oued for their help in the field.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. General scheme summarizing the method.

Figure 2. Location map of the Ouled Hannèche forest.

Figure 3. Land cover map of the Ouled Hannèche Mountains (coordinates are given in meters in the UTM31N projection system, DATUM WGS84).

Figure 4. Map of geomorphon classes in the Ouled Hannèche forest.

Figure 5. Seasonal changes in NDVI in the Ouled Hannèche Mountains in 2017.

Figure 6. Average NDVI by plant formation type and for each image date.

Figure 7. Maps of Rao’s quadratic entropy. (A,B) monodate indices, (C) multidate indice.

Figure 8. Shannon entropy map of geomorphon classes.

Figure 9. Favorability map of the potential ecosystem services supply of the Ouled Hannèche forest.

Table 1

Summary of the parameters and evaluation of the classification model.

	17 February	8 May	25 June	3 Dates
Model	SVM radial basis	SVM radial basis	SVM radial basis	SVM radial basis
Number of training samples	237	237	237	237
Number of predictors	7	7	7	21
Number of classes	7	6	6	10
Model evaluation	repeated K-fold CV	repeated K-fold CV	repeated K-fold CV	repeated K-fold CV
Number K-fold	10	10	10	10
Number repetition	10	10	10	10
Best sigma parameter	0.2	1	0.5	0.5
Best C parameter	0.1	1.25	15	0.9
Overall accuracy	0.66	0.76	0.78	0.63
Kappa	0.53	0.68	0.68	0.55

Table 2

Saaty AHP matrix.

	Diversity Geomorphon	Rao Entropy	Mean NDVI
Diversity geomorphon	1	0.5	0.25
Rao entropy	2	1	0.33
Mean NDVI	4	3	1

Table 3

Weights assigned.

Diversity Geomorphon	Rao Entropy	Mean NDVI
0.136	0.238	0.625

Table 4

Consistency ratio.

CR

0.017

Table 5

Main species identified in the forest of Ouled Hannèche.

Plant Species	R	Plant Species	R	Plant Species	R
Quercus ilex	42/68	Rosa canina	15/68	Pinus halepensis	9/68
Juniperus oxycedrus	29/68	Globularia alypum	14/68	Rhamnus lycioides	9/68
Juniperus phoenicea	27/68	Cedrus atlantica	13/68	Pistacia atlantica	8/68
Cistus albidus	26/68	Genista microcephala	12/68	Thapsia garganica	7/68
Pistacia lentiscus	20/68	Astragalus sp.	12/68	Lonicera impelaxa	7/68
Ampelodesmos mauritanicus	19/68	Jasminum fruticans	11/68	Teucrium polium	7/68
Phillyrea angustifolia	18/68	Crataegus monogyna	11/68	Bunium sp.	6/68
Santolina rosmarinifolia	17/68	Rosmarinus officinalis	10/68	Bupleurum sp.	6/68
Stipa tenacissima	16/68	Galium sp.	9/68	Eryngium sp.	6/68
Acer monspessulanum	16/68	Crataegus azarolus	9/68	Lomelosia stellata	6/68
Calicotome spinosa	15/68	Asparagus acutifolius	9/68	Aegilops geniculata	6/68

References

1. Mooney, H.A.; Cropper, A.; Reid, W. The Millennium Ecosystem Assessment: What Is It All About?. Trends Ecol. Evol.; 2004; 19, pp. 221-224. [DOI: https://dx.doi.org/10.1016/j.tree.2004.03.005] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16701257]

2. Mainka, S.; McNeely, J.; Jackson, B.; McNeely, J.A. Depend on Nature: Ecosystem Services Supporting Human Livelihoods; IUCN: Grand, Switzerland, 2005.

3. Viviroli, D.; Weingartner, R.; Messerli, B. Assessing the Hydrological Significance of the World’s Mountains. Mt. Res. Dev.; 2003; 23, pp. 32-40. [DOI: https://dx.doi.org/10.1659/0276-4741(2003)023[0032:ATHSOT]2.0.CO;2]

4. Grêt-Regamey, A.; Brunner, S.H.; Kienast, F. Mountain Ecosystem Services: Who Cares?. Mt. Res. Dev.; 2012; 32, (Suppl. 1), pp. S23-S34. [DOI: https://dx.doi.org/10.1659/MRD-JOURNAL-D-10-00115.S1]

5. Payne, D.; Spehn, E.M.; Snethlage, M.; Fischer, M. Opportunities for Research on Mountain Biodiversity under Global Change. Curr. Opin. Environ. Sustain.; 2017; 29, pp. 40-47. [DOI: https://dx.doi.org/10.1016/j.cosust.2017.11.001]

6. Grêt-Regamey, A.; Weibel, B. Global Assessment of Mountain Ecosystem Services Using Earth Observation Data. Ecosyst. Serv.; 2020; 46, 101213. [DOI: https://dx.doi.org/10.1016/j.ecoser.2020.101213]

7. Schirpke, U.; Scolozzi, R.; Dean, G.; Haller, A.; Jäger, H.; Kister, J.; Kovács, B.; Sarmiento, F.O.; Sattler, B.; Schleyer, C. Cultural Ecosystem Services in Mountain Regions: Conceptualising Conflicts among Users and Limitations of Use. Ecosyst. Serv.; 2020; 46, 101210. [DOI: https://dx.doi.org/10.1016/j.ecoser.2020.101210]

8. Cudlín, P.; Seják, J.; Pokorný, J.; Albrechtová, J.; Bastian, O.; Marek, M. Chapter 24—Forest Ecosystem Services Under Climate Change and Air Pollution. Developments in Environmental Science; Matyssek, R.; Clarke, N.; Cudlin, P.; Mikkelsen, T.N.; Tuovinen, J.-P.; Wieser, G.; Paoletti, E. Climate Change, Air Pollution and Global Challenges Elsevier: Amsterdam, The Netherlands, 2013; Volume 13, pp. 521-546. [DOI: https://dx.doi.org/10.1016/B978-0-08-098349-3.00024-4]

9. Masiero, M.; Pettenella, D.; Boscolo, M.; Kanti Barua, S.; Animon, I.; Matta, R. Valuing Forest Ecosystem Services: A Training Manual for Planners and Project Developers; Forestry Working Paper; Food and Agriculture Organization of The United Nations: Rome, Italy, 2019; 216.

10. Ten Brink, P. The Economics of Ecosystems and Biodiversity in National and International Policy Making; Taylor & Francis: Abingdon, UK, 2011; 390.

11. Évaluation des Ressources Forestières Mondiales 2020; FAO: Rome, Italy, 2020; [DOI: https://dx.doi.org/10.4060/ca8753fr]

12. De Groot, R.S.; Wilson, M.A.; Boumans, R.M.J. A Typology for the Classification, Description and Valuation of Ecosystem Functions, Goods and Services. Ecol. Econ.; 2002; 41, pp. 393-408. [DOI: https://dx.doi.org/10.1016/S0921-8009(02)00089-7]

13. Dixon, R.K.; Solomon, A.M.; Brown, S.; Houghton, R.A.; Trexier, M.C.; Wisniewski, J. Carbon Pools and Flux of Global Forest Ecosystems. Science; 1994; 263, pp. 185-190. [DOI: https://dx.doi.org/10.1126/science.263.5144.185]

14. Stocker, T. Climate Change 2013: The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2014.

15. Palahí, M.; Birot, Y.; Borges, J.G.; Bravo, F.; Pettenella, D.; Sabir, M.; Hassen, H.D.; Shater, Z.; Başkent, E.Z.; Kazana, V. A Mediterranean Forest Research Agenda—MFRA 2010–2020. European Forest Institute. Mediterranean Regional Forest Office–EFIMED. 2009; Available online: https://efi.int/sites/default/files/images/publications/Mediterranean%20Forest%20Research%20Agenda%202010-2020_.pdf (accessed on 7 February 2022).

16. Blondel, J. La Production Durable de Biens et Services en Forêt Méditerranéenne: Le Point de Vue de L’écologue. Méditerranéenne; 2009; 30, pp. 133-138. Available online: http://www.foret-mediterraneenne.org/fr/catalogue/id-1178-la-production-durable-de-biens-et-services-en-foret-mediterraneenne-le-point-de-vue-de-l-ecologue (accessed on 2 April 2019).

17. Anaya-Romero, M.; Muñoz-Rojas, M.; Ibáñez, B.; Marañón, T. Evaluation of Forest Ecosystem Services in Mediterranean Areas. A Regional Case Study in South Spain. Ecosyst. Serv.; 2016; 20, pp. 82-90. [DOI: https://dx.doi.org/10.1016/j.ecoser.2016.07.002]

18. Lovrić, M.; Da Re, R.; Vidale, E.; Prokofieva, I.; Wong, J.; Pettenella, D.; Verkerk, P.J.; Mavsar, R. Collection and Consumption of Non-Wood Forest Products in Europe. For. Int. J. For. Res.; 2021; 94, pp. 757-770. [DOI: https://dx.doi.org/10.1093/forestry/cpab018]

19. Nocentini, S.; Travaglini, D.; Muys, B. Managing Mediterranean Forests for Multiple Ecosystem Services: Research Progress and Knowledge Gaps. Curr. For. Rep.; 2022; 8, pp. 229-256. [DOI: https://dx.doi.org/10.1007/s40725-022-00167-w]

20. Benayas, J.M.R.; Newton, A.C.; Diaz, A.; Bullock, J.M. Enhancement of Biodiversity and Ecosystem Services by Ecological Restoration: A Meta-Analysis. Science; 2009; 325, pp. 1121-1124. [DOI: https://dx.doi.org/10.1126/science.1172460] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19644076]

21. Assessment, M.E. Ecosystems and Human Well-Being: Wetlands and Water Synthesis; World Resources Institute: Washington, DC, USA, 2005.

22. Médail, F.; Quezél, P. Biodiversity Hotspots in the Mediterranean Basin: Setting Global Conservation Priorities. Conserv. Biol.; 1999; 13, pp. 1510-1513. [DOI: https://dx.doi.org/10.1046/j.1523-1739.1999.98467.x]

23. Tatar, H. Production forestière, exploitation et valorisation en Algérie. Méditerranéenne; 2012; 33, pp. 361-368.

24. Maire, M.; Garavaglia, V.; Picard, N. Maximize the Production of Goods and Services of Mediterranean Forest Ecosystems in the Context of Global Changes; FAO/Plan Bleu: Rome, Italy, 2016; Available online: https://www.fao.org/3/i5887e/i5887e.pdf (accessed on 7 February 2022).

25. Khaznadar, M. Etude et Cartographie des Services Ecosystémiques du Parc National d’El Kala (Algérie). Ph.D. Thesis; Ferhat Abbes Sétif 1 University: Sétif, Algeria, 2016.

26. Luck, G.W.; Harrington, R.; Harrison, P.A.; Kremen, C.; Berry, P.M.; Bugter, R.; Dawson, T.P.; de Bello, F.; Díaz, S.; Feld, C.K. et al. Quantifying the Contribution of Organisms to the Provision of Ecosystem Services. BioScience; 2009; 59, pp. 223-235. [DOI: https://dx.doi.org/10.1525/bio.2009.59.3.7]

27. Andrew, M.E.; Wulder, M.A.; Nelson, T.A. Potential Contributions of Remote Sensing to Ecosystem Service Assessments. Prog. Phys. Geogr. Earth Environ.; 2014; 38, pp. 328-353. [DOI: https://dx.doi.org/10.1177/0309133314528942]

28. Andrew, M.E.; Wulder, M.A.; Nelson, T.A.; Coops, N.C. Spatial Data, Analysis Approaches, and Information Needs for Spatial Ecosystem Service Assessments: A Review. GIScience Remote Sens.; 2015; 52, pp. 344-373. [DOI: https://dx.doi.org/10.1080/15481603.2015.1033809]

29. Skidmore, A.K.; Pettorelli, N.; Coops, N.C.; Geller, G.N.; Hansen, M.; Lucas, R.; Mücher, C.A.; O’Connor, B.; Paganini, M.; Pereira, H.M. et al. Environmental Science: Agree on Biodiversity Metrics to Track from Space. Nature; 2015; 523, pp. 403-405. [DOI: https://dx.doi.org/10.1038/523403a]

30. Ayanu, Y.Z.; Conrad, C.; Nauss, T.; Wegmann, M.; Koellner, T. Quantifying and Mapping Ecosystem Services Supplies and Demands: A Review of Remote Sensing Applications. Environ. Sci. Technol.; 2012; 46, pp. 8529-8541. [DOI: https://dx.doi.org/10.1021/es300157u]

31. Cord, A.F.; Brauman, K.A.; Chaplin-Kramer, R.; Huth, A.; Ziv, G.; Seppelt, R. Priorities to Advance Monitoring of Ecosystem Services Using Earth Observation. Trends Ecol. Evol.; 2017; 32, pp. 416-428. [DOI: https://dx.doi.org/10.1016/j.tree.2017.03.003]

32. Burkhard, B.; Kroll, F.; Müller, F.; Windhorst, W. Landscapes’ Capacities to Provide Ecosystem Services—A Concept for Land-Cover Based Assessments. Landsc. Online; 2009; 15, 22. [DOI: https://dx.doi.org/10.3097/LO.200915]

33. Syrbe, R.-U.; Walz, U. Spatial Indicators for the Assessment of Ecosystem Services: Providing, Benefiting and Connecting Areas and Landscape Metrics. Ecol. Indic.; 2012; 21, pp. 80-88. [DOI: https://dx.doi.org/10.1016/j.ecolind.2012.02.013]

34. Lausch, A.; Schaepman, M.E.; Skidmore, A.K.; Catana, E.; Bannehr, L.; Bastian, O.; Borg, E.; Bumberger, J.; Dietrich, P.; Glässer, C. et al. Remote Sensing of Geomorphodiversity Linked to Biodiversity—Part III: Traits, Processes and Remote Sensing Characteristics. Remote Sens.; 2022; 14, 2279. [DOI: https://dx.doi.org/10.3390/rs14092279]

35. Abelleira Martínez, O.J.; Fremier, A.K.; Günter, S.; Ramos Bendaña, Z.; Vierling, L.; Galbraith, S.M.; Bosque-Pérez, N.A.; Ordoñez, J.C. Scaling up Functional Traits for Ecosystem Services with Remote Sensing: Concepts and Methods. Ecol. Evol.; 2016; 6, pp. 4359-4371. [DOI: https://dx.doi.org/10.1002/ece3.2201] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27386081]

36. Shen, J.; Chen, C.; Wang, Y. What Are the Appropriate Mapping Units for Ecosystem Service Assessments? A Systematic Review. Ecosyst. Health Sustain.; 2021; 7, 1888655. [DOI: https://dx.doi.org/10.1080/20964129.2021.1888655]

37. Pettorelli, N.; Schulte to Bühne, H.; Tulloch, A.; Dubois, G.; Macinnis-Ng, C.; Queirós, A.M.; Keith, D.A.; Wegmann, M.; Schrodt, F.; Stellmes, M. et al. Satellite Remote Sensing of Ecosystem Functions: Opportunities, Challenges and Way Forward. Remote Sens. Ecol. Conserv.; 2018; 4, pp. 71-93. [DOI: https://dx.doi.org/10.1002/rse2.59]

38. Hooper, D.U.; Chapin, F.S., III; Ewel, J.J.; Hector, A.; Inchausti, P.; Lavorel, S.; Lawton, J.H.; Lodge, D.M.; Loreau, M.; Naeem, S. et al. Effects of Biodiversity on Ecosystem Functioning: A Consensus of Current Knowledge. Ecol. Monogr.; 2005; 75, pp. 3-35. [DOI: https://dx.doi.org/10.1890/04-0922]

39. de Bello, F.; Lavorel, S.; Díaz, S.; Harrington, R.; Cornelissen, J.H.C.; Bardgett, R.D.; Berg, M.P.; Cipriotti, P.; Feld, C.K.; Hering, D. et al. Towards an Assessment of Multiple Ecosystem Processes and Services via Functional Traits. Biodivers Conserv.; 2010; 19, pp. 2873-2893. [DOI: https://dx.doi.org/10.1007/s10531-010-9850-9]

40. Lausch, A.; Erasmi, S.; King, D.J.; Magdon, P.; Heurich, M. Understanding Forest Health with Remote Sensing -Part I—A Review of Spectral Traits, Processes and Remote-Sensing Characteristics. Remote Sens.; 2016; 8, 1029. [DOI: https://dx.doi.org/10.3390/rs8121029]

41. Kc, A.; Wagle, N.; Acharya, T.D. Spatiotemporal Analysis of Land Cover and the Effects on Ecosystem Service Values in Rupandehi, Nepal from 2005 to 2020. ISPRS Int. J. Geo-Inf.; 2021; 10, 635. [DOI: https://dx.doi.org/10.3390/ijgi10100635]

42. Bord, J.P. Cartographie de L’utilisation du sol dans l’Est Algérien: Essai de Zonage Agricole. Ph.D. Thesis; Université Paul Valéry-Montpellier III: Montpellier, France, 1981.

43. Bureau National d’Études pour le Développement Rural (BNEDER). Etude D’aménagement et de Développement de la Forêt Domaniale de Ouled Hanneche Dans la Wilaya de Bordj Bou Arréridj; Bureau National d’Études pour le Développement Rural (BNEDER): Chéraga, Algérie, 2010.

44. Louail, A.; Messner, F.; Missaoui, K.; Djellouli, Y.; Gharzouli, R. Evolution of Local Bioclimates Facing Climate Change in Algeria: Mapping by Global Bioclimatic Classification and Clustering Method. Proceedings of the International Society for Ecological Modelling Global Conference; Salzburg, Austria, 30 September–4 October 2019.

45. Official Journal of the Algerian Republic (JORA). N°03. 18-01-2012. Executive Decree No. 12-03 of 10 Safar 1433 Corresponding to January 4, 2012 Establishing the List of Protected Non-Crop Species. 2012; 44.Available online: https://gazettes.africa/archive/dz/2012/dz-government-gazette-dated-2012-01-18-no-3.pdf (accessed on 21 June 2021).

46. Atbib, M. Etude phytoécologique de la réserve biologique de Mehdia (Littoral Atlantique du Maroc). 1ère partie-La végétation hygrophile de la merja Sidi Bou Ghaba. Bull. L’institut Sci.; 1980; 4, 112.

47. Masek, J.G.; Vermote, E.F.; Saleous, N.; Wolfe, R.; Hall, F.G.; Huemmrich, K.F.; Gao, F.; Kutler, J.; Lim, T.K. LEDAPS Landsat Calibration, Reflectance, Atmospheric Correction Preprocessing Code; ORNL DAAC: Oak Ridge, TN, USA, 2012; [DOI: https://dx.doi.org/10.3334/ORNLDAAC/1080]

48. Sola, I.; González-Audícana, M.; Álvarez-Mozos, J. The Added Value of Stratified Topographic Correction of Multispectral Images. Remote Sens.; 2016; 8, 131. [DOI: https://dx.doi.org/10.3390/rs8020131]

49. Hantson, S.; Chuvieco, E. Evaluation of Different Topographic Correction Methods for Landsat Imagery. Int. J. Appl. Earth Obs. Geoinf.; 2011; 13, pp. 691-700. [DOI: https://dx.doi.org/10.1016/j.jag.2011.05.001]

50. Sola, I.; González-Audícana, M.; Álvarez-Mozos, J. Multi-Criteria Evaluation of Topographic Correction Methods. Remote Sens. Environ.; 2016; 184, pp. 247-262. [DOI: https://dx.doi.org/10.1016/j.rse.2016.07.002]

51. Soenen, S.A.; Peddle, D.R.; Coburn, C.A. SCS + C: A Modified Sun-Canopy-Sensor Topographic Correction in Forested Terrain. IEEE Trans. Geosci. Remote Sens.; 2005; 43, pp. 2148-2159. [DOI: https://dx.doi.org/10.1109/TGRS.2005.852480]

52. Vázquez-Jiménez, R.; Romero-Calcerrada, R.; Ramos-Bernal, R.N.; Arrogante-Funes, P.; Novillo, C.J. Topographic Correction to Landsat Imagery through Slope Classification by Applying the SCS + C Method in Mountainous Forest Areas. ISPRS Int. J. Geo-Inf.; 2017; 6, 287. [DOI: https://dx.doi.org/10.3390/ijgi6090287]

53. Mountrakis, G.; Im, J.; Ogole, C. Support Vector Machines in Remote Sensing: A Review. ISPRS J. Photogramm. Remote Sens.; 2011; 66, pp. 247-259. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2010.11.001]

54. Roberts, D.A.; Gardner, M.; Church, R.; Ustin, S.; Scheer, G.; Green, R.O. Mapping Chaparral in the Santa Monica Mountains Using Multiple Endmember Spectral Mixture Models. Remote Sens. Environ.; 1998; 65, pp. 267-279. [DOI: https://dx.doi.org/10.1016/S0034-4257(98)00037-6]

55. Jasiewicz, J.; Stepinski, T.F. Geomorphons—A Pattern Recognition Approach to Classification and Mapping of Landforms. Geomorphology; 2013; 182, pp. 147-156. [DOI: https://dx.doi.org/10.1016/j.geomorph.2012.11.005]

56. Tadono, T.; Nagai, H.; Ishida, H.; Oda, F.; Naito, S.; Minakawa, K.; Iwamoto, H. Generation of the 30 M-mesh global digital surface model by ALOS PRISM. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2016; 41, pp. 157-162.

57. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring Vegetation Systems in the Great Plains with ERTS. Proceedings of the 3rd ERTS Symposium, NASA SP-351; Washington, DC, USA, 10–14 December 1973; pp. 309-317.

58. Tucker, C.J.; Townshend, J.R.G.; Goff, T.E. African Land-Cover Classification Using Satellite Data. Science; 1985; 227, pp. 369-375. [DOI: https://dx.doi.org/10.1126/science.227.4685.369] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17815712]

59. Wegmann, M.; Leutner, B.; Dech, S. Remote Sensing and GIS for Ecologists: Using Open Source Software; Exeter Pelagic Publishing Ltd.: Exeter, UK, 2016; 352.

60. Pettorelli, N. The Normalized Difference Vegetation Index; 1st ed. Oxford University Press: Oxford, UK, 2013; 194.

61. Krishnaswamy, J.; Bawa, K.S.; Ganeshaiah, K.N.; Kiran, M.C. Quantifying and Mapping Biodiversity and Ecosystem Services: Utility of a Multi-Season NDVI Based Mahalanobis Distance Surrogate. Remote Sens. Environ.; 2009; 113, pp. 857-867. [DOI: https://dx.doi.org/10.1016/j.rse.2008.12.011]

62. Paruelo, J.M.; Texeira, M.; Staiano, L.; Mastrángelo, M.; Amdan, L.; Gallego, F. An Integrative Index of Ecosystem Services Provision Based on Remotely Sensed Data. Ecol. Indic.; 2016; 71, pp. 145-154. [DOI: https://dx.doi.org/10.1016/j.ecolind.2016.06.054]

63. Jullian, C.; Nahuelhual, L.; Laterra, P. The Ecosystem Service Provision Index as a Generic Indicator of Ecosystem Service Supply for Monitoring Conservation Targets. Ecol. Indic.; 2021; 129, 107855. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.107855]

64. Rocchini, D.; Salvatori, N.; Beierkuhnlein, C.; Chiarucci, A.; de Boissieu, F.; Förster, M.; Garzon-Lopez, C.X.; Gillespie, T.W.; Hauffe, H.C.; He, K.S. et al. From Local Spectral Species to Global Spectral Communities: A Benchmark for Ecosystem Diversity Estimate by Remote Sensing. Ecol. Inform.; 2021; 61, 101195. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2020.101195]

65. Gillespie, T.W.; Foody, G.M.; Rocchini, D.; Giorgi, A.P.; Saatchi, S. Measuring and Modelling Biodiversity from Space. Prog. Phys. Geogr. Earth Environ.; 2008; 32, pp. 203-221. [DOI: https://dx.doi.org/10.1177/0309133308093606]

66. Boyd, D.S.; Foody, G.M. An Overview of Recent Remote Sensing and GIS Based Research in Ecological Informatics. Ecol. Inform.; 2011; 6, pp. 25-36. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2010.07.007]

67. Palmer, M.W.; Earls, P.G.; Hoagland, B.W.; White, P.S.; Wohlgemuth, T. Quantitative Tools for Perfecting Species Lists. Environmetrics; 2002; 13, pp. 121-137. [DOI: https://dx.doi.org/10.1002/env.516]

68. Shannon, C.E. A Mathematical Theory of Communication. Bell Syst. Tech. J.; 1948; 27, pp. 379-423. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x]

69. Rocchini, D.; Marcantonio, M.; Ricotta, C. Measuring Rao’s Q Diversity Index from Remote Sensing: An Open Source Solution. Ecol. Indic.; 2017; 72, pp. 234-238. [DOI: https://dx.doi.org/10.1016/j.ecolind.2016.07.039]

70. Rao, C.R. Diversity and Dissimilarity Coefficients: A Unified Approach. Theor. Popul. Biol.; 1982; 21, pp. 24-43. [DOI: https://dx.doi.org/10.1016/0040-5809(82)90004-1]

71. Botta-Dukát, Z. Rao’s Quadratic Entropy as a Measure of Functional Diversity Based on Multiple Traits. J. Veg. Sci.; 2005; 16, pp. 533-540. [DOI: https://dx.doi.org/10.1111/j.1654-1103.2005.tb02393.x]

72. Thouverai, E.; Marcantonio, M.; Bacaro, G.; Re, D.D.; Iannacito, M.; Marchetto, E.; Ricotta, C.; Tattoni, C.; Vicario, S.; Rocchini, D. Measuring Diversity from Space: A Global View of the Free and Open Source Rasterdiv R Package under a Coding Perspective. Community Ecol.; 2021; 22, pp. 1-11. [DOI: https://dx.doi.org/10.1007/s42974-021-00042-x]

73. Amatulli, G.; Domisch, S.; Tuanmu, M.-N.; Parmentier, B.; Ranipeta, A.; Malczyk, J.; Jetz, W. A Suite of Global, Cross-Scale Topographic Variables for Environmental and Biodiversity Modeling. Sci Data; 2018; 5, 180040. [DOI: https://dx.doi.org/10.1038/sdata.2018.40]

74. Doxa, A.; Prastacos, P. Using Rao’s Quadratic Entropy to Define Environmental Heterogeneity Priority Areas in the European Mediterranean Biome. Biol. Conserv.; 2020; 241, 108366. [DOI: https://dx.doi.org/10.1016/j.biocon.2019.108366]

75. Malczewski, J. Multicriteria Analysis. Comprehensive Geographic Information Systems; Elsevier: Amsterdam, The Netherlands, 2018; pp. 197-217. [DOI: https://dx.doi.org/10.1016/B978-0-12-409548-9.09698-6]

76. Evans, J.S. Spatialeco. R Package Version 1.3-6. 2021; Available online: https://github.com/jeffreyevans/spatialEco (accessed on 6 May 2022).

77. Malczewski, J. On the Use of Weighted Linear Combination Method in GIS: Common and Best Practice Approaches. Trans. GIS; 2000; 4, pp. 5-22. [DOI: https://dx.doi.org/10.1111/1467-9671.00035]

78. Saaty, T.L. The Analytic Hierarchy Process: Planning, Priority Setting, Resource Allocation; McGraw-Hill International Book Co.: New York, NY, USA, London, UK, 1980.

79. Saaty, T.L. Decision Making—The Analytic Hierarchy and Network Processes (AHP/ANP). J. Syst. Sci. Syst. Eng.; 2004; 13, pp. 1-35. [DOI: https://dx.doi.org/10.1007/s11518-006-0151-5]

80. Cho, H.T.F. Analytic Hierarchy Process for Survey Data in R (Vignettes). 2018; Available online: https://doi.org/10.13140/RG.2.2.25278.13125 (accessed on 6 May 2022).

81. Sirvent, L. Les Types Biologiques: État de L’art, Actualisation des Définitions et Mise en Place d’un Référentiel. 2020; Available online: https://doi.org/10.13140/RG.2.2.14673.89440 (accessed on 6 February 2022).

82. Stein, A.; Kreft, H. Terminology and Quantification of Environmental Heterogeneity in Species-Richness Research. Biol. Rev.; 2014; 90, pp. 815-836. [DOI: https://dx.doi.org/10.1111/brv.12135]

83. Badgley, C.; Smiley, T.M.; Terry, R.; Davis, E.B.; DeSantis, L.R.G.; Fox, D.L.; Hopkins, S.S.B.; Jezkova, T.; Matocq, M.D.; Matzke, N. et al. Biodiversity and Topographic Complexity: Modern and Geohistorical Perspectives. Trends Ecol. Evol.; 2017; 32, pp. 211-226. [DOI: https://dx.doi.org/10.1016/j.tree.2016.12.010]

84. Rocchini, D.; Luque, S.; Pettorelli, N.; Bastin, L.; Doktor, D.; Faedi, N.; Feilhauer, H.; Féret, J.-B.; Foody, G.M.; Gavish, Y. et al. Measuring β-Diversity by Remote Sensing: A Challenge for Biodiversity Monitoring. Methods Ecol. Evol.; 2018; 9, pp. 1787-1798. [DOI: https://dx.doi.org/10.1111/2041-210X.12941]

85. Khare, S.; Latifi, H.; Rossi, S. Forest Beta-Diversity Analysis by Remote Sensing: How Scale and Sensors Affect the Rao’s Q Index. Ecol. Indic.; 2019; 106, 105520. [DOI: https://dx.doi.org/10.1016/j.ecolind.2019.105520]

86. Lausch, A.; Bannehr, L.; Beckmann, M.; Boehm, C.; Feilhauer, H.; Hacker, J.M.; Heurich, M.; Jung, A.; Klenke, R.; Neumann, C. et al. Linking Earth Observation and Taxonomic, Structural and Functional Biodiversity: Local to Ecosystem Perspectives. Ecol. Indic.; 2016; 70, pp. 317-339. [DOI: https://dx.doi.org/10.1016/j.ecolind.2016.06.022]

87. Regos, A.; Gonçalves, J.; Arenas-Castro, S.; Alcaraz-Segura, D.; Guisan, A.; Honrado, J.P. Mainstreaming Remotely Sensed Ecosystem Functioning in Ecological Niche Models. Remote Sens. Ecol. Conserv.; 2022; [DOI: https://dx.doi.org/10.1002/rse2.255]

88. Wulder, M.A.; Coops, N.C. Satellites: Make Earth Observations Open Access. Nature; 2014; 513, pp. 30-31. [DOI: https://dx.doi.org/10.1038/513030a] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25186885]

89. Baret, F.; Guyot, G. Potentials and Limits of Vegetation Indices for LAI and APAR Assessment. Remote Sens. Environ.; 1991; 35, pp. 161-173. [DOI: https://dx.doi.org/10.1016/0034-4257(91)90009-U]

90. Baret, F.; Buis, S. Estimating Canopy Characteristics from Remote Sensing Observations: Review of Methods and Associated Problems. Advances in Land Remote Sensing: System, Modeling, Inversion and Application; Liang, S. Springer: Dordrecht, The Netherlands, 2008; pp. 173-201. [DOI: https://dx.doi.org/10.1007/978-1-4020-6450-0_7]

91. Jacquemoud, S.; Verhoef, W.; Baret, F.; Bacour, C.; Zarco-Tejada, P.J.; Asner, G.P.; François, C.; Ustin, S.L. PROSPECT+SAIL Models: A Review of Use for Vegetation Characterization. Remote Sens. Environ.; 2009; 113, pp. S56-S66. [DOI: https://dx.doi.org/10.1016/j.rse.2008.01.026]

92. Berger, K.; Atzberger, C.; Danner, M.; D’Urso, G.; Mauser, W.; Vuolo, F.; Hank, T. Evaluation of the PROSAIL Model Capabilities for Future Hyperspectral Model Environments: A Review Study. Remote Sens.; 2018; 10, 85. [DOI: https://dx.doi.org/10.3390/rs10010085]

93. Weiss, M.; Baret, F. S2toolbox Level 2 Products: Lai, Fapar, Fcover; Institut National de la Recherche Agronomique (INRA): Avignon, France, 2016.

94. Malczewski, J.; Jankowski, P. Emerging Trends and Research Frontiers in Spatial Multicriteria Analysis. Int. J. Geogr. Inf. Sci.; 2020; 34, pp. 1257-1282. [DOI: https://dx.doi.org/10.1080/13658816.2020.1712403]

95. Baskent, E.Z.; Borges, J.G.; Kašpar, J.; Tahri, M. A Design for Addressing Multiple Ecosystem Services in Forest Management Planning. Forests; 2020; 11, 1108. [DOI: https://dx.doi.org/10.3390/f11101108]

96. Landuyt, D.; Broekx, S.; D’hondt, R.; Engelen, G.; Aertsens, J.; Goethals, P.L.M. A Review of Bayesian Belief Networks in Ecosystem Service Modelling. Environ. Model. Softw.; 2013; 46, pp. 1-11. [DOI: https://dx.doi.org/10.1016/j.envsoft.2013.03.011]