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

: Regardless of their biogeographic origins or degree of artiﬁcialization, 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 speciﬁc 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 difﬁcult. Our work presents an introductory approach, allowing us to set up methodological and scientiﬁc 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. Speciﬁcally, 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 classiﬁcation of topographically corrected LANDSAT images and a geomorphic-type classiﬁcation 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.


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 wellbeing 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 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).

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].

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 2 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)).

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: where ρ SCSC,λ is the ground reflectance for the λ-band corrected by the SCS + C method, ρ λ is the uncorrected ground reflectance for the λ band, β is the slope computed from the DSM, θ s is the solar zenith angle at the date of image acquisition, and c λ is the empirical coefficient of the λ band, calculated as the ratio of the intercept and slope of the illumination cosine regression (cos γ i ) by the ground reflectance of each band. The coefficient c λ is determined individually for each regular 5 • slope interval [50,52].

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 photointerpretation 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).

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 (DφL) is defined as the difference between the vertical and the maximum elevation angle in any of the eight directions. The angle at nadir (DψL) 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∆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: 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 • .

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: ρ N IR − ρ Red ρ N IR + ρ Red where ρ Red and ρ N IR 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

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: 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.

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: 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.

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: where v(a ik ) is the value for pixel i of the kth criterion, a, normalized to a score from 0 to 1. min i {a ik } and max i {a ik } are the minimum and maximum values of the kth 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 (1980Saaty ( , 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(A i ) = ∑ 1 k w k v(a ik ). V(A i ) denotes the overall evaluated value associated with each pixel i; w k is the weight relative to the kth attribute; and v(a ik ) 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.

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 drygrass-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. 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.

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).

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.

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.

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.

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.

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 2 for LANDSAT) reveals the largest spatial landscape scale. The Rao entropy of the 3 × 3 focal window (8100 m 2 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 2 ). 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 2 ) 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.