Anthropogenic Activities and Climate Variability Impacts on Juniperus thurifera Forest of El Houanet, Morocco

The current case analysis aims to explore the spatio-temporal evolution of thuriferous juniper ( Juniperus thurif-era L.) forest formations on the El Houanet plateau within the Biological and Ecological Interest Site (BEIS) of Boutferda in Morocco, focusing on their vulnerability to extinction due to various complex reasons. Our research focused on the period from 1986 to 2022, conducting a meticulous assessment of the spatial changes observed over this chronology. To achieve these objectives, we adopted a diachronic approach, combined with a methodology supported by field examinations. We derived data, including climatic analysis and vegetation indices (NDVI and VCI), from LANDSAT 5 and 8 satellite imagery, enabling a detailed synthesis of the forest’s spatiotemporal changes. This analysis was further deepened by considering a range of parameters, including components of the physical environment, climatic variables, as well as anthropogenic influences, with the aim of refining our understanding of the dynamics governing these populations. Findings reveal significant impacts on various strata within this forest, reflecting both regressive and progressive changes that reshaped the landscape’s composition and structure.


INTRODUCTION
The forest formations within the Mediterranean region serve as vital repositories of biodiversity essential for sustaining ecological equilibrium (Quézel & Médail, 2003).Despite their significance, these ecosystems confront imminent threats from both natural and human-induced degradation, endangering emblematic species such as the Thuriferous juniper (IUCN, 2008).In Morocco, the Juniperus thurifera populations face challenges compounded by local communities' heavy reliance on these forests for sustenance and livelihoods, leading to mounting pressures (Gauquelin et al., 1999;Auclair, 1993).This study presents a dedicated examination of the El Houanet juniper forest within the Boutferda Biological and Ecological Interest Site in Morocco's Central High Atlas from 1986 to 2022.Employing a diachronic approach, our research meticulously analyzes Landsat satellite imagery to extract and evaluate NDVI data for key years (1986,1995,2004,2013,2022).
Through adept Change Detection techniques, our focus is to discern and delineate areas exhibiting significant alterations -be it degradation or regenerative growth -shedding light on the complex spatiotemporal dynamics.Embedded within our methodological framework is a deep exploration aiming to uncover the fundamental forces steering the vegetation dynamics within the El Houanet juniper forest.This endeavor not only delves into the spatiotemporal variations of precipitation and humidity but also scrutinizes the evolution of vegetation hydric conditions (1986-2022) amidst factors such as grazing impacts and excessive firewood harvesting.By addressing these multifaceted aspects, our study seeks to bridge existing gaps in understanding while offering insights crucial for devising sustainable management strategies essential for preserving this ecologically sensitive area in the Central High Atlas.

Study area presentation: the thuriferous woodland of the Boutferda BEIS
The Boutferda BEIS is located in the central High Atlas of the Azilal province, which falls administratively under the Beni Mellal Khenifra region, as illustrated on the map below (Fig. 1).It covers approximately 27.419 hectares.

Eco-climatic context of the study area: precipitation
The evolution of annual precipitation recorded at representative stations, notably Tizi N-Isly and Ait Ouchène, reveals that in 1995/96, considerable amounts were recorded, exceeding 1000 mm at the Tizi N-Isly station and 700 mm at the Ait Ouchène station (Fig. 2).In contrast, for these two stations, the years 1983/84 and 1993/94 experienced the lowest levels of precipitation, not exceeding 220 mm in either station (Fig. 2).The precipitation regime follows an WSPS type, characterized by a high concentration of rainfall in winter.Indeed, the first peak of precipitation occurs during winter for the entire study area.To perform a comprehensive spatial analysis of precipitation, we utilized the nearest climatic stations surrounding the study area.Incorporating data from these stations, three additional sources were integrated for spatial interpolation.The comparison between the deterministic method (IDW) and the geostatistical method (kriging) distinctly favors the latter, showcasing RMS values of 56.12 and 24.52, respectively.Figure 3 provides a visual representation of the precipitation's spatial distribution throughout the study area.Derived from detailed observations from five distinct climatic stations rather than two used to contextualize the variability in annual precipitation.This depiction illustrates a noticeable gradient stretching from north to south.

Temperatures
Due to a lack of station data, we utilized thermometric open source data from Climate Engine, covering the period from 1979 to 2021.The observation of the temperature trend line (Fig. 4 and Fig. 5, Table 1 and Table 2) indicates an upward trend, which could be an indicator of ongoing climate changes.However, it is important to note that this trend should be interpreted with caution, as the study period is relatively short compared to the scale of climatic timescales.Figure 5 reveals that the average maximum temperature reaches its peak in July, with a value of 36.95°C.On the other hand, January records the lowest average minimum temperature, at -10.6°C.The average temperatures are 1.68°C in January and 23.03°C in July.Overall, the annual average of maximum temperatures is increasing, indicating that summers in the study area are hot.Additionally, the climate in the area is cool and cold at the peaks and high altitudes.

Bioclimatic summary -rainfall quotient and emberger climagram
The rainfall quotient categorizes the study area into a bioclimatic range from semi-arid with a very cold variant to subhumid with a very cold variant (Table 1).

Geology
The dominant geological formation is the marine Dogger, covering an extensive area of 260.8 km², which represents a significant portion of 93.8% of the total (Fig. 6).Middle and recent Quaternary formations, primarily composed of alluvium and scree, occupy an area of 6.3 km², accounting for 2.3% of the total.On the other hand, ancient Palaeo-Quaternary, Miocene-Pliocene Tertiary, Upper Triassic, and Jurassic formations are less widespread, each occupying areas ranging from 0.5 to 4.2 km², representing respectively from 0.2% to 1.5% of the total area (IBOUH, 2004).

Altitude, slopes, and exposures
The dominant altitudes range from 2010-2250 m, covering an area of 116.9 km², or 42.4% of the total area.These characterize the central and southern parts of the study area.As for the lowest altitudes, they are below 1500 m, accounting for 7.2%, and are located in the northern part.Altitudes higher than 2250 m cover an area of 26.6 km², or 9.5% of the total area.Furthermore, the study area is characterized by rugged terrain with varied topography.Indeed, the slope of the slopes varies from 5% to over 36% in different directions, giving rise to varying exposures of the slopes (Flat, North, Northeast, East, Southeast, Southwest, ...).

Vegetation
The BEIS of Boutferda exemplifies the Oromediterranean stage, particularly in its lower series characterized by the presence of the Juniperus thurifera.This stage corresponds to the climax vegetation, similar to that of the highest altitudes of the Middle, High, and Anti-Atlas, as well as partially in the Rif (Benabid, 1982).In this  semi-arid and very cold context, species sensitive to low temperatures, such as the Evergreen Oak, are absent (Donadieu et al., 1976).The floristic composition is distinguished by a gradual replacement depending on altitude, transitioning from species like Lavandula depressus, Nepeta atlantica, Genista florida var.maroccana, and Thymus pallidus, to xerophytes such as Cytisus balansae, Alyssum spinosum, Ormenis scariosa, and Bupleurum spinosum.It constitutes a pre-steppe formation (Barbero et al., 1981).The stands formed by this Cupressaceae tree are generally considered pre-steppic (Benabid, 1985), described as "sparse tree formations" where the substrate is practically devoid of exclusively forest elements, but rather invaded by perennial species with steppe affinities, on poorly evolved soils, often superficially truncated (Barbero et al., 1990).The climatic analysis of the BEIS of Boutferda reveals a regime characterized by a precipitation variation between 383 and 436 mm/year, as well as a considerable thermal amplitude, with an extreme average value of 36.95°C.The wettest season occurs in winter, representing 40% of the precipitation, while summer is hot and dry.
In contrast to the recent analysis based on current data, which led to characterizing the pluviometric regime as WSAS type, the literature, notably Achahal et al. (1979), describes this regime as WASS type, with a strong concentration of rains in the winter period.Indeed, it is a climate with a balanced regime, where the dry period extends over 4 to 5 months, from May to September.As for the average minimum temperatures, these are very low for all stations, while the maximum temperatures remain high.This analysis leads us to conclude that the BEIS of Boutferda is situated in a bioclimatic environment ranging from semiarid to humid, with a cool to cold variant (Achhal et al., 1979), which differs significantly from our analysis, emphasizing a bioclimatic environment ranging from semi-arid to very cold variant, up to subhumid with a very cold variant.

Anthropogenic characteristics
From an anthropogenic perspective, this mountainous massif is characterized by a multitude of administrative divisions, totaling four territorial communes distributed throughout the massif.All of these territorial communes have experienced evolutionary trends during both periods (Table 2).The anthropogenic aspect holds great importance in this type of study because, in general, the impact of natural resource extractions often takes on an abusive character.

Methodology
In the context of this study, we undertook a meticulous analysis of the dynamics of the thuriferous juniper forest formation in El Houanet over a period of thirty-six years.To achieve this, we employed a rigorous methodology centered around two main components: land cover mapping and the assessment of changes occurring over time in the study area.
To do so, we opted for the use of openly accessible data.We acquired a set of five acquisitions, one acquisition per year.This included three Landsat 5 acquisitions with the TM sensor and two Landsat 8 acquisitions with the OLI sensor (see Annexes 1, 2, and 3).The images were acquired in May or April for two major reasons: (i) to ensure cloud-free images to facilitate their preprocessing, (ii) to capture the optimal vegetation period of the trees, which allows for a clearer and more precise spectral signature.The obtained data were processed using open-source software, namely: QGIS, SCP, and Google Earth.To meet the objective of our study, it was necessary to establish a land cover map following an approach aiming for maximum precision.
The aforementioned methodology consists of three main steps: preprocessing, classification, and validation.Preprocessing aims to ensure data reliability through radiometric corrections (converting pixel digital values to radiance) and atmospheric corrections (cloud, shadow, snow, and water masks using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) chain).This preliminary step is essential to guarantee the quality of processed data and the accuracy of subsequent analyses.
Furthermore, the classification aims to determine land cover classes using the Random Forest (RF) algorithm.Three land cover classes were identified, and 40 Regions of Interest (ROIs) were created for each class.Subsequently, we proceeded with an evaluation and validation phase to assess the classification accuracy.This evaluation relied on validation data (20% of the training sample), NDVI classes with desktop missions (high-resolution Google Earth), and field missions to validate spectral signatures.The ultimate goal of these steps is to produce an accurate map of variations in land cover class distribution over time.Given the relatively limited temporal dimension compared to the ecosystem scale, the introduction of the Vegetation Condition Index (VCI) is crucial to complement the analysis of vegetation dynamics, in conjunction with the NDVI.Developed by Kogan (2003), the VCI improves the meteorological component of our study.It is independent of underlying ecological conditions and relies solely on meteorological parameters (Ramesh et al., 2003).
The statistical treatment of temporal variability in NDVI data was structured around two aspects.The first aimed to determine the existence of a trend in each chronological series formed by pixels of an image.The second aspect sought to identify the presence of breakpoints in these time series.For these analyses, we employed the Mann-Kendall test, which takes into account the seasonality of a time series.Finally, it will be essential to highlight the anthropogenic pressures exerted on vegetation and soils.To achieve these objectives, we employed geostatistical methods, remote sensing, and GIS, thus enabling the analysis of climate and land cover evolution, as well as detecting correlations between climate, anthropogenic activities, and observed changes in the environment.

Land cover classification
After completing the classification, a validation phase was undertaken to assess the accuracy of the process.This was done based on validation data (20% of the training sample), NDVI classes, and spectral signatures.

Confusion matrix (model evaluation)
The confusion matrices provide a quantitative evaluation of the results of satellite image classification for different years and sensors (Landsat 5 TM and Landsat 8 OLI) (Tables 3, 4, 5, 6, and 7).

NDVI and ground truth
The NDVI vegetation index extracted and analyzed from satellite images enabled the delineation of distinct classes within this index.To fortify this initial process, we opted for the classification of three specific ranges (below 0.15 / 0.1-0.9/ above 0.9).Subsequently, field verification was carried out to validate these classes based on the actual forest conditions.This validation involved systematic sampling surveys, employing a grid measuring 2.5×2.5 kilometers and establishing 5 ares plots (Fig. 7).These surveys were focused on evaluating the density and coverage of the forest population by projecting the canopy onto a 100-meter by 20-meter transect.These field procedures have substantiated and fortified our initial image processing endeavors.The resultant findings, rooted in the NDVI index, offer empirical affirmation for our land cover classification, thereby augmenting the robustness of our methodology.Moreover, the nexus between coverage and density, coupled with the minimal standard deviations evident in the NDVI across each class, further buttresses the statistical underpinning of our approach.Note: The classification rates are generally high, with user's accuracy ranging from 95.21% to 99.26%.This suggests good performance in classifying the different categories.

Note:
The classification appears to have been highly accurate, with user's accuracy rates ranging from 81.11% to 100%.The 'Dense forest' and 'Open forest' classes were particularly well identified.

Dense forest class (DF)
We established coverage levels at or above 40%, encompassing densities spanning from 100 to 220 trees per hectare (ha), with an average density of 150 trees/ha.Concurrently, the average coverage stands at 55.25% (Fig. 8 and 9).A closer examination of Fig. 8 and 9 reveals patterns in the coverage and density data, such as the presence of clusters of high-density areas and the relationship  between coverage and density.These patterns may be indicative of underlying ecological processes, such as the influence of climate change, deforestation, and afforestation activities on forest dynamics (Verma et al., 2013).

Open forest class
It was agreed for this class to have a coverage ranging from 10% to 40%.The average density is 69 trees/ha, and the average coverage is 22.52% (Fig. 10 and 11).The search results provide information on the classification of forests based on their density and coverage, as well as on the methods for measuring forest cover density.The sources indicate that the "Open forest" class has a coverage ranging from 10% to 40%, while the "Dense forest" class has a coverage greater than 70%.Those results suggest that forest cover density can vary significantly depending on factors such as geography, climate, topography, and forest management practice (Verma et al., 2013).

Bare soil class
The term "bare soil" does not exclude the presence of subjects such as thuriferous juniper and other secondary species, generally pioneer species indicating significant degradation.As a result, the coverage does not exceed 10%, and the average density is 14 trees /ha (Fig. 12 and 13).
That being said, a detailed understanding of land cover provides more meaningful insights than analyzing it solely through statistics.Density alone does not convey the state of a population, but when combined with coverage, the insights gained better express the pressure exerted and the dynamics of the vegetation.Furthermore, a study of the structure and architecture of this forest is unquestionably warranted.

Spectral signatures and spectral distances
This stage represents a refinement in the classification process, essential to ensure the practical reliability of data for scrutinizing thuriferous juniper forest dynamics.The outcomes validated a spectral signature aligned with the expected reflective properties of the species, consistent with established scientific literature in this domain.Visual representations in Figure 14 delineate the comparative analysis of spectral signatures and spectral distances among the designated classes

Spatial evolution of El Houanet's thuriferous juniper forest between 1986 and 2022
This essentially involves comparing the relative areas of different cover density zones between the years 1986, 1994, 2003, and 2022, studying and evaluating the variations (Fig. 15).Since we do not have ground truth for the reference years, we estimated it based on the 1986 and 2003 images.The basis of this diachronic analysis relies on the preliminary assumption that the respective cover density percentages of the three zones are the same across the five dates.Given that the five corresponding maps are based on similar documents, both in terms of acquisition method and the nature of the information they contain, we can consider a pixel-by-pixel comparison.
Indeed, we will focus on the three major cartographic units taken as a whole, as well as the evolution of the cover in particularly well-identified sectors.The comparison of the total area of the thuriferous juniper forest between the five dates highlights a significant dynamic of the population, combining both growth and regression, requiring a cross-analysis of the data from different perspectives.Initially, it is necessary to compare between two successive periods before conducting, in a second step, a comparison between the beginning and the end of the observation period.The intersection of two maps from two successive periods allowed the establishment of a third map depicting the evolution between these two periods.This is a diachronic map that synthesizes the recent spatial evolution of this forest over 36 years (Fig. 15).This overlay reveals a remarkable spatial change and reflects the following observations: • stable strata correspond to areas that have not experienced any change; • strata with evolutionary dynamics correspond to areas where vegetation has become denser; • strata with regressive dynamics correspond to areas where vegetation density has decreased.
Ground truth verification confirmed that changes in the density of these thuriferous individuals are primarily changes in the canopy structures of the trees, consequently affecting canopy coverage.It is evident that decline is irreversible, and this is evident in the strata that have transitioned to the "bare soil" stratum, indicating no change in density due to the loss of stumps.Indeed, we will refer to changes in coverage for all trends, except for decline, which will be characterized simply as pure regression.

Area evolution of el houanet's thuriferous juniper forest between 1986 and 2022
The statistical data applied in this step is based on estimates calculated in hectares from the change maps established for each period (Fig. 16).The diachronic analysis of the surface evolution of strata has allowed us to visualize the flow of changes from one stratum to another for each period A (1986-1995), B (1995-2004), C (2004-2013), and D (2013-2022).There are a total of nine flows, which represent combinations of three types of land occupancy density, given that this is a pure stand of thuriferous juniper.Table 8 represents the corresponding areas for each flow, along with the characterization of the dynamics and its percentage relative to the overall area of the stand.In general terms, this involves observing the changes in the areas of the three classes, abstracting from the conditions that dictated them.Table 9 and Figure 17 provide a total of these areas by class and by year of observation.The various fluxes illustrate the dynamics of this population.To better interpret this data, we can summarize it according to three trends: evolutionary, regressive, and then static, as depicted in Figure 18.As a summary, the evolution of the different strata within this Juniperus thurifera formation during the studied period was characterized by the following notable facts: • the coverage experienced a significant densification across the entire formation over 36 years.
• there was a significant resurgence of individuals in the stratum with very low density belonging to the land occupancy class known as "bare soil".This percentage decreased from 42% in 1986 to 16.62% in 1995.Then, there was a slight regression without returning to the initial year (27.3% in 2004; 30% in 2013 and 30.68% in 2022).• the clear and dense forest strata exhibited a similar dynamic:  − dense populations went from 9.30% in 1986 to 16.45%, 13.58%, 12.95%, and 14.73% respectively in 1995, 2004, 2013, and 2022.
More precisely, the establishment of evolution maps based on the cross-referencing of maps from different observations interspersed by nine years allowed for the distinction of two trends spanning eighteen years each: 1985-2004 and 2004-2022.These are diachronic maps that synthesize the recent spatial evolution of these Juniperus thurifera populations over these two consecutives 18-year periods.This superimposition reveals a remarkable spatial change and reflects the following observations for the period 1985-2004: • the trend in vegetation dynamics is clearly progressive; • a significant proportion of the central and northern parts of the plateau, with altitudes between 2000 and 2200 meters and characterized by very low densities or vast areas devoid of vegetation, have seen their density evolve towards clear populations; • dense populations, confined between 1800 and 2000 meters in altitude on the banks of watercourses, notably the Attache wadi to the east and the tributaries of the Al Abid wadi to the north, have expanded through densification of adjacent clear populations; • the areas that have remained stable in dynamics are the most predominant.
However, for the second period (2004-2021), the spatiotemporal evolution revealed a rather regressive trend.The overlay of the two maps highlighted the following observations: • a decrease in density of high-density formations, resulting in an increase in low-density areas; • a shift in areas devoid of vegetation due to a loss in density in the central plateau, accompanied by the appearance of Juniperus thurifera on the northern cliffs.

DISCUSSIONS
The ensuing discussion rigorously scrutinizes the intricacies of the studied forest population's evolution, encompassing the timeline from 1986 to 2022.This comprehensive analysis centers on the pivotal factors outlined below, each playing a defining role in delineating the observed dynamics.An effective approach to analyze the dynamics gleaned from prior studies could involve employing a blend of statistical tools and remote sensing techniques within Geographic Information Systems (Lausier & Shaleen, 2018).This method might entail statistically analyzing hydrological conditions within the study area, leveraging precipitation data, vegetation indices, and surface temperatures extracted from satellite imagery (Bevacqua et al. 2022).Subsequently, establishing correlations between these hydrological conditions and the dynamics of forest populations would offer valuable insights into the pronounced impact of climatic variability on forest health and evolution (Speich et al., 2019).

Main evolutionary factors
Through this diachronic analysis, we have endeavored to evaluate the recent evolution of the El Houanet thuriferous forest, uncovering discernible trends influenced by a confluence of anthropogenic and natural factors, and at times, an intertwined interplay between the two.It's imperative to underscore that our study does not aim to exhaustively catalogue all contributors to this regression but rather to identify and spotlight the most prominent factors.Notably, the primary influencers significantly impacting the vegetation dynamics in this locale are linked to:

Climatic factors
The analysis of linear regression between this vegetation dynamic and the evolution of rainfall quantities from 1986 to 2022 allows us to understand that rainfall is not directly responsible for the vegetation dynamic.However, indirectly, it is attributed to the cumulative effects of recurrent droughts observed in Morocco since the 1980s.Indeed, previous research (Laouina, 2004;Aderghal, 2011) explains this vegetation dynamic by two factors: the removal of woody material and increased pastoral pressure in forested areas as a result of droughts.This observation remains concerning, especially considering the significant decline in rainfall.According to the results obtained, the relationship between NDVI dynamics and climatic parameters (Precipitation, temperature) processed through successive decade averages is readily apparent for all study years where the correlation is positive or negative and significant at the 5% error level (the significance threshold being r=0.16).Nevertheless, one can distinctly note: • the expansion of BS area is negatively correlated with precipitation (CC = -0.54);• the enhancement of OF area is positively correlated with average precipitation (CC = 0.59); • and quite remarkably, the improvement of DF area is positively correlated with the increase in average minimum temperature, and slightly positively correlated with average precipitation (CC = 0.50).
As confirmed by Thi et al. (2023), a robust relationship is evident between NDVI dynamics and climatic parameters (precipitation, temperature) analyzed across successive decades for all study years, showing statistically significant correlations at a notable level of significance.

Spatiotemporal dynamics of precipitation
Morocco, in general, has experienced a decrease in precipitation since the mid-1970s and early 1980s (Bijaber et al., 2018).The time frame for analyzing the temporal dynamics of precipitation (1986-2022) falls within a period where rainfall quantities have shown significant fluctuations between declines and increases.Moreover, this new analysis covering a span of 36 years enables the examination of interannual variations in precipitation to highlight various trends: • Near-zero or highly positive rates of change in evolution.• Notable increase in precipitation.
• Decreasing trend in precipitation quantities.
To verify the presence or absence of trends in these precipitation series, the Mann-Kendall test is applied.This is a non-parametric trend test, initially studied by Mann (1945) and later refined by Kendall (1975) and improved by Hirsch (1982Hirsch ( , 1984) ) to account for a seasonal component (Yue and Pilon, 2004).The Kendall's S statistic represents a balance, where a positive value indicates an upward trend in the studied phenomenon, and a negative value indicates a decrease over time.However, to determine if these increases and decreases are significant enough to be considered a meaningful trend, the p-value of the Mann-Kendall test must be lower than the predetermined level of significance (which is set at 5% here).Thus, the values of the Mann-Kendall test (S) regarding the search for significant trends in these time series confirm this general observation of interannual precipitation evolution.Furthermore, the observation of the coefficients of variation helps to understand that the trends in precipitation, whether increasing or decreasing, are not linear; it stands at 30%.The deviation from the annual mean is 152 mm.This, in fact, reflects the complexity in modeling this interannual precipitation dynamic (Fig 19) In general, the interannual dynamic of precipitation shows a moderate upward trend, summarizing an initial phase of increase and decrease from 1986 to 1988.The segments from 1988 to 1996 indicate an upward trend in precipitation.This is followed by a slight decrease until the end of the series.The year 1996 appears as the rainiest year of the entire series.As for the minimum, it was recorded in 2017.However, in comparison to the lower limit of the chronological series of vegetation dynamics (1986), the least humid year is 2017.Moreover, it is worth noting the absence of a series breakpoint from 1986 to 2022, according to the Pettitt test with 99% confidence intervals (Fig. 20).

Seasonal evolution of precipitation
In general, the study area is subject to a semiarid Mediterranean climate regime.This precipitation pattern is unimodal, with a single wet season from September to May.The average annual precipitation is around 500 mm.The spatial distribution of precipitation exhibits a North-South and East-West gradient, with the northern and western areas of the zone being the wettest.The significant interannual variation in precipitation is also indicative of marked intra-annual variability.The wettest month of the year is November, while July is the driest.In typically dry months, the amount of precipitation is usually minimal.
Although the Mann-Kendall test (Fig. 21 and Fig. 22) rejected the hypothesis of the existence of significantly remarkable interannual trends in monthly precipitation, with a confidence level of 95% for all months except for August, we can still discern slight evolving trends which can be categorized into classes, distinguishing between wet and dry months: • Among the wet months, a decrease was observed in February, April, May, and December; while January, March, October, and November showed an increase.
• As for the dry months, a decrease was recorded for June.However, July and August experienced an increase, with the latter being statistically significant according to the Mann-Kendall test for August.Prior research supports the significant variability in precipitation, both on a monthly and seasonal basis, over similar time frames.While the observed trends in monthly precipitation aren't globally significant, these studies also highlighted notable variations in specific wet and dry months (Zhang et al., 2021).Moreover, additional studies have affirmed seasonal regression dynamics, particularly during the winter season, traditionally associated with the highest precipitation levels of the year (Anderson et al. 2015).These observations substantiate the possibility that changes in precipitation patterns could potentially influence the observed developments in forest populations, especially during the wet seasons, where significant regressions have been identified (Moraglia et al. 2022).Building upon these established patterns, the examination of seasonal precipitation trends over this period fails to reveal definitive patterns.However, noteworthy observations emerge: autumn, spring, and summer portray minimal positive slopes, while winter, recognized as the wettest season, consistently exhibits a negative slope (Table 10, Fig. 23).Those slopes indicate slight yet consistent variations in seasonal precipitation, with slight positive trends for autumn, spring, The meticulous study of spatiotemporal dynamics of precipitation provides a window into climate variability.These analyses unveil seasonal and monthly nuances, significant trends, and nonlinear interannual evolutions, highlighting the inherent complexity in modeling these climatic patterns.However, this exploration of precipitation extends beyond a mere assessment of rainfall trends.It opens the door to a deeper reflection on annual humidity dynamics.By leveraging indices such as the Standardized Precipitation Index (SPI), we delve into a more holistic analysis of annual humidity.This reveals both dry and wet years, offering a comprehensive insight into the evolution of humidity over a 36-year period.

General characterization of annual humidity dynamics
The period from 1986 to 2022 experienced a deficit in rainfall compared to the previous decades.The identification of whether a year is wet or dry can be determined through the standardized precipitation index (SPI) (Ali et al., 2008).The SPI was calculated for the period from 1986 to 2022.This interval provides an insight into the temporal evolution of humidity in the area.Thus, over this 36-year period, the study area experienced 11 wet years (SPI > 0.5), which accounts for 30% (Fig. 24).Although some years showed a decrease in precipitation, the SPI indicates an equal number of genuinely dry years (SPI < -0.5), also totaling 11 years, representing 30% of the period.However, the majority of these wet years (8 out of 11) occurred after 1995, aligning with the general trend in precipitation evolution.Years with a deficit in rainfall are more significant (23 out of 36).The overall trend of the SPI suggests an area trending towards rainfall deficits, contrary to the general observation of precipitation dynamics.Nevertheless, the analysis of the interannual dynamics of rainfall in this area reveals that this humid characteristic is increasingly threatened, given the significant decreases in precipitation.
The synthesis of seasonal, monthly, and SPI analyses collectively confirms a notable shift in humidity patterns between 1986 and 2022.Seasonal trends showcased slight variations, displaying consistent slopes across different seasons, while monthly analyses underscored intermittent fluctuations in specific months.This complements the SPI findings, solidifying a noticeable change, highlighting both wet and genuinely dry years over the 36-year period.These collective observations underscore a significant deviation in precipitation dynamics, indicating an emerging trend towards decreased precipitation and an increasingly threatened humid characteristic within the region.This concerning shift in precipitation may further influence the unstable dynamics of forest populations, leaning towards regression, coinciding with the findings of (Chen & Ford, 2023)water resources, and natural ecosystems.Recently, there is growing attention to the transitions of precipitation extremes, or shifts between heavy precipitation

Evolution of vegetation hydric conditions: 1986-2022
The vegetation hydric conditions from 1986 to 2022 are analyzed through the vegetation condition index (VCI) and surface temperatures.The VCI is an index strongly linked to meteorological conditions, particularly precipitation.Expressed as a percentage, it indicates various degrees of drought severity (Ramesh et al., 2003;Kogan 1997Kogan , 2003)).The analysis of VCI variation reveals low levels of vegetation conditions (generally below 30% on average) within the study area (Fig. 25).
Although forested areas are known for maintaining better vegetation conditions (averaging between 50 and 70%), the interannual evolution of the VCI indicates very strong variations in these vegetation conditions from one year to The observed trend of intensified severe droughts, as indicated by the analysis of the VCI, poses a notable threat to the resilience of thuriferous juniper formations.Despite the historical adaptability of this species to fluctuating climatic conditions, the recurrent periods of extended and severe droughts throughout the examined timeframe (1986-2022) present a significant concern.The inherent resilience of juniper forests faces gradual challenges under sustained adverse conditions, particularly the increased frequency of severe droughts observed in the latter decades.These climatic shifts, resonating with findings from related studies, accentuate the vulnerability of thuriferous juniper populations, potentially impeding their continued vitality and progression (Tallieu et al., 2023).
Integrating the nuanced understanding of spatiotemporal precipitation dynamics, the profound insights into annual humidity shifts, and the detailed analysis of vegetation hydric conditions converge to present a holistic portrayal of the ecosystem's response to climatic variations.These intertwined facets, examined collectively, serve as a robust foundation to comprehend the intricate interplay between climate fluctuations and ecosystem resilience.The amalgamation of broad climatic trends with specific vegetation responses unveils a comprehensive narrative, emphasizing the vulnerabilities and adaptations within these forest ecosystems.

Human activities
The forest formations in the El Houanet plateau have undergone significant changes between 1986 and 2022, influenced by anthropogenic, natural, and interactive factors.The primary disturbances caused by human activities have had a significant impact on vegetation dynamics in this region.Anthropogenic pressure has been intensified by sedentarization, leading to the emergence of new forms of forest resource exploitation to meet the growing needs of neighboring populations.Human activities are the main drivers of the regression of the El Houanet juniper forest, particularly through grazing and overgrazing.

Grazing and overgrazing
The forest pastures on the El Houanet plateau are subject to heavy pastoral pressure due to the presence of numerous herds that find the necessary forage resources there.However, this high pastoral load leads to an estimated overgrazing rate of 83%, directly impacting the degradation of forest pastures by reducing their standing biomass and affecting their structure.The practices of trimming and pruning frequently used by herders to feed their livestock, coupled with significant herd mobility, also result in advanced mutilation of the Juniperus thurifera stands, affecting the dynamics of vegetation formations in general, and the natural regeneration of the thuriferous species in particular.Furthermore, users have made significant transformations in the management of the system.The management of pastures no longer adheres to customary laws, leading to individualized space management.This increased weakness in customary laws is the cause of overgrazing and, consequently, the degradation of vegetation cover and soil.The agricultural production system, through the cultivation of collective pasture lands, further increases the pressure of livestock on the permanent pastures represented by forest formations (Abaab, 1995;Bourbouze, 2000;Benbrahim, 2004).

Illegal firewood harvesting
In the El Houanet region, illegal firewood harvesting has a considerable impact on forest resources.Although pastoral practice is common, unauthorized wood cutting accounts for nearly 82% of recorded forest offenses between 1990 and 2019, according to the reports of the Tagueleft Water and Forestry Services.Firewood collection is an essential daily activity for households in the region, as there is a significant energy need, especially in winter.The use of butane gas is not sufficient to meet the energy needs of the rural population, especially when there is snowfall contributing to the isolation of certain areas.Collecting deadwood is not enough, leading to the cutting of living trees in the forest.Rural families are supposed to spend an average of 3,000 to 4,500 DH per year to buy firewood legally, but this situation is often inaccessible or unbearable due to their limited means, forcing them to turn to the black market where the purchasing cost is much lower.Moreover, the illegal sale of firewood creates an income-generating activity, especially for the young people of the region during winter (Moufaddal, 2008).

CONCLUSION
Our comprehensive analysis spanning from 1986 to 2022 has unveiled intricate insights into the evolution of the Juniperus thurifera forest formations on the El Houanet plateau within the Boutferda BEIS.Employing a diachronic approach complemented by rigorous fieldwork, this study delves deep into the complex interplay of factors shaping these ecosystems.
The integration of LANDSAT 5 and 8 satellite imagery has granted us a meticulous spatial-temporal understanding, revealing a mosaic of change across the landscape.These visual representations, accompanied by quantified surface area changes, vividly illustrate both regressive and progressive trends in forest coverage.Notably, these shifts permeate all identified strata, signaling a profound impact on the entire forest ecosystem.Nestled in the central High Atlas, the Boutferda BEIS.epitomizes the oromediterranean stage, characterized by the resilient presence of Juniperus thurifera.A precise altitudinal floristic transition underscores its unique character, culminating in a pre-steppe formation.Contrary to prior characterizations, our climatic scrutiny reveals a regime marked by pronounced precipitation ranges and significant thermal amplitude.
The meticulous dissection of spatial evolution between 1986 and 2022 unravels a dynamic narrative.Our methodology, anchored in diachronic mapping, illuminates areas of stability, evolutionary vigor, and regressive tendencies.Ground truth verification reaffirms the pivotal role of canopy structure alterations in reshaping canopy coverage density.As we aimed to identify primary evolutionary factors, our study elucidated the multifaceted interplay between anthropogenic and natural influences.Among these, intricate precipitation dynamics oscillate in a non-linear fashion, highlighting the complexity in modeling interannual precipitation patterns and the need for a comprehensive approach.Human activities, particularly grazing and overgrazing, emerge as pivotal agents in the regression of the El Houanet juniper forest.The cumulative impacts of pastoral pressures, alongside emerging exploitation paradigms and illegal firewood harvesting, underscore the urgent need for sustainable management strategies.
In summation, our longitudinal investigation provides a comprehensive account of the juniper thuriferous forest formations' evolution and the manifold factors influencing their transformation.The insights gleaned hold salient implications for conserving and sustainably managing this vital ecosystem within the broader context of ecological preservation.

Fig. 1 .
Fig. 1.Location of the Boutferda biological and ecological interest site

Fig. 14 .
Fig. 14.Comparison of spectral signatures of selected classes

Fig. 15 .
Fig. 15.Maps depicting the spatial distribution of strata over time

Fig. 16 .
Fig. 16.Maps depicting the flow of changes between strata over time, land use change map of Boutfera BEIS during the period 1986-2022

Fig. 17 .Fig. 18 .
Fig. 17.Evolution of areas between 1986 and 2022 . The phases of this dynamic can be summarized as follows: • a first phase of both increases and decrease in interannual precipitation from 1986 to 1988; • the segments from 1988 to 1996 show an increasing trend; • the segments from 1996 to 2001 exhibit a decreasing trend; • the segments from 2001 to 2010 demonstrate an increasing trend; • the segments from 2010 to 2017 indicate a decreasing trend; • the segments from 2017 to 2022 show an increase in 2018, followed by a decrease in 2019, and then a rise until 2022.

Fig. 21 .
Fig. 21.Summary of the Mann-Kendall trend test and p-value diagram

Fig. 22b .
Fig. 22b.Interannual evolution of monthly precipitations in the study area

Fig. 23 .
Fig. 23.Interannual evolution of seasonal precipitations in the study area

Fig. 24 .
Fig. 24.Characterization of the interannual evolution of precipitation through the standardized precipitation index (SPI)

Table 1 .
Emberger-sauvage pluviothermal quotient in the study area

Table 2 .
Annual growth rates

Table 3 .
Confusion matrix for the year 1986 (Landsat 5 TM Sensor) Dense forest -the classification correctly identified dense forest with a user's accuracy of 100%, this indicates that all regions identified as dense forest are indeed so; open forest -similarly, open forest was correctly identified with an accuracy rate of 99.80%; bare soil -bare soil was also well identified with an accuracy rate of 98.53%. Note:

Table 6 .
Confusion Matrix for the Year 2013 (Landsat 8 OLI Sensor)

Table 7 .
Confusion Matrix for the Year 2022 (Capteur Landsat 8 OLI) Note: User's accuracy rates are high, ranging from 87.01% to 99.81%.Overall, the classification results seem to indicate a high performance of the classification model, with precision rates (producer's accuracy) and overall accuracy (overall accuracy) well above 80%.This suggests that the model has been effective in classifying the images into different categories such as dense forest, open forest, and bare soil for the various years and sensors considered.

Table 8 .
Changes in area flux for strata by observation period Note: DF -dense forestt, OF -open forest, BS -bare soil.

Table 9 .
Areas of strata by year of observation

Table 10 .
Summary of the Mann-Kendall Test