Contrasting ecosystem vegetation response in global drylands under drying and wetting conditions

2.1 Aridity

In this study, aridity is defined as 1 − AI. The AI depends on P and PET and thus the index itself, as well as AI-based projections of dryland extent and expansion, depend on the input data, particularly the estimates of PET (Lian et al., 2021). We use P and PET data from ERA5-Land reanalysis (Table S1) to calculate aridity. P is provided hourly/monthly at 0.1° spatial resolution, and the PET is a recently developed global hourly product (hPET) (Singer et al., 2021) at 0.1° spatial resolution. (Singer et al., 2021). We resampled P and PET to 0.25° to match the most common resolution of all datasets used within the study, calculated annual averages from 2000 to 2020 and the aridity was calculated subsequently.

2.2 Study area

Our analysis focuses on dryland ecosystems, specifically arid, semiarid, and dry subhumid areas with aridity between 0.35 and 0.95 (Cherlet et al., 2018). Polar and cold areas with PET <400 mm are not included in the definition of a dryland (Cherlet et al., 2018) and thus, we exclude polar areas as well as areas classified as having dry and cold summers and winters (classes 19, 20, 23, and 24) based on the Köppen–Geiger climate classification (Beck et al., 2018). We further masked out irrigated cropland and wetland based on the European Space Agency (ESA) Climate Change Initiative (CCI) land cover data (classes 20, 150, 160, and 170) (Table S1).

2.3 Temporal cluster analysis

First, the aridity time series was normalized by the mean and standard deviation and a 5-year moving average (Jiao et al., 2021) was applied to remove the interannual fluctuations. Next, we used an agglomerative hierarchical cluster analysis to identify dominant temporal oscillations in the aridity time series and map their spatial distribution (Fraley & Raftery, 2002; Johnson, 1967; Radovanović et al., 2020; Wu et al., 2009). To measure the similarity between two time series, we use dynamic time warping (DTW). Although this method is more complex than others (e.g., simple Euclidean distance), it accounts for warping or delays in the time series present in the data. DTW thus seeks for the temporal alignment of two time series by minimizing the Euclidean distance between the two aligned time series and stores the output (Xi et al., 2006) in a distance matrix (Berndt & Clifford, 1994; Jeong et al., 2011). In agglomerative clustering, each cluster (i.e., one pixels' aridity time series) is initialized as its own individual cluster. The algorithm then iteratively merges the least dissimilar clusters (i.e., another pixels' aridity time series) into larger clusters, relying on the distance matrix. We use complete linkage to measure dissimilarity, which is the maximum distance between observations across all pairs of observations in two clusters. As a result, complete linkage tends to generate clusters with similar members. Hierarchical clustering does not require a specific number of clusters prior to the analysis, instead the final clusters are visually identified afterwards using a dendrogram. We also correlated the per-pixel time series of aridity to the mean time series of the cluster it belongs to. From this, the R² and p-values provide a measure of probability of the belonging (and uncertainty) of each pixel to its assigned cluster. We calculated the linear trend based on the Theil–Sen estimator on the mean aridity time series for each cluster to quantify the average aridity increase/decrease per cluster.

2.4 Essential ecosystem state and climate variables

We selected nine ECVs that describe components of the land (vegetation cover, aboveground biomass, surface reflectivity, soil water availability, and land cover and burned area) and atmosphere (wind speed, cloud cover, temperature, and vapor-pressure deficit) available from satellite observations covering the full study period from 2000 to 2020. The data were downloaded at their original resolutions (Table 1) and were subsequently harmonized (spatially aggregated) to a spatial resolution of 0.25°. The temporal resolution of the data varies between daily and 10-daily, but were averaged into annual observations. Additionally, to include a metric for vegetation functioning, we used the Sequential linear Regression Slope (SeRGS, Abel et al., 2019). SeRGS describes the vegetation rainfall relationship using the slope of an ordinary least square (OLS) regression for short sequential periods (temporally moving window). The trend over these slopes then represents the change in the dependent variable (vegetation productivity) per unit change in the independent variable. It can therefore be used as an indicator of vegetation sensitivity to changes in rainfall, a proxy for vegetation functioning in drylands (Abel et al., 2020), since vegetation productivity is tightly linked to the availability of water (mainly through rainfall). Some of the selected (atmospheric) variables are closely correlated with, or embedded in the calculation of, aridity. These variables will still be considered in the analysis to explore potential differences in their response to different dynamics of aridity. The focus of the analysis is, however, on the land-related variables.

We calculate the pixel-wise Pearson correlation between the time series of aridity and each ecosystem state variable. For each pixel, we sum the variables that have a significant (p ≤ .05) correlation, which provides an indication of the general relationship between the ecosystem state and aridity. We define three categories representing the number of significantly correlated variables: few (0–3), moderate (4–6), and majority (more than 7). For each category, we check which variables were the ones most often significantly correlated as well as the distribution of the aridity clusters within.

2.5 Ancillary data on land cover

Based on the ESA CCI land cover products (Table S1), we first merged individual classes of related types to the following overall classes: tree cover (CCI land cover classes: 50–90), mosaic tree cover (100–110), shrubland (120–122), grassland (130), sparse vegetation (150–153), cropland (10–40), and bare (including urban) areas (190–202). These classes were used throughout the analyses to potentially link observed patterns to specific land cover classes. We also estimated the change in land cover between 2000 and 2020 toward more or less green/natural land cover classes by considering changes toward the tree cover class being positive changes and changes toward bare soil/urban area being negative changes (Figure S2). This was used to create a dynamic land cover (change) variable that could then be related to changes in aridity over time.

2.6 Model to explain LAI trends under stress/non-stress conditions

Increasing VPD and radiation (less cloud cover) increase transpiration, which causes plants to close their stomata to minimize water loss, and ultimately lead to a reduction in photosynthesis (Tagesson et al., 2021). The probability of plants to frequently experience water-related stress (called stress hereafter) further increases with water deficit in the soil from reduced precipitation (Barkhordarian et al., 2019). As such, we identified areas with potentially stressed (non-stressed) vegetation as those where VPD increased (decreased) and soil moisture, precipitation, and total cloud cover decreased (increased). Pixels that do not fall within these classes were merged in a third class describing “other conditions.” We then limited the areas potentially experiencing stress (non-stress) to the drying (wetting) clusters and analyzed which factors affected positive and negative Theil–Sen trends in LAI in these areas. As an ECV, LAI is defined as a measure of foliage amount and known to partly control important mass and energy exchange processes, such as photosynthesis and rain interception, which are crucial to the ecosystem as they couple vegetation with the climate system (GCOS, 2022). LAI captures both structural and biophysical (e.g., foliage chlorophyll content) attributes of the vegetation and is used here as a proxy for vegetation productivity/growth, which is a good indicator for the general state of the ecosystem (Fang et al., 2019; Smith et al., 2019). To explain LAI trends, we fitted Gradient Boosted Decision Tree (GBDT) models (xgboost in Python) using anthropogenic, climatic, topographical, and soil variables as covariates. Those include organic carbon (OC) density in soil, root-zone (RZ) water storage capacity, canopy height, fire frequency, topographical attributes (land cover and elevation), human influence (human footprint and population density), and a set of comprehensive climate attributes (trend, average, and interannual variability for aridity, precipitation, temperature, soil moisture, and cloud cover) (see Table S1). Gradient boosting models build an additive—and in this case—regression model in a forward stagewise fashion, enabling the user to relate a response to its predictors. Thus, GBDT models are effective in modeling complex ecological systems because they accommodate continuous and categorical predictors, missing data, and outliers, do not require prior data transformation, and capture complex nonlinear relationships. In total, we constructed four GBDT models including: positive LAI trends under stress and non-stress as well as negative LAI trends under stress and non-stress conditions. We do not apply a variable selection procedure and instead use all available variables to parameterize the models. This will ensure models with the highest possible explanatory power, and overfitting is no concern, as our aim is to explain and not to predict LAI trends. The data were split into a training and test dataset (80% and 20% of all data, respectively) and we parameterized each GBDT model using an automated random search. The random search entails a threefold cross-validation, that is, the dataset was shuffled randomly and split into three groups. The random search exhaustively generates candidates from a range of given parameter values, and for all possible combinations of parameter values models are fit with the training data. Model performance for all parameter combinations is evaluated on the test datasets, and the best parameter combination is retained. To interpret the GBDT output, we used SHAP (SHapley Additive exPlanations) values (Lundberg & Lee, 2017). The SHAP value is the average marginal contribution of a feature value across all possible coalitions and is a way of interpreting machine learning models based on a concept originating from game theory (Lundberg & Lee, 2017). To get the global importance, absolute SHAP values per feature are averaged across the data.

3.1 Global spatiotemporal variability of aridity between 2000 and 2020

We determined five clusters with distinct spatial patterns—using the agglomerative hierarchical cluster analysis, which represent the dominant temporal oscillations over the past two decades (Figure 1a,b). Two clusters with increasing (drying) and two with decreasing (wetting) trends in aridity, as well as one cluster with no temporal trends were found. Overall, we observe that 31.6% of all pixels within the study area experienced wetting (clusters 1 and 2), 44.5% experienced drying (clusters 3 and 4), and 23.8% exhibited no temporal trends in aridity (cluster 5) between 2000 and 2020. Drylands are characterized by strong interannual variability in precipitation, which becomes an inherent characteristic of aridity and is visible in all clusters (Figure 1b). Average linear trends in aridity based on the smoothed moving average time series are significant for clusters 2, 3, and 4, but not for clusters 1 and 5. Cluster 1 represents 14.0% of the global dryland area (4.7% considering only pixels with a significant trend) and is characterized by a fast decrease in aridity until ca. 2010 from where the decrease stopped and slightly reversed. The average aridity in the cluster decreased by −0.0006 year⁻¹ but is not significant (p = .38). In cluster 2, aridity showed a stronger and significant (p = .0007) average decrease (−0.003 year⁻¹), which accelerated toward the end of the study period. Cluster 2 covers 17.6% of the global dryland area (24.4% considering only pixels with a significant trend). The largest significant increase in aridity was observed in cluster 3 (on average 0.004 year⁻¹, p = .000), which also represented the largest absolute change in aridity across all clusters. Here, aridity increased strongly in the first 4 years of the time series and thereafter it stabilized and increased moderately. Cluster 3 covers 26.2% of the global dryland area (55.2% considering only pixels with a significant trend). Cluster 4 characterizes areas with a relatively slower yet significant (p = .04) average increase in aridity (0.001 year⁻¹) and represents 18.3% of the global dryland area (8.7% considering only pixels with a significant trend). The strongest increase occurred again in the beginning of the study period, and we observe a distinctive drop in aridity around November 2010 and further recovery until 2015. To check for uncertainties in the classification of individual pixels into a certain cluster, we calculated the strength of association for a given pixel to actually belong to the cluster it has been classified in (Figure S1, see methods). The cluster with highest uncertainties, that is, pixels that were attributed to a specific cluster, but the pixels aridity was not significantly correlated with average aridity in this cluster, was cluster 5 representing no temporal trends in aridity, followed by clusters 1 and 4 representing moderate wetting and drying trends, respectively (Figure S1). Overall, spatial patterns and statistics were consistent regardless of weakly associated pixels being included or excluded (Figure S1).

Generally, we observed that aridity trends are strongest on the wetter side of the aridity gradient and the magnitude of the trends decreased along that gradient for all clusters in a similar manner (Figure 1c). Areas that experienced the strongest drying trends (cluster 3) are predominantly located in the southern hemisphere as well as in large parts of Kazakhstan. Strongest wetting trends were observed in the central parts of the Sahel as well as in central East Africa and South Asia. In most areas, a clear gradient of the clusters (from drying to wetting trends or vice versa) is apparent, for example on the Iberian Peninsula, southern Africa, and western Sahel. We found drying and wetting trends across all land cover classes and continents (Figure 1). However, cropland, grassland, and bare soil were more often present in the clusters that got wetter (or did not have a trend) over the past 20 years, whereas shrubland, tree cover, and sparse vegetation were more frequent in the clusters affected by drying trends (Figure 1d).

3.2 Ecosystem response to changes in aridity between 2000 and 2020

We calculated the Pearson correlations between the interannual variability of each ecosystem state variable (EV) and aridity, and subsequently summed the number of significantly correlated variables across all EVs at the pixel level. Among the EVs, there is only a low degree of correlation (Figure S3). Overall, we found a clear gradual pattern where the number of pixels in drying clusters increased along with the amount of EVs significantly correlated with aridity (Figure 2b). The reversed pattern was found for wetting clusters: These were mostly characterizing areas with only few significantly correlated EVs, and their proportion decreased with an increase in the number of significantly correlated EVs (Figure 2b). However, there are regional differences. Areas with most significantly correlated variables (7–11) are predominantly located in the southern hemisphere (Figure 2a). These areas are dominated by clusters 3 and 4, that is, they are mainly characterized by an increase in aridity in the past 20 years (Figure 2b). On the other hand, areas with fewer significantly correlated variables (0–6) are primarily located in the northern hemisphere (Figure 2a). These areas have a relatively even share of all clusters, that is, they show both positive and negative trends in aridity with the majority of the negative/wetting trends found here (Figure 2b).

While the number of EVs that are significantly correlated with aridity varies at global scale, soil moisture and cloud cover are always significantly and overall strongest correlated with aridity (median Pearson values around −0.8 and −0.65, respectively). This was to be expected as both variables are strongly correlated with the input parameters used to calculate aridity. Among the land variables, LAI and biomass often have a significant correlation with aridity, with median values around −0.64 and −0.60, respectively (Figure S4; Figure 3c). Overall, soil moisture, cloud cover, biomass, tree cover, non-tree vegetation and LAI are negatively correlated whereas VPD, wind speed, albedo, and vegetation functioning (as captured by SeRGS) are generally positively correlated with aridity (Figure 3a). Changes in land cover and burned area are positively correlated with the exception of clusters 1 and 4 as well as cluster 2, respectively; here these EVs are negatively correlated (median values) with aridity (Figure 3c). However, only very few pixels (below 3%) show significant correlations between burned area, changes in land cover, and vegetation functioning (Figure S4). Both atmospheric- and land-related EVs—except tree cover and non-tree vegetation—have weaker correlations with aridity in wetting clusters as compared to drying clusters (Figure 3). We found largest differences in median correlations for cloud cover and albedo.

3.3 Key variables related to LAI trends under conditions of stress and non-stress

We analyzed linear trends in LAI between 2000 and 2020 as a function of ecosystems experiencing stress or non-stress related to water availability (Figure 4a). As expected, most of the areas where plants are potentially stressed were found in drying clusters (73.44%), few in wetting clusters (7.68%), and the remaining were located in the cluster with no temporal trends (18.89%). Within the drying clusters, we found areas which exhibited both positive (67.78%) and negative (30.65%) trends in LAI under stress (Figure 4b). On the contrary, within the wetting clusters, we found that the majority of trends were positive (78.98% under no stress and 14.45% under stress, respectively) (Figure 4b). Within the cluster with no temporal trend, most trends were positive, the majority (71.09%) under stress and fewer (18.24%) under no stress conditions (Figure 4b).

Regardless of the trend in LAI, areas not characterized as being stressed are dominated by cropland, whereas shrubland dominates areas where the system is potentially stressed. Areas characterized by both stress and decreasing LAI trends have a large share of cropland as well. Independent of stress, LAI is mainly increasing in areas classified as grassland and sparse vegetation. There are no clear dominating land cover classes identified for areas of decreasing LAI, yet there is a larger share of tree cover as compared to that of areas of increasing LAI.

We fitted GBDT models to understand which anthropogenic, climatic, topographical, and soil variables are associated with positive and negative trends in LAI under stress and non-stress conditions. The models exhibited moderate to high accuracy on classifying either positive (66.2%) and negative (75.4%) trends in LAI for areas where plants are stressed, and high accuracy on positive LAI trends (77.9%) where plants were not stressed. Results of a model on decreasing LAI trends in cases where plants did not experience stress were not produced due to the small sample size of 82 pixels. Overall, we show a high relevance of the amount of cloud cover on LAI trends according to the SHAP values. Among the climatic variables, temperature (average and variability) was found to be the most important for LAI trends. We further observed that more anthropogenic variables were related to LAI trends when plants did not experience stress, whereas soil variables were more often of importance under non-stress conditions. Besides cloud cover, the main driving factors for the decreasing trend in LAI under stressed conditions (Figure 5a) are RZ water storage capacity, and annual average temperature. Meanwhile, the increasing trend in LAI (Figure 5b) is driven by canopy height and average burned area. On the other hand, the increasing trend in LAI under non-stressed conditions is driven primarily by elevation and organic carbon density (Figure 5c).

Depending on whether there is water-related stress (or not), the same variables may have different effects on positive LAI trends (Figure 6). For example, the ability of RZ water storage capacity and organic carbon (OC) density to increase LAI depends on whether vegetation is under stress. That is, in non-stressed systems, LAI trends keep increasing along the gradients of these variables, but under stress further increases in the magnitude of these variables have no effect on the increasing LAI trends (Figure 6a,b). For canopy height, we observe the opposite: LAI trends in stressed systems increase strongly and linearly with increasing canopy height, whereas in non-stressed systems canopy height only has a moderate positive impact on LAI trends (Figure 6c). Average population density impacts LAI trends negatively in non-stressed systems, while in stressed systems it has a positive impact on LAI trends (Figure 6i). Looking at the impact of land cover on LAI trends, tree cover, shrubland, and cropland emerge as diverging classes (Figure 6o). In shrubland areas a positive impact on LAI trends is observed when the system is stressed, however, the relationship is reversed (i.e., negative impact on trends in LAI in shrublands) when the system does not experience stress. The opposite pattern was observed for cropland, where LAI trends are negatively/positively impacted when the system is stressed/not stressed.

4.1 Spatial distribution of aridity clusters

Over the past two decades, we observed more drying (44.5%) than wetting (31.6%) trends. The southern hemisphere is dominated by drying trends whereas we found a more balanced distribution of drying, wetting, and no temporal trends in the northern hemisphere. This is in line with other studies on historical trends in precipitation and drought indices reporting drying trends (although for different periods) for South and East Africa (Dai, 2013; Huang et al., 2016), as well as a significant positive trend in the land area under drought, higher drought frequency, and amplitude in the southern hemisphere making droughts more variable (Damberg & AghaKouchak, 2014). Studies also found large parts of South America, northern Australia, and Central Asia were characterized by wetting trends over six decades spanning ~1948–2010 (Dai, 2013; Huang et al., 2016). We found that these trends have been reversed during the period from 2000 to 2020, and the areas in question are characterized by some of the strongest drying trends (Figure 1), which is possibly due to accelerated warming in the past 20 years (Daramola & Xu, 2022; Maharana et al., 2021). Such phenomena are not uncommon in climate time series data and referred to as trend turning (Alexander et al., 2006; Zuo et al., 2019). We compared aridity trends during the period 1980–2000 to the trends during 2000–2020 (Figure S6) to indicate the extent to which trends in preceding decades relate to observed patterns during the study period. We found a reversal in 13.57% of the trends, which spatially match the identified clusters. This suggests that the majority of pixels are characterized by no change in trend direction or a change from and to a nonsignificant trend (86.42%, see Figure S6). The spatial patterns reported here are thereby related to past conditions only to a limited extent. The aridity clusters also show clear regional gradients of the aridity trends from dryer to wetter and vice versa; for example, in the Iberian Peninsula, southern Africa, and western Sahel. These patterns reflect the spatial distribution of rainfall that are possibly due to large-scale climatic factors such as Mediterranean Oscillation and the Angola Low (Müller-Plath et al., 2022; Munday & Washington, 2017; Rodrigues et al., 2021). Despite notable interannual variability of aridity, the smoothed time series have clear trend directions as the wetting clusters show a decrease in aridity and the drying cluster shows an increase (Figure 1b). Cluster 4 shows a distinct drop in aridity around 2010/2011, and further inspection of its spatial coverage shows that it is the dominant cluster in eastern Australia and northeastern Brazil. It is therefore likely that the drop in aridity is related to the El Niño-Southern Oscillation (ENSO) (Figure S5) that has caused heavy rains and flooding in those regions during 2010/11.

4.2 Relationship between aridity and ecosystem state variables

Generally, we found stronger links between aridity and EVs in clusters with increasing trends of aridity (Figure 2). This is somewhat expected because shifts in ecosystem functioning due to persistent reduction in water availability have been shown to result in a systemic acclimatization to water stress or a change in plant composition (Berdugo, Vidiella, et al., 2022). Aridity showed a strong positive correlation with VPD and a coinciding strong negative correlation with soil moisture (Figure 3), while the negative correlation between aridity and biomass and non-tree vegetation increased from wet to dry clusters. This pattern is expected because in drylands the proportion of non-tree vegetation such as grasses and shrubs are large relative to forest (D'Odorico et al., 2019) and have been shown to be sensitive to both increases in VPD (Konings et al., 2017) and reductions in soil moisture (Liu et al., 2020). The strong positive correlation that aridity has with albedo (especially in drying clusters) could be linked to two possible scenarios: (1) change in vegetative composition due to the combined increase in VPD and decrease in soil moisture, particularly in more humid environments (Lian et al., 2021), and (2) a persistent increase in the fraction of bare ground due to shrubification, which reduces the herbaceous cover (Schreiner-McGraw et al., 2020). It is interesting to note that burned areas did not exhibit any strong correlation with aridity as fire activity is enhanced in xeric drylands due to the prevalence of non-tree vegetation (Abatzoglou et al., 2018). Similarly, changes in land cover as well as vegetation functioning (as captured by SeRGS) did not show strong correlations with aridity, suggesting other drivers of their respective dynamics, for example, population pressure or land abandonment (Abel et al., 2020; Bernardino et al., 2020). In general, however, the limited number of pixels that showed a significant correlation between vegetation functioning and aridity, were positively correlated, demonstrating that vegetation functioning is increasing with increasing aridity and vice versa.

4.3 Understanding the key variables related to LAI trends under stress and non-stress

The GBDT models that were used to identify key variables related to LAI trends under conditions of stress and non-stress showed overall good performance (66%–78%). Yet, we recognize some degree of unexplained variance, and we acknowledge that model accuracy could have been improved if gridded time series of other factors that affect LAI such as plant species richness (Isbell et al., 2013) or nutrient availability (Reich et al., 2014; Reich & Hobbie, 2013) were available. For example, areas where aridity and EVs were less often significantly related are primarily located in the northern hemisphere (Figure 2), which is also where there are relatively fewer vegetation types (Dixon et al., 2014) and higher soil nitrogen availability (Nendel et al., 2019) than in the southern hemisphere. However, our analysis does succeed in providing insights into important environmental factors explaining trends in LAI under conditions of stress and non-stress, which we discuss in more detail below.

Plants in drylands are adapted to limitations in resource and water availability and are characterized by slow growth rates and a well-developed capacity for reserve storage (Berdugo, Vidiella, et al., 2022; Chapin, 1991; Lian et al., 2021). We found evidence for this adaptation as we did not observe high SHAP values and thus high impact of aridity (both average and trends) on LAI trends in general (Figure 5). Furthermore, the parabolic nature of the relationship between LAI trends and meteorological drought (Figure 6k) suggests that there is an inherent resistance to water deficit in stressed ecosystems. This is supported by experimental research that shows mild stress, for example, low availability of a resource, already activates the plants' stress-response system, which is a decline in leaf growth, photosynthesis, and in the rate of acquisition of all resources (Chapin, 1991; Zhang et al., 2020)—to prepare them for the possibility of more severe stress in the future. It is also evident that an increase in both RZ water storage capacity and OC density (soil properties) (Figure 6a,b) beyond a certain threshold does not further increase positive trends in LAI, suggesting that biophysical parameters (e.g., canopy height) and external effects (e.g., cloud cover, average temperature, burned area, and population density) may play a larger role in determining LAI trends, especially under stress. Burn frequency, for instance, showed high SHAP values for both stressed and non-stressed systems and an increase in burn frequency is linked to a corresponding increase in LAI trends. As drylands comprise many of the world's fire-prone ecosystems, such as xeric savannas (Abdi et al., 2022), postfire vegetation recovery in these ecosystems is rapid due to the increased availability of nutrients (Kellman et al., 1985) and the reproduction strategies of different species (Zirondi et al., 2021). It is thus reasonable to assume that this rapid recovery manifests itself as an increase in LAI, which significantly alleviates the negative effect of fire on the carbon balance (Zhou et al., 2022). Canopy height was further found to have a positive (linear) relationship to LAI trends (Figure 6c), which is stronger in stressed than in non-stressed systems. This is best understood in combination with land cover (Figures 4c and 6o): stressed systems are primarily found in grasslands and shrubland (relatively low canopy height), which suggests that under stress conditions there could be a transformation from grass to shrub and a corresponding increase in LAI. This is further supported by the fact that many of the areas with increasing LAI under stress match areas of abrupt positive changes in vegetation (NDVI) found in Berdugo, Gaitán, et al. (2022), which suggests the role of land use changes such as land abandonment (Bernardino et al., 2020). Non-stressed systems where LAI is increasing were mostly found in croplands (Figure 4c), so a moderate impact of canopy height also makes sense.