Forests are ecosystems, with the concept of ecosystems being that of ecological systems. An ecosystem is defined as a system made up of biotic (animals and plants) and abiotic (bedrocks, soils, topography, microclimate) components that are in interaction with one another  . External factors significantly affect a forest ecosystem. Such factors include climate and anthropic actions, which are the core driving forces acting on Maamora forest in Morocco.
With a current total area of about 133,000 ha (including enclaves), Maamora forest used to be considered as the world’s largest one stand oak forest  . Indeed, the potential area of cork oak in the Maamora plain is considered to be around 300,000 ha. However, the actual cover is considered to be effectively around 50,000 ha, with a highly degraded, almost absent climax. In fact, cork oak currently occupies mostly the western section of the forest where the ecosystem is more resilient due to an influence from the Atlantic Ocean and the sandy nature of soil. The eastern section experiences a drier climate and is characterized by harder substrates, where cork oak has almost completely disappeared  . These limiting conditions, coupled with anthropic factors such as overgrazing and acorn harvesting  , result in a complex system that is the Maamora forest, with parameters that are increasingly difficult for managers to control.
This forest has been subject to several management plans since 1951, and the area occupied by cork oak has been considerably reduced in favor of pine, acacias and eucalyptus species that progressively occupied an important area of the forest   . Fortunately, some regeneration from seed started to be successful in recent times, giving new hope to the recovery of cork oak in Maamora forest  .
Some recent studies have addressed the issue of site suitability for cork oak regeneration in Maamora forest   . However, no work has been conducted to assess the vulnerability of forest ecosystems to global changes. Taking into account the changes in climate in the last decades, including the global warming of the climatic system  , the evolution of population, and thus the joint consequences of both climate and human activity on Maamora forest, one question ought to be posed: Can the consideration of vulnerability of Maamora forest ecosystems to global changes help in recovering its oak forest to its initial state? It might be possible that some plots of the forest consisting of other species are more favorable to the development of cork oak and that oak regeneration failure on other plots can be avoided simply by accounting for the vulnerability to global changes (climate, site, human activity).
Therefore in this work, we first aim to better understand forest dynamics before assessing the vulnerability of forest ecosystems to global changes in the second part, taking into account the most influencing factors in the Maamora forest context.
2. Materials and Methods
2.1. Study Area
Maamora forest is located in North-eastern part of Morocco (Figure 1). Its 5 blocs subdivided in 33 groups and 460 parcels range over 68 km (West-East) and 38 km (North-South)  .
The climate in Maamora forest is Mediterranean, with an influence from the Atlantic Ocean. The bio-climate ranges from semi-arid with temperate winters in the eastern part of the forest to sub-humid with warm winters in its western part  . The average seasonal rainfall regime is W.A.S.S (Winter-Autumn- Spring-Summer). The topography of Maamora forest is characterized by a not so rugged relief, with slopes ranging between 0.6% and 0.8% in average, except the eastern part where slope is steep enough to provoke significant erosion  . The soils characterizing the forest are generally sand lying on clayey layer. The type of soil is important in water retention. In fact, with thick sand layers, water can become a limiting factor due to its unavailability for use by plants  .
The vegetation of Maamora is mainly comprised of cork oak, even though it has decreased in area lately, giving space for eucalyptus, pine and acacia plantations. Throughout its history, and with the products and services that it provides, Maamora forest is one of the most important forests in Morocco. It sig-
Figure 1. Geographic location of Maamora forest.
nificantly accounts for the needs in fuelwood, grazing for livestock, recreation, etc. for more than 10 communes  . Human activity is quite notable and mainly includes acorn collection, grazing and woodcutting for energy  .
2.2. Forest Dynamics Study
In this section, data in the form of Landsat images was obtained at different dates, from the http://glovis.usgs.gov/ website. Available forest maps were also obtained from forest administration. Satellite images for the period between June and September, downloaded as individual layers, were stacked and sharpened to produce multispectral images. That period of the year (summer) was chosen in order to prevent the confusion that could be caused by shrubs during classification. All the available images for the summer period were stacked together to produce a composite image for each year. This was essential in order to reduce radiometric noise that would be notable if a single image were used per period.
The Figure 2, below, summarizes the methodology used for studying forest dynamics.
The classification of images for the three years considered (1987, 2000 and 2014) for studying the dynamics using the SVM (Support Vector Machine) algorithm made it possible to create a map showing the distribution of the main tree species of the forest. After classification and before obtaining the final map, sieving was done, considering a minimum mappable area   of 2 ha which represented 20 pixels. A control for accuracy was performed using the Kappa index to validate the classification  . Ground truth was collected for the year 2014. However, for the years 1987 and 2000, we performed a visual interpreta-
Figure 2. Methodology for analyzing forest dynamics.
tion, assisted by available maps, of at least 60 points per density-species class which we compared to the result of the classification in order to obtain the Kappa index. Afterwards, change detection was done by comparing the land cover maps. Finally, some drivers of change were identified.
2.3. Assessment of Vulnerability
The assessment of vulnerability was done using spatial multi-criteria analysis  . Eleven individual factors grouped in four groups of factors (biophysical, climatic, anthropic and silvicultural) were considered as result of forest dynamics study:
・ Biophysical factors: Natural terrain slope, slope of the clay layer, type of soil and depth of sand;
・ Climatic factors: Water stress and isohyets (proxy for continentality);
・ Anthropic factors: need in fuelwood, grazing pressure and population density;
・ Silvicultural factors: Age of stand and health state.
Data was collected from different sources (previous studies, Direction of national meteorology and forest administration). Water stress (WS) was determined from raw data and computed using the following formula:
P is Precipitation (in mm),
AWC is the Available water capacity  ,
ETP is the Evapotranspiration  .
WS was computed for each month and then averaged for months with a water stress, i.e. when WS < 0. Each individual factor was spatialized, categorized (thresholds), averaged by group (unit of the forest) and combined at 2 levels (factor and final vulnerabilities) using the formula of Equation (2) in order to obtain the final vulnerability at the reference year, 2010 (year chosen for availability of all data needed).
V is vulnerability,
Wi are the weights,
Fi represents each individual factor or group of factors.
The method is summarized in Figure 3. The weights used to linearly combine the different factors were determined by Analytic Hierarchy Process (AHP) developed by Saaty  . The intensity of importance required for computing the weights are given by experts.
In order to forecast vulnerability at years 2045 and 2070, we projected factors considered changing over the period of study. Biophysical factors for instance were not projected given the fact that they wouldn’t significantly change over the
Figure 3. Methodology for assessing vulnerability (adapted from  ).
considered period. The climatic factor was derived by projecting temperatures and precipitations to the 2 different years considering 2 different scenarios, one we called optimistic (scenario RCP 4.5, RCP being Representative Concentration Pathway) and the other pessimistic (scenario RCP 8.5). Those scenarios were conceived based on possible trajectories of greenhouse gases  and an interface available on internet ( http://climexp.knmi.nl/plot_atlas_form.py ) was used to project climatic data. Anthropic factors were projected taking into account average growth of rural population in Morocco  . Stand age was recalculated supposing no regeneration would occur during that time.
3.1. Results of Dynamics Study
The classification of the 3 images was followed by a control. Each time, the Kappa index (written on each map) was found to be higher than 80%, thus allowing us to validate the classification. The analysis of the obtained maps (Figure 4) showed that cork oak was currently mainly present in Blocs A, B and C, with high density only in bloc A. In the other blocs (D and E), it had almost disappeared, giving space to mainly eucalyptus and pine. Of all the forest species, acacia was the least represented.
To better interpret the maps shown above, a synthesis of area was necessary in order to obtain a cover percentage per species, with all density classes put together. Figure 5 presents the trends for each species. Percentages are relative to the total area of the forest calculated excluding all enclaves, i.e. 131,020 ha.
The graph shows that:
・ The area covered by cork oak (49.2% or 64,461.38 ha in 1987) declined by
Figure 4. Species-Density maps for 1987, 2000 and 2014, with Qs = Quercus suber (cork oak), Eu = Eucalyptus, density codes 1, 2, 3, 4 = dense, less dense, clear, sparse.
Figure 5. Percentage cover of each species.
11.7% during the first period, between 1987 and 2000 (43.46% of total area or 56,937.4 ha in 2000), before increasing by 7.96% during the second period, between 2000 and 2014 (46.92% of total area or 61,471.65 ha in 2014);
・ The area covered by eucalyptus (34.13% or 44,719.49 ha in 1987) increased by 7.5% during the first period and then decreased by 11.31% during the second period (32.54% or 42,635.13 ha in 2014);
・ The area under pine plantations (4.4% or 5769.73 ha in 1987) increased significantly by 73.54% over all the studied period (7.64% or 10,012.84 ha in 2014);
・ The area under acacia (2.94% or 3849.94 ha in 1987) decreased by 33.67% during the first period before increasing by 41% during the second period (2.75% or 3600.9 ha in 2014).
That said, it is important to point out that the decrease in area of cork oak resulted in available space for the other species, including eucalyptus and pine. The increase in cover between 2000 and 2014 was most likely a result of successful cork oak regeneration by seed  . This was obtained after the decision of forest administration, in 2003, to restore the Maamora forest with its natural species which is cork oak.
These results are partly in tandem with those obtained in a study, done few years ago by the forest administration, as presented in Figure 6. According to that study, cork oak area decreased from 100,000 ha to 60,000 ha between 1951 and 1992 followed by an increase of up to 65,000 ha in 2006. In our study, we obtained an area of almost 57,000 ha in 2000. Hence, it is possible (if we place this value in the trend shown in Figure 6) that cork oak cover continued to decline until 2000, before increasing again.
However, the trend in cork oak area between 2000 and 2014, in our study does not match well with the study above mentioned. In fact, the evaluated cork oak area was 70,400 ha in 2011, whereas the classification in this work resulted in an area of about 61,500 ha in 2014. This difference could be partly explained by the omission that might be a result of very young plantations of cork oak not having sufficient crown growth. This is essential for reflecting light as trees, and helps in the identification of cork oak. Consequently, plantations for the last years might not be well accounted during the classification of the 2014 image.
The changes that occurred during the studied period are presented in the graph of Figure 7 below:
This chart illustrates the importance of the dynamics of forest ecosystems in Maamora forest. In fact, for both periods, only about 40% of the total area (or about 60% if we do not consider shifts to other density classes as changes) remained unchanged. Despite all the efforts for restoring cork oak ecosystems, a slight increase in the area shifting to lower density classes was still noticeable.
Figure 6. Area of cork oak in Maamora from 1951 to 2011  .
Figure 7. Changes for 1987-2000 and 2000-2014. Codes Qs: cork oak, Eu: Eucalyptus, symbol: =, −, +: unchanged density, shift to lower density and shift to higher density classes respectively. Unch.: Other unchanged class, Conv.: conversion to other species (plantation with change of species), defor.: Deforestation, Affor.: Afforestation.
Figure 8. Drivers of change in Maamora forest.
However, increased conversion and minor afforestation activities could, for the most part, be accounted for cork oak regeneration.
This dynamics was mainly a result of management actions, although other anthropic activities such as eucalyptus woodcutting for making charcoal, fuelwood, grazing, etc. could have contributed.
Below is a schema summarizing the most important drivers of change in the Maamora forest.
From the analysis that led to the elaboration of schema in Figure 8, 4 groups of important drivers of change in Maamora forest were identified, as follows:
・ Climatic conditions: aridity was more pronounced as we moved from West to East, with the influence of Atlantic Ocean being reduced. This factor would be more important with current trends in climate  .
・ Biophysical conditions: in Maamora, the most important biophysical factor could be considered to be the type of soil (sand lying on clay layer) as the slope of the clayey layer and the thickness of the sand layer could considerably influence water retention in the soil.
・ Health state: parasite attacks combined with aging of the forest weaken the vegetation and lead to its decline.
・ Anthropic factors: population increase coupled with increased livestock presence resulted in increased acorn picking, illegal logging and overgrazing.
These 4 groups of factors were used to assess the vulnerability of forest ecosystems to global changes.
3.2. Results of Vulnerability Assessment
The weights used to combine individual factors to obtain factor vulnerabilities are given in the table below:
The linear weighted combination of individual factors maps with the weights of the Table 1 produced the maps in Figure 9. The biophysical factor vulnerability map shows two groups (BG2 and BG7) as non-vulnerable, with no group of the forest being of high vulnerability considering only this factor. This implies that biophysical factors, even though not very favorable, are not the most limiting in Maamora. The climatic factor vulnerability map underlines the effect of continentality on vulnerability, thus highlighting that most continental ecosystems are the most vulnerable. As for the anthropic factor, unlike the biophysical factor, we do not have any non-vulnerable group, thus confirming the importance of human activity in the forest. The silvicultural factor vulnerability, on the other hand, is on average fairly low. This implies that even though Maamora forest might be old, especially in reference to cork oak, its age might not be the most important factor with regards to vulnerability of forest ecosystems to global changes.
Table 1. Weights of each individual factor, obtained by AHP.
Figure 9. Factor vulnerabilities for reference year (2010).
Table 2. Aggregation weights of factor vulnerabilities, obtained by AHP.
The final vulnerability map did not highlight any non-vulnerable group, confirming the fact that all the ecosystems were, to a certain extent, vulnerable to global changes. It is important to note that only the DGI group was highly vulnerable in 2010, whereas all the other groups were lowly to moderately vulnerable.
After projection of changing factors (climatic, anthropic and silvicultural) final vulnerabilities obtained for each scenario and depending on the considered years (2045 and 2070) are presented in Figure 11.
Final vulnerability predicted for 2045 using the optimistic scenario (for pro-
Figure 10. Final vulnerability at reference year.
Figure 11. Final vulnerability maps for the 2 scenarios and 2 considered years (2045 and 2070).
jection of climatic factors) showed that 3 groups of bloc A and 1 group of bloc B would be of low vulnerability whereas the remaining groups of blocs A and B, as well as the groups of bloc C would be moderately vulnerable. The groups of blocs D and E would be of high vulnerability. When considering the pessimistic scenario, all the groups of blocs C, D and E would be highly vulnerable whereas those of blocs A and B moderately vulnerable.
In 2070, with the optimistic scenario, the obtained vulnerability map of forest ecosystems was similar to that of 2045 with the pessimistic scenario. The groups of blocs A and B would be moderately vulnerable, whereas those of blocs C, D and E highly vulnerable. Additionally, when considering the pessimistic scenario, some groups of bloc B would shift to high vulnerability.
・ Statistics of vulnerability per category at reference year 2010 and for the horizons 2045 and 2070 considering the 2 climate scenarios (RCP 4.5 and RCP 8.5).
The graph shown in Figure 12 presents a synthesis of percentage of total forest area per vulnerability class for the reference year (2010). It also illustrates the scenarios for horizons 2045 and 2070, with reference to the 2 climatic scenarios.
The graph shows that the vulnerability of forest ecosystems would increase with time as we go from 2010 to the horizons 2045 and 2070, with this being more or less significant depending on the considered scenario. In fact, forest ecosystems seem more vulnerable to global changes with the pessimistic scenario (RCP 8.5).
One important thing to note is the fact that no group presented an insignificant vulnerability, be it at reference year or at horizons 2045 and 2070. In 2010, 53% of the entire forest area had a low vulnerability whereas only 11% in 2045 would remain with low vulnerability with scenario RCP 4.5. The rest of the forest would shift to moderate and high vulnerabilities. At horizon 2045, with scenario RCP 8.5, and at horizon 2070, regardless of the considered scenario, all the groups of the forest would be moderately to highly vulnerable.
It should be noted that a certain level of subjectivity in the thresholding (low, moderate, high vulnerability) of the individual factors remains. This could be a source of significant variability in the resulting final vulnerability, from author to author. However, with the thresholds being defined rigorously based on our best knowledge of Maamora forest, the obtained results could be considered a reflection of the vulnerability of Maamora forest ecosystems to global changes.
Figure 12. Percentage per category of vulnerability and scenario.
The study of Maamora forest dynamics pointed to the fact that the forest was experiencing significant changes. This dynamics which is mainly a result of the replacement of cork oak by other species such as pines, acacias and eucalyptus is, moreover, the result of various management actions, most of which were driven by previous regeneration failures.
The assessment of vulnerability helped highlight the forest’s continental blocs (blocs D and E) as being more vulnerable to global changes. It is, therefore, not surprising that oak almost disappeared in these forest blocs. That said, the methodology applied appears to be reliable given that it helped identify sections of the forest which were characterized by vulnerability. As such, it could be recommended, but with adaptations, in other studies of vulnerability assessment under different ecological contexts. Also, it clearly appears that ecosystems will become more vulnerable, going by projections at horizons 2045 and 2070, affecting even the currently more resilient blocs (A, B and C).
Nevertheless, this doesn’t mean that Maamora forest ecosystems will disappear by 2070. Instead, the findings of this study aim to draw the attention of forest managers on the future of Maamora forest if no actions were taken. Possible actions include the creation of income generating activities and the development of alternatives to the use of fuelwood, in order to reduce anthropic pressure on the forest. Moreover, reforestation with species other than cork oak should be limited, and artificial regeneration by seed encouraged in the western blocs (A, B and C) where cork oak ecosystems are still more resilient. Reforestation of cork oak (if keeping the other species is not to be considered as an option) ought to be more intensive, with more tending activities in the eastern blocs (D and E) where not only biophysical but also climatic conditions are harsher.
This research was financially supported by the FFEM (Fond Français pour l’Environnement Mondial) in collaboration with FAO (Food and Agricultural Organization) and HCEFLCD (Haut Commissariat aux Eaux et Forêts et à la Lutte Contre la Désertification).
 El Boukhari, E.M., Brhadda, N. and Gmira, N. (2015) Un regard rétrospectif sur la gestion antérieure et sur les principaux résultats acquis de la régénération du chêne liège (Quercus suber L.) dans la Maamora (Maroc occidental). [A Retrospective Insight on the Management History and the Main Results of Cork Oak (Quercus suber L.) Regeneration in Maamora (Western Morocco).] International Research Journal of Multidisciplinary Science, 1, 30-37.
 Laaribya, S., Alaoui, A., Gmira, N. and Nassim, G. (2014) Contribution à l’évaluation de la pression pastorale dans la forêt de la Maamora. Parcours forestiers et surpaturage. [Contribution to Grazing Pressure Assessment in the Maamora Forest. Forest Pastures and Overgrazing.] Nat. Technol., 10, 39-50.
 Aafi, A. (2007) Etude de la diversité floristique de l’écosystème de chêne-liège de la forêt de la Mamora. [Study of Plant Diversity of Cork Oak Ecosystem in the Maamora Forest.] Thèse de Doctorat Institut Agron. et Vétér. Hassan II, Rabat, 190 p.
 Bagaram, B.M., Mounir, F., Lahssini, S. and Ponette, Q. (2016) Site Suitability Analysis for Cork Oak Regeneration Using GIS Based Multicriteria Evaluation Techniques in Maamora Forest-Morocco. OALib, 3, 1-9.
 Lahssini, S., Lahlaoi, H., Alaoui, H.M., Hlal, E.A., Bagaram, M. and Ponette, Q. (2015) Predicting Cork Oak Suitability in Maamora Forest Using Random Forest Algorithm. Journal of Geographic Information System, 7, 202-210.
 IPCC (2013) Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Summary for Policy Makers. Cambridge University Press, Cambridge, New York, 33.
 Mounir, F. (2002) Conception d’un système d’information géographique pour l’aménagement et la gestion forestière au Maroc; Intégration des critères et indicateurs du dévelopement durable; Application à la forêt de la Maamora. [Designing a Geographic Information System for Forest Planning and Management in Morocco; Integration of Criteria and Indicators of Sustainable Development; Application to Maamora Forest.] Thèse de Doctorat, UCL, Belgique.
 Abourouh, M., Taleb, M., Makhloufi, M., Boulmane, M. and Aronson, J. (2005) Biodiversité et dynamique de la végétation dans la subéraie de la Maamora (Maroc) Effet de la durée de cloture. [Biodiversity and Dynamics in the Maamora (Morocco) Cork Oak Forest; Effect of Fence Duration.] Forêt Méditerranéenne, XXVI, 275-286.
 Lepoutre, B. (1967) Régénération artificielle du chêne-liege et équilibre climacique de la subéraie en forêt de la Mamora. [Artificial Regeneration of Cork Oak and Climax of Maamora Cork Oak Forest.]
 Saura, S. (2002) Effects of Minimum Mapping Unit on Land Cover Data Spatial Configuration and Composition. International Journal of Remote Sensing, 23, 4853-4880.
 Chen, J., et al. (2014) Global Land Cover Mapping at 30 m Resolution: A POK-Based Operational Approach. ISPRS Journal of Photogrammetry and Remote Sensing, 103, 7-27.
 GIZ (2013) Guide méthodologique Approche spatiale multifactorielle d’analyse de la vulnérabilité des écosystèmes face au changement climatique Cas de la subéraie en Tunisie. [Methodological Guide on Spatial Mutli-Criteria Analysis Approach to Assess the Vulnerability of Ecosystems to Climate Change: Case Study of Cork Oak in Tunisia.]
 Kili, M., El Mansouri, B. and Chao, J. (2008) Bilan hydrique des sols et recharge de la nappe profonde de la plaine du Gharb (Maroc). [Water Balance of Topsoils and Recharge of the Deep Groundwater on the Gharb Plain (Morocco).] Sécheresse, 19, 145-151.
 Bouteldjaoui, F., Bessenasse, M. and Guendouz, A. (2011) Etude comparative des différentes méthodes d’estimation de l’évapotranspiration en zone semi-aride (cas de la région de Djelfa). [Comparative Study of the Different Methods for Estimating Evapotranspiration in Semi-Arid Areas (Case Study of the Djelfa Region).] Nature Biotechnology, 109-116.
 RGPH (2004) Caractéristiques démographiques et socio-économiques de la population. Recensement Général de la Population et de l’Habitat, 2004. [Demographic and Socio-Economic Characteristics of the Population. General Census of Population and Habitat, 2004.] Haut Commissariat au Plan, Rabat-Maroc.
 Belghazi, B., Badouzi, M., Belghazi, T. and Moujjani, S. (2011) Semis et plantations dans la forêt de chêne-liège de la Maamora (Maroc). [Sowing and Planting in the Maamora Cork Oak Forest (Morocco).] Flora Mediterr., XXXII, 301-314.
 DREFLCD (2012) Evolution de la superficie du chêne liège dans la Maamora. Unpublished. [Area Change of the Maamora Cork Oak Forest.] Direction Régionale des Eaux et Forêts et à la Lutte Contre la Désertification, Rabat-Maroc.