## 2 Materials and Methods
### Data Analysis  
#### 2.6.1. Phenological Analysis of RS Variables  
We utilized linear mixed models (LMMs) to capture the spatial and temporal variations of spectral indices across the eight levels of forest conservation (Table 1). The LMM approach was chosen for its ability to explicitly handle both fixed and random effects, making it well-suited for addressing temporal heterogeneity across semi-arid seasonal ecosystems. Specifically, to account for spatio-temporal heterogeneity in LMMs, we included forest type as a fixed effect, as it represents the primary factor influencing spectral variation and conservation status, while transect ID and month were included as random effects, accounting for variability at different spatial locations and across time. The interaction term (Transect ID \(\times\) Month) was introduced to account for site-specific seasonal patterns, ensuring that temporal changes were analyzed within each site rather than being treated as uniform across the study area. This allowed us to explicitly model seasonal variation at each site, ensuring that different forests could exhibit distinct temporal dynamics.  
As for seasonal and geospatial differences, the random effect for Transect ID captures the inherent spatial variability among sampling sites, while the random effect for Month accounts for differences across time points, allowing us to model phenological changes in spectral indices. The inclusion of these random effects improves model transparency and prevents pseudo-replication, ensuring that repeated measurements from the same sites are not treated as independent observations. The denominator degrees of freedom were calculated using the Kenward-Roger approach. The intra-class correlation coefficient (ICC) was presented as an indicator of the extent of random variance in the models. The quality of fit for mixed models was assessed using conditional and marginal determination coefficients (R\({}^{2}\)c and R\({}^{2}\)m) to quantify unbiased assessments of variance explained by fixed effects alone and by fixed plus random effects, respectively [70].  
After analyzing the models, we computed the marginal mean and the accompanying 95% confidence interval for each spectral index. We then conducted a Tukey test with Bonferroni correction to determine whether there were statistically significant differences between the forest classes. We used the R programming language (version 4.3.2 [71]) for all analyses, implementing the following packages: 'ggplot2' (version 3.5.1) [72] for data visualization, 'lme4' [73] for mixed-effects modeling, 'rstatix' (version 0.7.2) [74] for  
\begin{table}
\begin{tabular}{c c c c c} \hline
**Acronym** & **Name** & **Formula** & **Proxy** & **Reference** \\ \hline \multirow{2}{*}{MCARIS\({}_{\text{cert}}\)} & Modified Chlorophyll & \multirow{2}{*}{\(((REDEDGE-RED)-0.2(REDEDGE-GREEN))\)} & \multirow{2}{*}{Leaf chlorophyll content} & \multirow{2}{*}{[67]} \\  & Absorption in Reflectance & & & \\  & Index & & & \\ \hline \multirow{2}{*}{NDWI\({}_{\text{Cao}}\)} & Normalized Difference & \multirow{2}{*}{\(\frac{\text{NIR-SWIB}}{\text{NIR+SWIR}}\)} & \multirow{2}{*}{Leaf water content} & \multirow{2}{*}{[68]} \\  & Water Index & & & \\ \hline BI\({}_{2}\) & Brightness Index 2 & \(\sqrt{\frac{\text{RED}^{2}+\text{GREEN}^{3}+\text{NIR}^{2}}{3}}\) & Presence of bare surfaces & [69] \\ \hline \end{tabular}
\end{table}
Table 3: Spectral indices describing forest characteristics derived from Sentinel-2 imagery, retained for further analysis due to low multicollinearity (VIF \(<\) 5).  
statistical tests,'multcomp' (version 1.4-28) [75] for multiple comparisons,'sjPlot' (version 2.8.17) [76] for visualizing model outputs, and 'usdm' (version 2.1-7) [77] for addressing multicolllinearity in ecological data.