In the 21st century, global climate change becomes more severe which is due to greenhouse gas emissions, which are recognized as one of the key drivers of ecosystem degradation and climate disruption (Scott et al. 2020). This phenomenon has had serious consequences, including global warming, ocean acidification, accelerated glacier melt, and an increase in the frequency and intensity of extreme weather events (Calvin et al., 2023).

In the context of climate change, the uptake of carbon dioxide by forest ecosystems is precarious for regulating it (Friedlingstein et al., 2022). They play a key role since maintaining and increasing the sink capacity of forests is essential to reduce growing greenhouse gas emissions into the atmosphere (Friedlingstein et al., 2022; R. B. Myneni et al., 2001; Pan et al., 2024; Schilling et al., 2012). In December 2015, the COP 21 in Paris led to an agreement within the United Nations Convention on Climate Change with the purpose of keeping the increase in global surface temperature well under 2°C, while pursuing efforts to limit the rise to 1.5°C. In this concern, each party involved in the agreement has to establish a national goal to limit greenhouse gas emissions (Erickson & Brase, 2020; Ourbak & Magnan, 2018). Biomass carbon pools act as a sink for atmospheric CO2 and, in the Mediterranean region, carbon sequestration by forests ranges between 0.01 and 1.08 t C ha⁻¹ annually (Merlo & Croitoru, 2005). The ability to quantify forest carbon stock at the regional and local levels is expected to support compliance with the treaty and its goals.

Forests are a vital source of ecosystem services for humans and mainly act as carbon sinks (FAO, 2020). Nonetheless, forest improvement activities and changes in land and forest use emanate directly from forests and account for all emissions from agriculture and other related uses (Laaribya et al., 2024; Nourelbait et al., 2016; Rudel et al., 2005). In addition, activities linked to deforestation, reforestation, and forest conservation are important. Combined with the effects of deforestation and acceptable sustainable harvesting, forests can also act as a source of carbon long before the agreement. In this context, the reduction of greenhouse gas emissions from deforestation and increased forest degradation is part of a sustainable development approach and enhances carbon storage (Alaoui et al., 2021; Forsell et al., 2016; García et al., 2010; Laaribya et al., 2021; Sinha, 2022).

Although much research has been carried out on the Atlas cedar forest to assess its state of conservation, much remains to be discovered about its capacity to store carbon in biomass and the long-term sustainability of this emblematic Moroccan ecosystem (Boulmane et al., 2015; El Mderssa, 2022; El Mderssa et al., 2019; Laaribya, 2024; Laaribya et al., 2024; Linares et al., 2011; Terrab et al., 2006; Zaher et al., 2020a). This work has highlighted the need to improve conservation strategies to preserve this ecosystem, as its ability to act as a carbon sink is highly dependent on its sustainability and maintenance. Indeed, this remarkable ecosystem plays an essential role in carbon sequestration, helping to mitigate climate change.

The aim of this study is to use the available machine learning techniques, adapted inside the GEE environment, to assess the cover dynamic and to predict the above-ground biomass of the Azrou Cedar Atlas forest in the Middle Atlas Mountains in Morocco.

Satellite imagery, coupled with the power of Artificial Intelligence (AI) and cloud-based platforms like GEE, has revolutionized the way environmental monitoring is conducted, making it possible to analyze vast forest landscapes over extended periods efficiently (Laaribya & Alaoui, 2025; Mutanga & Kumar, 2019; Zhao et al., 2021). Indeed, given that traditional methods are difficult to meet the requirements in this field due to the long period of experimentation in the field, the availability of timber cuttings, and the high cost. Nowadays, machine learning (ML) is emerging as a new research paradigm to facilitate research in the field of machine learning for forest biomass prediction.

The Azrou forest is located on the northern edge of the Middle Atlas plateau (Morocco) and covers an area of 17,807 ha. Contrasting relief characterizes it, with altitudes ranging from 1,100 m to 2,100 m. Precipitation is relatively high and comes in the form of rain or snow. Annual rainfall varies between 563 mm and 1122 mm, while maximum temperatures range from 30.3°C to 43°C, with July and August being the hottest months (Laaribya, 2024).

The bioclimate is humid Mediterranean with a cold variant and subhumid with a temperate variant. Atlas Cedar (Cedrus atlantica Manetti) is the main species in this forest, and depending on the nature of the substrate, it forms pure stands or a mixture with holm oak (Quercus ilex) (Laaribya, 2024). The topographic characteristics of the study area are shown in Fig. 2.

Referring to the International Soil Classification System (WRB 2014), the study area offers three main soil groups (Fig. 3).

In our study area, analysis of the Gaussen Index (Bagnouls & Gaussen, 1953) (Fig. 4) reveals a dry period lasting approximately four months, from June to September (1985-2022). This prolonged aridity significantly affects vegetation cover and tree growth in the Azrou forest.

To achieve our objectives, we used a dual approach to analyze environmental changes and biomass evolution over time (Fig. 5). This study relies on various data from reliable and verified sources. All thematic maps were produced using software tools ArcGIS 10.8.

The accuracy of the forest stand types classification is assessed using a confusion matrix, which compares the stand type results to a set of reference data (ground truth or other high-quality datasets).

User’s accuracy: Measures the accuracy of classification from the user’s perspective (correctly classified instances out of all instances predicted for a given class) (Eq. 1).

(1)

Producer’s accuracy: Evaluates the accuracy of the classification from the producer’s perspective (correctly classified instances out of all instances belonging to a given class) (Eq. 2).

(2)

Kappa coefficient (K): A statistical measure that assesses the overall accuracy of the classification, accounting for random chance (Eq. 3).

(3)

Monitoring and change detection for indices used throughout the year (mean for 4 seasons) all over the 2001-2024 period.

For analyzing the vegetation conditions across time, global MODIS vegetation indices (NDVI, EVI, LAI and FPAR) were used (Table 1). The two indices provide insights into vegetation health and productivity.

Normalized Difference Vegetation Index (NDVI): Used to assess vegetation density and health, where higher values correspond to denser vegetation.

The Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) is the most commonly used vegetation index for observe greenery globally (Eq. 4):

(4)

where NIR - Near-Infrared reflectance, R - Red reflectance

Enhanced Vegetation Index (EVI): Similar to NDVI but includes corrections for atmospheric and soil variations, making it particularly useful in areas with dense vegetation.

The Enhanced Vegetation Index (Huete, 1997) is an improved version of the NDVI, designed to minimize the influence of atmospheric conditions and soil reflectance, particularly in areas with dense vegetation (Eq. 5):

(5)

where: NIR: Near-Infrared reflectance, R: Red reflectance, Blue: Blue reflectance, G: Gain factor, C1: Coefficient for the red band, C2: Coefficient for the blue band, L: Canopy background adjustment

Leaf Area Index (LAI) : LAI (Eq. 6) is broadly defined as the amount of leaf area (m²) in a canopy per unit ground area (m²) (Watson, 1958). Leaf area index (LAI) is one of the most frequently used parameters for the analysis of canopy structure and is an important structural characteristic of crop monitoring and crop productivity (Behera et al., 2010).

(6)

Variables: LA: Leaf area m²), P: Ground area (m²)

Note also that if LAI is the mean leaf area per plant, and n is the plant density, then also (Eq. 7)

(7)

Variables: LA: Leaf area of a single plant (in m² or cm²),

n: Plant density the number of plants per unit ground area (e.g., plants/m²)

Fraction of photosynthetic active radiation (FPAR): Photosynthetic active radiation used by plants in the photosynthesis process. PAR knowledge can provide key inputs for modeling biomass and forestry production (Aguiar et al., 2012; García-Rodríguez et al., 2021).

The two indices LAI and FPAR were used from the MOD15A2H V6.1 (MODIS product) combining leaf area index (LAI) and fraction of photosynthetically active radiation (FPAR) in an 8-day composite dataset at 500 m resolution (R. Myneni et al., 2021).

To detect trends and changes in vegetation conditions, the following statistical methods were applied especially to NDVI and EVI indices:

Sen’s slope estimator: A non-parametric method for estimating the slope of a trend in time series data. It is widely used for trend analysis when dealing with datasets that may contain non-linear trends or outliers (Sen, 1968).

To to apply biomass prediction over 3 years (2015, 2020 and 2024), we have used the summer period (month 5 to month 9) to calculate the biomass explanatory indices NDVI, EVI, LAI, and FPAR. The median was used to perform all those parameters. Indeed, the summer period is generally the best time to calculate the values of these indices, making it easier to identify patterns, assess vegetation health, and monitor changes.

Given the model’s robustness in prediction, we have used the Random Forest Machin Learning algorithm. The dependent variable is biomass 2010 (ee.ImageCollection(‘NASA/ORNL/biomass_carbon_density/v1’). This is the carbon stock density of the above-ground living biomass of the combined woodland and herbaceous cover in 2010. This includes carbon stored in living plant tissues above the earth’s surface (stems, bark, branches, and twigs) (Spawn et al., 2020).

The random forest is an ensemble learning method mainly used for modeling. Its principle is to build a multitude of decision trees during training and merge their results to improve overall accuracy and control overfitting (Schonlau & Zou, 2020). Random forests are a combination of tree predictors such that each tree depends on the values of a random vector sampled independently and with the same distribution for all trees in the forest (Breiman, 2001). The model parameters and their characteristics are presented in Table 1 below. Other parameters (excluding indices) are also included in the Random Forest model as independent variables.

Biomass = f (NDVI, EVI, LAI, FPAR, LULC, Elevation)

Var dataset = ee. Image. cat([NDVI, EVI, LAI, FPAR, LULC, Elevation])

Analysis of the LULC map shows that our study area is marked by a diversity of vegetation cover, mainly grassland, which accounts for more than half the surface area (57%). Forest cover appears to be open and in a state of degradation all over the study area (Fig. 6 and Table 2).

To deepen the analysis and prepare data for the prediction of forest biomass, we prepared a map of forest stand types based on data from the National Forest Inventory. An accuracy assessment was carried out to determine the validity of the classification of the results of this inventory in the field.

The composition of the forest species in our study area includes pure stands of Atlas cedar (Cedrus atlantica) (8.4%), Atlas cedar mixed mainly with holm oak (Quercus ilex) (40.3%), pure holm oak stands (24.8%) and other areas (24.7%) (Secondary species and non-wooded areas) (Fig. 7 and Table 3).

The Atlas cedar is a noble Moroccan species with a much more majestic and imposing appearance than other species.

The higher Kappa (81.55%) coefficient obtained in our analysis (Table 4) is a strong validation of the classification accuracy, allowing us to confidently focus our study on Forest stand. This robust classification framework forms the basis for assessing spatio-temporal trends in the main indices and corresponding land cover classes, in particular trees, crops, and pasture, over the selected study period (2001-2024).

The vegetation assessment parameters NDVI and EVI are widely used to analyze the condition of forest areas. According to the results obtained for the period 2001-2024, NDVI values are generally higher than EVI values over time in the study area (Figs. 9 and 10). In addition to the NDVI index, the use of the EVI index offers additional benefits by mitigating the effects of saturation and correcting for soil and atmospheric influences. The two vegetation indices complement each other and improve the detection of changes in vegetation.

Analysis of the descriptive statistics for the two series (2001-2024) confirms the results of the LULC classification, where vegetation cover is generally sparse and in a degraded state. The coefficient of variation varies by 13 and 15% for NDVI and EVI, respectively (relatively low variability), with relatively low mean values of 0.53 and 0.27 (Table 5).

The coefficient of variation varies from 32% to 17% for LAI and FPAR indices, respectively, with relatively low mean values of 0.99 and 4.13 (Table 5).These results show that the LAI index is the most variable, reflecting direct changes in leaf area over time or space. FPAR is slightly more variable than NDVI and EVI but less than LAI, representing small fluctuations in vegetation productivity (Figs. 11 and 12).

In conclusion, overall vegetation cover and greenness in the study area remain relatively low and stable in space and time in the period 2001-2024.

The differences in dynamics between the two indices (NDVI and LAI) are normal, as they are sensitive to different vegetation characteristics. NDVI reflects chlorophyll content and greenness, but it reaches saturation in dense or mature vegetation. However, LAI continues to increase with leaf growth and vegetation cover stratification, linking it more directly to leaf areas and biomass. NDVI reacts more quickly to greening at the beginning of the season, while LAI shows more gradual and sustained growth. During senescence, NDVI decreases more rapidly, while LAI continues to increase until significant leaf loss occurs.

For further statistical evaluation, we applied Sen’s slope spatio-temporal trend analysis to both the NDVI and EVI series (2001-2024). This method was chosen for its robustness in detecting monotonic trends, making it particularly suitable for analyzing vegetation dynamics over time. The results of this analysis, detailed below, offer an explanation for the spatial evolution of vegetation over the study period. A summary of results is presented in the following Table 6.

Spatio-temporal analysis carried out over the entire study area reveals both positive and negative trends in vegetation dynamics (NDVI and EVI) (Figs. 13 and 14). These trends vary and cover the entire study area. The decreasing values of Sen’s slope in the study area confirm the findings of forest degradation and the impact of climate change in the area. The two vegetation indices complement each other and improve the detection of changes in the study area.

Degradation is occurring mainly in forest ecosystems conquered by Atlas cedar (Cedrus atlantica), as well as in mixed stands of Atlas cedar and holm oak (Quercus ilex). These forest ecosystems are predominantly vulnerable due to a combination of natural and anthropogenic pressures.

Spatio-temporal analysis carried out over the entire study area reveals both positive and negative trends in vegetation dynamics (NDVI and EVI) (Figs. 13 and 14). These trends vary and cover the entire study area. The decreasing values of Sen’s slope in the study area confirm the findings of forest degradation and the impact of climate change in the area. The two vegetation indices complement each other and improve the detection of changes in the study area.

Degradation is occurring mainly in forest ecosystems conquered by Atlas cedar (Cedrus atlantica), as well as in mixed stands of Atlas cedar and holm oak (Quercus ilex). These forest ecosystems are predominantly vulnerable due to a combination of natural and anthropogenic pressures.

Biomass estimation models based on remote sensing data (NDVI, EVI, LAI, FPAR) are sensitive to changes in vegetation structure and vigor, which can decrease without any visible change in land cover type. Biomass modelling provided an assessment of the mass in the forest area studied, expressed in dry weight, of the woody parts (stem, bark, branches and twigs) of all living trees, excluding stumps and roots (Spawn et al., 2020). The Random Forest model designed for our prediction (correlation = 0, 7 with a p-value < 0, 05) has enabled us to obtain the first results by period (2010, 2015, 2020 and 2024) in the forest study area for data based on the satellite dataset (Fig. 15).

The results obtained showed a decrease in value (-8%) between 2010 and 2024, with a biomass of 501,172 tons/ha in 2010 versus 461,587 tons/ha predicted by our model for 2024.

In 2024, the biomass median values by species were 29 tons/ha for pure Atlas cedar, 24 tons/ha for pure holm oak, and 31 tons/ha for a mixture of Atlas cedar and holm oak (Table 7, Figs. 16 and 17).

Generally, between holm oak (Quercus ilex) and Atlas cedar (Cedrus atlantica), above-ground biomass potential depends on several factors such as region, ecological conditions (soil type, climate, elevation), stand density and tree age.

These results further confirm that Atlas cedar produces a higher above-ground biomass than holm oak, particularly under favorable conditions. These results provide a comprehensive approach to mapping biomass estimation in forestry and suggest guidelines for forest planning.

In addition to vegetation condition over time and space, this research work examines the assessment of forest biomass by machine learning algorithms in GEE. This innovative approach replaces the use of costly field investigations. The biomass values obtained are reference values for the main forest species in the area, namely Atlas cedar and holm oak.

A negative evolution was highlighted, in biomass values, between 2010 and 2024, materializing the negative trend in vegetation parameters studied in the area. In 2024, the predicted biomass is 461,587 tons, compared with 501,172 tons in 2010. This measurement is the carbon stock density of the above-ground living biomass of the combined woodland and herbaceous cover. The biomass median values by species were 29 tons/ha for pure Atlas cedar, 24 tons/ha for pure holm oak, and 31 tons/ha for a mixture of Atlas cedar and holm oak. According to the FAO (2006) in (Oubrahim et al., 2016), carbon stocks in forests in North Africa (the total carbon in biomass, dead wood, forest floor and the first 30 cm of the soil profile) were on average 64.9 tons/ha.

In the Middle Atlas cedar area, in four reservoirs, i.e., aboveground biomass, belowground biomass (roots), necromass (litter and deadwood) and the soil, carbon stocks were esteemed at 395.37 Mg/ha for the natural cedar Atlas and 76.05 Mg/ha for the cleared area. Analysis of the carbon stock distribution in the ecosystem discovered that soil was the largest reservoir. Indeed, the soil carbon stock varies from 46.4% to 93.5%, that of the biomass (aboveground and belowground) fluctuates between 4.3% and 52.7% and in the necromass, it is between 0.8 and 2.2% (Zaher et al., 2020b).

The highest carbon stocks are found in the most densely wooded areas (dense forests). This finding is confirmed by other studies on the subject (Le Clec’h et al., 2013; Oubrahim et al., 2016). In addition to aboveground biomass, assessing the contribution of forest soils makes it possible to estimate the total biomass level of the ecosystem. Forest soils are a significant reservoir of carbon; more than 40% of the total organic carbon in terrestrial ecosystems is stored in forest soils (Wei et al., 2014; Weston & Whittaker, 2004). In the banj oak forests (Quercus leucotrichophora) of the Central Himalaya, tree biomass declined by 62% from undisturbed to degraded forests, the carbon sequestration rate decreased by 73%, peaking in moderately disturbed-A forests, while total soil carbon fell by 79% (Pandey et al., 2020).

The decline in biomass values in our increasingly fragile ecosystem is attributed to several interdependent processes and factors that do not necessarily involve a change in LULC classification. Firstly, we can note the degradation of forest areas, such as the excessive logging of precious Atlas cedar wood and overgrazing that exceeds the carrying capacity, which can significantly reduce biomass even though the overall forest cover appears unchanged. Secondly, reduced tree density and stress can also lead to lower biomass estimates. Known climatic stress factors in recent decades (droughts and rising temperatures, etc.) have limited tree growth and health, thereby reducing biomass accumulation.

The increased stress on vegetation in the area was highlighted by analyzing spatial and temporal variations in vegetation indices (NDVI, EVI, LAI and FPAR). These indices are reliable indicators of vegetation health and are sensitive to changes in vegetation cover and structural properties (González‐Alonso et al., 2006; Shammi & Meng, 2021).

The negative trends observed for NDVI and EVI indices reflect a reduction in photosynthetic activity and vegetation density in the forest study area. Shortened vegetation affects carbon sequestration, biodiversity, and ecosystem services in the study area.

Models based on remote sensing and machine-learning techniques have made it possible to detect subtle changes in biomass, even in areas where LULC cover does not appear to have changed visibly. We can therefore conclude that these tools are powerful for monitoring and assessing the state of forest ecosystems beyond simple changes in land use.

The downward trends observed in biomass, particularly in cedar forests and mixed oak and cedar forests, reflect both local degradation processes and regional environmental pressures. In our area, carbon stocks vary considerably depending on the type of forest. The ecosystem is vulnerable to degradation, which reduces its carbon sequestration potential. Overgrazing and deforestation not only reduce above-ground biomass but also lead to soil erosion and loss of organic matter, contributing to a decrease in soil carbon stocks. In a regional context marked by human pressures and climate change (Del Río et al., 2017; Gómez et al., 2012; Vayreda et al., 2012), intensified land use, and difficulties in natural regeneration, similar trends in biomass decline and carbon loss are observed, suggesting that these trends may be regional. Globally, these findings are consistent with broader concerns about the declining carbon storage capacity of dry Mediterranean forests, pointing to the importance of sustainable management strategies.

It would be interesting to take into account local data validation (forest inventories or biomass measurements in the field) for a more accurate comparison. Assessment of biomass by species would be more interesting if we focused on results by station type, taking into account the main ecological (soil, climate, stand age, elevation) and local socio-economic factors. In fact, human and pastoral pressure on the environment would have a negative impact on the forest ecosystem in question.

Today, the adoption of innovative approaches offered by Google Earth Engine (GEE) combined with GIS and remote sensing tools is playing an increasingly central role in the analysis, monitoring, and management of forest ecosystems. Indeed, the use of these platforms in our work has provided very useful results for assessing the evolution of canopy dynamics and the prediction of aerial biomass in the study area.

This study reveals that between 2010 and 2024, biomass values in the Azrou forest studied showed a decline over time and space. This negative trend reflects a more general deterioration in vegetation vigor and health indicators. The main vegetation indices studied in the model, notably NDVI, EVI, LAI, and FPAR, followed descending trends. These trends are due to natural and human factors that have caused environmental stress. Overall, the results confirm a marked degradation of the ecosystem during the study period.

Significant spatio-temporal negative trends in vegetation indices and biomass levels underline the need for adaptive management strategies in the context of climate change. Future research should focus more on field investigations and the integration of socio-economic data to better understand the interactions of the studied forest ecosystem. Assessment of canopy and biomass dynamics would benefit from the integration of other environmental factors related to local sites (soil type and physico-chemical characteristics, stand structure, age, density).

Future research should also focus on integrating local socioeconomic data to better understand human-environment interactions and develop predictive models that promote effective mitigation and adaptation measures.