Assessment of forest dieback on the Moroccan Central Plateau using spectral vegetation indices

Cork oak forests in Morocco are rich in resources and services thanks to their great biological diversity, playing an important ecological and socioeconomic role. Considerable degradation of the forests has been accentuated in recent years by significant human pressure and effects of climate change; hence, the health of the stands needs to be monitored. In this study, the Google Engine Earth platform was leveraged to extract the normalized difference vegetation index (NDVI) and soil-adjusted vegetation index, from Landsat 8 OLI/TIRS satellite images between 2015 and 2017 to assess the health of the Sibara Forest in Morocco. Our results highlight the importance of interannual variations in NDVI in forest monitoring; the variations had a significantly high relationship (p < 0.001) with dieback severity. NDVI was positively and negatively correlated with mean annual precipitation and mean annual temperature with respective coefficients of 0.49 and − 0.67, highlighting its ability to predict phenotypic changes in forest species. Monthly interannual variation in NDVI between 2016 and 2017 seemed to confirm field observations of cork oak dieback in 2018, with the largest decreases in NDVI (up to − 38%) in December in the most-affected plots. Analysis of the influence of ecological factors on dieback highlighted the role of substrate as a driver of dieback, with the most severely affected plots characterized by granite-granodiorite substrates.


Introduction
Forests cover about a third of the Earth's land surface (Keenan et al. 2015) harbor as much as 80% of all terrestrial biodiversity, and the associated animals, plants, fungi and bacteria provide critical ecosystem services and invaluable economic benefits. However, many threats and stress factors degrade forest ecosystem and thus the overall forest health (Gauthier et al. 2015;Mate and Deshmukh 2016), which reflects the relative condition and resilience, productivity and sustainability of the forest ecosystem. Varying in magnitude, severity, origin, and frequency (Sugihara et al. 2006;Schoene et al. 2007), forest disturbances can be from factors such as forest fires, storms, pests, diseases, air pollution, and changes in climate (Larson and Franklin 2010;Zhou et al. 2012;Marini et al. 2017) and from logging and fuelwood collection, agriculture, overgrazing, and other human pressures. Inappropriate forest management strategies and socioeconomic factors have also led to problems such as forest encroachment, thus accelerating forest degradation.
Abstract Cork oak forests in Morocco are rich in resources and services thanks to their great biological diversity, playing an important ecological and socioeconomic role. Considerable degradation of the forests has been accentuated in recent years by significant human pressure and effects of climate change; hence, the health of the stands needs to be monitored. In this study, the Google Engine Earth platform was leveraged to extract the normalized difference vegetation index (NDVI) and soil-adjusted vegetation index, from Landsat 8 OLI/TIRS satellite images between 2015 and 2017 to assess the health of the Sibara Forest in Morocco. Our results highlight the importance of interannual variations in NDVI in forest monitoring; the variations had a significantly high relationship (p < 0.001) with dieback severity. NDVI was positively and negatively correlated with mean annual Moroccan forests are no exception to degradative pressures of varying intensity depending on the region. Two thirds of the rural population living in or near the forest exert immense pressure through harvesting of firewood, overgrazing, clear-cutting for agriculture (4500 ha a − 1 ) and urbanization (1000 ha a − 1 ) (Machouri 2010). Among the most affected ecosystems are the economically important cork oak (Quercus suber) forests. Its suberous bark, the cork, are in high demand as raw materials by local craftsmen and especially by the wine industry (Roula 2005); however, the cork oak forest cover has been successively degraded by human and natural factors including inappropriate management.
Forest health may be monitored by focusing on the number of forested acres, tree growth rate, plant condition and diversity, and animal species supported by the forest ecosystem. However, measuring the many components of forest ecosystems for each forest is impractical, so specific indicators of forest health are assessed. Among the most popular indicators, and fundamental for analyzing forest development, are vegetation indices, which are derived from remotely sensed data, collected by satellite sensors that measure the wavelengths of light absorbed and reflected by green plants (Hmimina et al. 2013;Riva et al. 2017;Bolton et al. 2020;Elhag et al. 2021).
Satellite remote sensing has been proven to offer a more practical and cost-effective alternative to field surveys in forest health monitoring studies, including detection of interannual variation in forest species phenology (Melaas et al. 2013;Shen et al. 2014). Stress to trees from various disturbances alters biochemical and structural properties, causing changes in leaf reflectance that can be recorded by imaging sensors to monitor forest health. With the continued development of imaging sensors, this technique has become a useful and practical tool for characterizing forest health in a variety of ways and with high accuracy. With the availability of remote sensing data from satellite sensors such as moderate resolution imaging spectroradiometer (MODIS), advanced very high-resolution radiometer (AVHRR), Landsat, Sentinel, among others, the technology is even more useful for developing vegetation indices such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), especially in broadleaf forests (Heumann et al. 2007(Heumann et al. , 2013Anees et al. 2013;Zhang et al. 2017;Barka et al. 2019;Ebinne et al. 2020).
For arid and semi-arid ecosystems, the soil adjusted vegetation index (SAVI) has been shown to provide more accurate results than the other two indices (Rodríguez-Moreno and Bullock 2014; Tariq et al. 2020). Despite the popularity of vegetation indices, their use, particularly in semi-arid and arid lands, can lead to uncertainties and potential biases in monitoring forest health due to the sparse but complex vegetation (Walker et al. 2015) and land heterogeneity (Matsushita et al. 2007).
Given the importance of Moroccan forests, monitoring their health is imperative to ensure their integrity and that of the associated ecosystems. Identifying vulnerable areas will help guide management strategies, including those for rehabilitation and full protection. To this end, we used remote sensing and meteorological data to monitor the health of the Sibara Forest in Morocco. We also sought to verify the applicability of using Landsat-derived NDVI and SAVI indices with meteorological indices to describe the dynamics that influence the health of the cork oak stands.

General description
The study area includes the cork oak stands in the Sibara Forest on the Moroccan Central Plateau (Fig. 1). It covers 8436 ha between the High Meseta of the plateau region and the Low Meseta of the Atlantic plains along an altimetric gradient of 50-600 m. The bioclimate characterizing the study area is semi-arid with temperate winters, with a mean annual precipitation of 390-600 mm. The mean monthly temperature varies from 12 to 25 °C, the maximum temperatures especially in the months of July and August range from 26 to 34.8 °C, and the minimum range from 4 to 6 °C. Geologically, the Sibara forest is in an area comprising a mountain chain composed of quartzites, shales, sandstones and dating from the Paleozoic era and a mostly granite plateau dating from the post-Stephanian era. The dominant soil types are brown soils and weakly evolved soils. The vegetation is dominated by cork oak in relatively well-preserved young to adult stands relative to cork oaks in other forests of the Central Moroccan Plateau. In addition, stands of Barbary thuja (Tetraclinis articulata) are found in the northern part of the forest, farther down on steep slopes, particularly on south and southeast aspects.

Drivers of degradation and decline
The Sibara Forest is considered one of the most spectacular cork oak forests in Morocco. However, the health of cork oak stands a worrying decline of varying severity in recent years. Much like other cork oak forests in Morocco, the decline of Sibara Forest is manifested by dieback and degradation resulting from human action, through land clearing, overgrazing, and firewood removal, combined with water stress due to recurrent droughts. In addition, attacks by pathogens such as Hypoxylon mediterraneum, which causes charcoal disease, aggravates the degradation of cork oaks, furthering their general decline.
A dieback event in December 2017 devastated considerable areas of the Sibara Forest. Table 1 presents information on trees that were killed and those with partial dieback in selected plots of the forest.
The mortality rate of cork oak trees exceeded 66% in these plots, demonstrating the severity of the dieback and the need to identify the main factors responsible for the decline of the cork oak in the Sibara Forest.

Data acquisition
Based on the locations of the observed dieback in the Sibara Forest, 83 plots of cork oak were identified for collecting relevant data for monitoring the health of the forest and elevation, slope and substrate were recorded for each plot. In addition, dendrometrics including stand age, stand density, height, and circumference at breast height were determined. Based on the observations and stand data for the plots, three classes of dieback severity were categorized: low severity (% dead trees < 10%), moderate severity (10% < % dead trees < 50%), and high severity (% dead trees > 50%).
The remote sensing data used in this study were obtained from the NASA Landsat 8 OLI/TIRS satellite, which covers the Earth every 16 days, producing 185 km × 185 km, 16-bit images with 11 spectral bands: 9 in the visible (8 multispectral with 30 m resolution; 1 panchromatic at 15 m) and 2 thermal (60 m). These data are available on the Google Earth Engine (GEE) platform (Gorelick et al. 2017). Working  with satellite imagery often involves retrieving, selecting, and downloading large-scale images, which is challenging because manipulating and analyzing these mega datasets is a time-and energy-consuming process. GEE offsets this difficulty by providing the ability to access a large satellite image data set and cloud computing, allowing users to analyze changes on the Earth's surface in real time in a continuous, dynamic, and free manner. Accordingly, image data from Landsat image collections corresponding to Landsat 8 OLI/TIRS covering all months for the years 2015-2017 were analyzed to extract vegetation indices. The respective bands used for extraction for the vegetation indices used in this study are presented in Table 2.
Monthly meteorological data from 2015 to 2017 were obtained from the Laghwalem weather station near the study area. Mean annual precipitation (P recip ), mean minimum (T min ), maximum (T max ) and monthly temperatures (T mean ), relative humidity (RH) and reference evapotranspiration (ET 0 ) were calculated and used for trend analysis in relation to the selected vegetation indices.

Data processing and analysis
Landsat images are prone to distortion due to the effects of the sensor, sun, atmosphere, and topography (Young et al. 2017). Preprocessing the image data using layer stacking and geometric, radiometric, and atmospheric corrections and other methods can minimize these effects to the extent desired for a particular application. But the process is tedious and can introduce error. The GEE platform used in this study facilitates these tasks. The first step consisted of geometric correction involving the georeferencing and orthorectification components to ensure the exact positioning of an image. Then absolute correction was applied to obtain surface reflectance, which allows for vegetation indices to be compared between images. These steps include conversion to radiance to obtain the radiance at the sensor, solar correction to give the reflectance at the top of the atmosphere, and atmospheric correction to provide the surface reflectance (Young et al. 2017). In the GEE environment, the calibrated Radiance and TOA methods were used for sensor radiance and TOA reflectance processing of the Landsat images, respectively, while the quality assessment (QA) band (pixel_qa) was leveraged for cloud masking.
In this study, the normalized difference vegetation index (NDVI: Rouse et al. 1974) and the soil adjusted vegetation index (SAVI: Huete 1988) were used to assess the health of the Sibara Forest. NDVI is widely used to analyze vegetation health and has been shown to be a relevant indicator of plant greenness (Xu et al. 2012), while EVI is designed to enhance the vegetation signal with better sensitivity, especially in high biomass areas. SAVI is often used in arid and semiarid regions to help reduce the influence of soil reflectance. NDVI and SAVI are based on the visible and near-infrared Landsat 8 spectral bands and calculated as: where, L1 is set at 0.5 and is a soil brightness correction factor.
We computed interannual variation in NDVI and SAVI at annual, seasonal, and monthly scales. To evaluate the empirical relationships between these vegetation indices and the meteorological variables, we used a multivariate linear regression approach. R version 2 (R Core Team 2021) in conjunction with GEE was used for extracting and mapping the selected vegetation indices.
To identify the relationship between the decline of cork oak and the vegetation indices, we analyzed the change in NDVI between the same months in 2015-2016 and in 2016-2017 (the year in which the decline was observed). The NDVI was chosen because it is a reliable ecological indicator for assessing forest and non-forest vegetation worldwide (Adole et al. 2016;Huang et al. 2021;Soubry et al. 2021).
For testing for a relationship between dieback and ecological, pedological, dendrometric and climatic factors that influence the emergence or progression of the dieback, maps of the different factors were superimposed on the map of dieback in the plots. Figs. S1-S5 show the maps for these factors in the Sibara Forest.

Statistical analyses
The nonparametric Kruskal-Wallis one-way ANOVA test (Kruskal and Wallis 1952) was used to test the effect of interannual (2016-2017) variation in NDVI on dieback severity in the Sibara Forest. This choice was based on the small number of plots selected for the study and because the data did not satisfy the assumptions of normality (Table S1) and homogeneity of variance (Table S2). Subsequently, Dunn's tests were used for multiple comparisons of the means for the dieback severity classes. The significance level was p < 0.05. In addition, the multiple correspondence analysis (MCA) method was used to determine any relationships between site characteristics (stand age/density, slope, elevation) and interannual variation in NDVI and cork oak dieback. The R environment was used for all analyses.

Monthly variation
The results of the mean monthly variation in the studied vegetation indices and the meteorological variables are presented in Table 3. The spatial variation in NDVI and SAVI over the study area is shown in Figs. 2 and 3, respectively. NDVI varied between 0.28 (October) and 0.53 (February) in 2015 (Fig. S6). The immediate observation was that NDVI  SAVI followed a trend similar to that of the NDVI, with 2016 values marginally lower from the values in to 2015 (Fig. S7). They increased in 2017 from the beginning of the year through October, then fell in November and December with values of 0.15 and 0.16, respectively, compared to 0.19 and 0.22 recorded in 2016. The general observation was that SAVI values were generally low throughout the study period, never exceeding 0.32.
The mean monthly variation in meteorological variables (Figs. S8-S11) showed that the mean monthly precipitation (P recip ) across the months increased between 2015 and 2016, then decreased in 2017 (Fig. S8). Generally, it was the higher in 2016 compared to the other years, particularly in the last three months, with P recip values of 49.80 mm (October), 63.33 mm (November), 89.67 mm (December). 2017 was the driest year, with the highest P recip value recorded at 73.80 mm (February). Mean monthly temperature (T mean ) showed little variation (Fig. S9). Nevertheless, both the highest (T max ) and lowest (T min ) T mean values were recorded in 2015 at 29.26 °C (July) and 7.18 °C (February), respectively.

Seasonal variation
The mean seasonal variations in both vegetation indices and meteorological variables are presented in Table 4. The seasonal spatial variations in NDVI and SAVI for the study area are shown in Figs. S12 and S13, respectively. Due to lack of meteorological data for the last 4 months (autumn) of each year, the respective data were computed only for winter, spring and summer.
The immediate observation for each year is that NDVI is highest for the winter and spring seasons and lowest during the warmer, drier summer and autumn periods. The highest value (0.46) was obtained for spring 2015, the lowest (0.28) in summer 2016 and summer and fall 2017. Across all study years, NDVI values were significantly lower in 2016 than in 2015, with marginally higher values in 2017. SAVI followed the same trend; values were higher in 2015 and 2017 than in 2016. The highest SAVI (0.29) was recorded for spring 2015, the lowest (0.17) for autumn 2017. Like NDVI, SAVI was generally higher at the beginning in 2015 compared to the last year.
Precipitation tended to increase between 2015 and 2016 and was highest in winter 2016 (60.18 mm), then decreased substantially in 2017. Changes in T mean were marginal over the years. The lowest T mean was recorded for winter 2015 (9.91 °C), and the highest was in summer 2017 (26.17 °C).

Annual variations
The mean annual variations in the vegetation indices and meteorological variables studied is presented in Table 5. The annual spatial variation of NDVI over the study area are shown in Fig. 4. NDVI is highest in 2015, with a mean value of 0.39. Conversely, it was lowest in 2016 at 0.32. A marginal increase was observed in 2017 at 0.34.
The spatial distribution of SAVI in the Sibara Forest for the study years is illustrated in Fig. 5. A similar pattern was observed for SAVI, which decreased slightly from 0.24 in 2015 to 0.20 in 2016, followed by marginal increases in 2017. As with NDVI, the highest SAVI value was observed in 2015 at 0.24, while the lowest was obtained for 2016.
Mean annual precipitation was highest in 2016 with 422.67 mm and lowest in 2017 with 234.73 mm. Its overall trend was a significant increase between 2015 and 2016, followed by a substantial decrease in 2017. On the other

Relationship between vegetation indices and meteorological variables
The results of the multiple linear regression analysis to assess the relationship between vegetation indices and meteorological variables are presented in Table 6. NDVI correlated strongly with SAVI, with a correlation coefficient of 0.89. In addition, it correlated relatively well with mean annual precipitation and relative humidity, with correlation coefficients of 0.49 and 0.53, respectively. Conversely, NDVI had the lowest negative correlation with mean monthly temperature, with a coefficient of − 0.67. Similar to its relationship with temperature, NDVI recorded a negative correlation with ET 0 , with a coefficient of − 0.46. Similar to NDVI, but to a lesser degree, SAVI was negatively correlated with mean monthly temperature and ET 0

3
with coefficients of − 0.31 and − 0.05, respectively. With respect to mean annual precipitation and relative humidity, it was positively correlated, with coefficients of 0.22 and 0.35, respectively.

NDVI change detection between 2016 and 2017
Based on the author's observations during the field survey in early 2018 of cork oak dieback in several plots in Sibara Forest, an analysis of monthly NDVI changes was performed, focusing mainly on the period between 2016 and   Table S3.

Relationship between NDVI change and dieback
To determine the ability of NDVI to predict dieback of oak trees in the forest, we used a nonparametric one-way ANOVA. The results presented in Table 7 show that the monthly interannual variation in the NDVI was very highly (p < 0.001) influenced by dieback severity.
Consequently, the Dunn's pairwise multiple comparisons test highlighted very highly significant differences (p < 0.001) in NDVI variation between the high/low and low/ moderate dieback severity plots (Table 8, Fig. 6). Consistent with the former, but to a lesser extent, NDVI variation was statistically different (p < 0.05) between plots with moderate and high dieback instances of cork oak.

Monthly interannual variations in NDVI
The percentage change in NDVI for each plot is presented in Table S4, and the spatial presentation is shown in Fig. 7. The largest positive changes in NDVI between months from 2016 to 2017 occurred in February and March. The three largest increases in February were 73.19%,76.69%,and 77.11% for plots 70,60,and 62,respectively,and in March were 49.84%,50.71%,and 58.92%,in plots 7,6,and 5,respectively. These values indicate potential recovery in vegetation in the earlier months of 2017. Conversely, the greatest decreases in NDVI occurred at the end of the year, in November and even more so in December (Fig. 8), indicating another probable case of dieback. For November, the three most affected plots were 40, 41, and 42 with NDVI decreases of − 32.86%, − 31.20%, and − 30.77%, respectively. The greatest negative monthly change in NDVI between 2016 and 2017 was in December in plots 33, 40, 41, 70, and 71 (− 38.64%, − 38.49%, 38.37%, − 38.01%, and − 37.93%, respectively), thus indicating these were the most

Ecological factors influencing the observed dieback event
The MCA to assess the influence of ecological factors on dieback showed that substrate and NDVI change were the factors most closely associated with dieback ( Fig. 9). Slope, stand age, and stand density apparently did not influence the health of cork oak. Substrate seemed to be the most influential; the most-affected plots were located on granitegranodiorite. On the other hand, the least affected plots were dominated by schist and granodiorite substrates.

Discussion
In this analysis of changes in the vegetation indices NDVI and SAVI from 2015 to 2017, extracted using Landsat 8 images and GEE, we demonstrated that these methods and indices can be used to assess the health of cork oak stands in the Sibara Forest in Morocco. For the years studied, the overall trend in the monthly NDVI and SAVI was a significant decrease in 2016 compared to 2015, which was followed by marginal increases in subsequent years, although NDVI in the final two months of 2017 were still lower than in the same months in 2016. These findings confirmed our field observations of dieback in various plots of the Sibara Forest, potentially as a result of water stress between 2015 and 2017. Dry periods lead to a decrease in canopy water content, and a sustained water deficit can reduce photosynthetic activity and lead to a loss in canopy greenness (Liu et al. 2021). This loss in greenness seems to be best indicated by greenness indices such as NDVI and SAVI, used in this study, as shown by Moreno-Fernández et al. (2021) in the study of pine stand dieback in Spain. The susceptibility assessment of dieback in the forest showed that a significant portion of the Sibara Forest is susceptible to significant dieback. Among the factors that most significantly influence the model's prediction are substrate and NDVI, thus confirming NDVI as a viable variable for predicting and thus assessing stand health. Substrate was the most significant factor in modelling, with findings on the most affected plots by dieback to be over a granite-granodiorite substrate, perhaps because the soils that develop on granite are highly prone to erosion (Deng et al. 2019;Sun et al. 2021) and the substrate has low water-holding capacity. These factors inhibit plant growth and lead to their failure in the stand, reflected by the decrease in NDVI values. Among the vegetation indices studied, NDVI had the strongest correlation with mean annual precipitation and temperature. Indeed, NDVI is effective for estimating plant biomass because it is based on the reflectance of photosynthetic tissues, which can be strongly influenced by interannual variations in these two climatic variables, in turn leading to damaging tree stress (Gitelson et al. 2002;Chen et al. 2006). The effect of the climatic variables on NDVI can also vary depending on the terrain and land cover (Thavorntam and Tantemsapya (2013). NDVI and precipitation over cropland and mixed deciduous forest were strongly positively related, but no correlation was evident for dry evergreen forest. However, increased precipitation does not always correlate with vegetation indices; Ebinne et al. (2020) reported a negative correlation between monthly precipitation and NDVI. Temperature has also been shown to have a strong negative correlation with NDVI (Wang et al. 2003;Pei et al. 2019); warmer temperatures cause more evapotranspiration, which limit plants growth. On the other hand, higher temperatures in colder regions in China, where mean monthly temperature was positively correlated with NDVI, were shown to increase plant metabolism (Duan et al. 2011).
Despite its documented benefits, remote sensing cannot completely replace traditional field studies. Passive optical remote sensing is sometimes inadequate for reliably estimating biophysical properties of the canopy and stand structure. Additionally, most studies to assess forest health are limited in their geographic scope, focus on experimental sites where the complexity of the forest environment is low, and largely use empirical approaches that require some form of site-specific adaptation, which can be tedious.

Conclusion
The investigated vegetation indices, particularly NDVI, were adequate indicators of dieback in the cork oak stand in the Sibara Forest; the monthly variation in NDVI between the years 2016 and 2017 confirmed field observations of low, moderate and high severity dieback in the forest. Although limited in scope, with respect to the area covered, our study demonstrated the potential of using vegetation indices to assess spatiotemporal changes in forest health in the region. Our findings on the decline of cork oak stands in Sibara Forest are expected to provide a basis for developing suitable regional protection and conservation practices and another tool for decision-makers to ensure optimal management of one of the most important cork oak forests in Morocco and aid future research on the drivers of forest degradation in similar ecosystems. Also needed are comparative studies that focus on dieback factors in different ecological contexts and/or explore other forest species to improve the techniques used in this study. The use of sophisticated tools such as unmanned aerial objects and the incorporation of the associated data could also prove to be useful for monitoring Moroccan forests in a changing climate. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.