Intensified and Extended Growing Seasons in Abies marocana Forests (2000–2024): A Robust Seasonal Trend Analysis Using 16-Day MODIS EVI Time Series

We modelled, for the first time, the seasonal dynamics and long-term trends of Abies marocana forests (Rif Mountains, northern Morocco) using remote-sensing-derived vegetation indices. Using the MODIS Terra Vegetation Indices product MOD13Q1 (enhanced vegetation index, EVI; 16-day frequency; 250 m spatial resolution) from 2000 to 2024 (575 images over 25 years), we applied a robust seasonal trend analysis (RSTA) workflow, representing an inferential extension of classical seasonal trend analysis (STA) through the explicit control of Type I error under serial and spatial correlation. This approach combined: (i) harmonic regression to capture the annual and semi-annual cycles of A. marocana forests, estimating seasonal amplitudes and phases while filtering out low-frequency noise; (ii) an iterative trend-free prewhitening (TFPW) procedure following Wang and Swail, applied only to time series with significant serial autocorrelation according to the Durbin–Watson test; (iii) the Theil–Sen slope (TS) estimator, a robust non-parametric method, to quantify the magnitude and direction of seasonality trends; (iv) the contextual Mann–Kendall (CMK) test to assess the statistical significance of seasonality trends, while correcting for spatial autocorrelation and accounting for cross-correlation among neighbouring pixels; (v) the Benjamini–Hochberg (BH) procedure to control the false discovery rate (FDR), ensuring that only statistically robust seasonality trends were retained; and (vi) reconstruction of seasonal curves representing the beginning and end of the study period and derivation of phenological metrics from the statistically significant seasonal trends retained after inferential filtering. After applying the complete analytical workflow, statistically significant trends were detected in 79.2% of pixels within A. marocana forests, compared with 86.4% when prewhitening and false discovery rate control were not applied. All Theil–Sen slopes retained by the RSTA workflow were positive, with a mean slope of approximately 0.00175 EVI year−1, corresponding to an average annual increase of roughly 0.7% and an overall increase of approximately 15% over the 2000–2024 study period relative to the initial mean EVI conditions. Browning trends identified by classical STA were not supported after inferential filtering and FDR control, indicating that all these patterns were spurious or only marginal, and confined to limited areas and edge zones. The reconstructed seasonal trend curves were consistent with a longer growing season, although this inference is based on land-surface vegetation dynamics rather than direct phenological observations. The long-term ecological consequences of these changes in seasonal vegetation activity will hinge on the interactions among warming, rising water demand, and potential disturbance regimes under future climatic conditions.