Cropland expansion drives vegetation greenness decline in Southeast Asia

Land-use and land-cover change (LUCC) is a key factor in determining regional vegetation greenness, impacting terrestrial carbon, water, and energy budgets. As a global LUCC hot spot, Southeast Asia has experienced intensive cropland and plantation expansion over the past 50 years; however, the impacts of these changes on the regional greenness have not been elucidated. Here, we harmonized multiple land-cover datasets, and used the satellite-derived leaf area index (LAI) in combination with a machine learning approach to quantify the impacts of LUCC on vegetation greenness in insular Southeast Asia (i.e. Peninsular Malaysia, Sumatra, and Borneo). We found that the regional LAI shows almost no trend (0.04 × 10⁻² m² m⁻² yr⁻¹) from 2000 to 2016: the net effect of an increased LAI (+5.71 × 10⁻² m² m⁻² yr⁻¹) due to CO₂ fertilization is offset by a decreased LAI mainly due to cropland expansion (−4.46 × 10⁻² m² m⁻² yr⁻¹). The impact of croplands on greenness in Southeast Asia contrasts with that in India and China. Meanwhile, oil palm expansion and climate change induced only small decreases in the LAI in Southeast Asia (−0.41 × 10⁻² and −0.38 × 10⁻² m² m⁻² yr⁻¹, respectively). Our research unveils how the LAI changes with different LUCC processes in Southeast Asia and offers a quantitative framework to assess vegetation greenness under different land-use scenarios.

1 Introduction

Terrestrial vegetation plays a pivotal role in regulating ecosystem services, conserving biodiversity, and mitigating climate change impacts. Over recent decades, long-term satellite records of the leaf area index (LAI) have disclosed a notable increase in vegetation greenness on Earth (Chen et al., 2019; Zhu et al., 2016). While elevated atmospheric CO₂ concentrations and climate change are regarded as the driving factors for vegetation greening at the global scale (Zhu et al., 2016), land-use and land-cover change (LUCC) can also substantially impact greenness at the regional scale. Previous studies have found that cropland intensification and afforestation are the primary drivers of greening in India and China (Chen et al., 2019; Kuttippurath and Kashyap, 2023). Meanwhile, other studies have reported that deforestation for croplands or pastures serves as a key driver of decreasing greenness in the Amazon (Chen et al., 2019; Querino et al., 2016). Thus, LUCC can either increase or decrease vegetation greenness, depending on the prior land use (e.g. forests and pastures), the subsequent land use (e.g. croplands, pastures, and plantations), and the intensity of these land uses (Wang and Friedl, 2019).

Southeast Asia harbours diverse biodiversity and ecosystems. However, the trends and drivers of regional greenness remain largely underexplored. Previous studies have revealed that CO₂ fertilization is a primary driver of the greening trend in Southeast Asia within a global context (Chen et al., 2022; Chen and Guestrin, 2016). Thus, climate change, especially temperature rise, could reduce vegetation growth in the tropics (Piao et al., 2020a) or drive “green-up” in maritime Southeast Asia during El Niño (Satriawan et al., 2024). More importantly, land-use change, particularly deforestation, is the predominant factor causing the decline in greenness in tropical countries like Indonesia (Chen et al., 2019; Piao et al., 2020a). However, these prior studies have primarily focused on a global scale and, thus, likely oversimplified regional dynamics, such as complex land-use processes.

Since the 1950s, Southeast Asia has been a global LUCC hot spot (Houghton and Nassikas, 2017), with maritime countries such as Indonesia and Malaysia having the greatest deforestation rates in the world (Harris et al., 2012). An increase in food crops and export-oriented crop production has driven a significant transformation of tropical forests into commodity plantations like oil palms or croplands for food (Fagan et al., 2022; Zeng et al., 2018). Indonesia and Malaysia are the largest producers of palm fruit, with 250×10⁶ and 97×10⁶ t of palm fruit produced in 2020 (FAOSTAT, 2022), respectively. Meanwhile, Indonesia and Malaysia have experienced tree cover losses of approximately 29.4×10⁶ and 8.92×10⁶ ha, respectively, over the past 2 decades, equivalent to 18 % and 30 % of their respective tree cover in 2000 (Global Forest Watch, https://www.globalforestwatch.org/, last access: 28 November 2024).

Despite the substantial LUCC in the past years, we lack a clear understanding of the impacts of LUCC on vegetation greenness in Southeast Asia. This is partly due to the complexity of the recent land-use history of the region. For example, in Indonesia, the dominant LUCC types in the 2000s were the conversion of forests to oil palm plantations and cropland in lowland regions; however, in the 2010s, the conversion of forests to oil palm slowed down, shifting more towards highland croplands and the rotation of plantations (Descals et al., 2021; Xu et al., 2020; Zeng et al., 2018). These LUCC processes can differentially affect vegetation greenness and biogeochemical cycles (Ito and Hajima, 2020), with further feedback on vegetation greenness (Wang and Friedl, 2019). However, current studies on the assessment of LUCC impacts do not often distinguish these individual land-use processes; instead, they categorize the loss of forests under one “deforestation” category (Sitch et al., 2015) or regard plantation as similar to natural forests (Hansen et al., 2013).

In this study, we aim to assess the impact of LUCC on vegetation greenness in Southeast Asia and quantify the contributions of the different LUCC processes to the changes in greenness. We collected and harmonized various types of land-cover datasets (Chini et al., 2021; Hansen et al., 2013; Sulla-Menashe et al., 2019; Xu et al., 2020) to build a detailed land-use history for Southeast Asia from 2000 to 2016. We further used a machine learning approach to quantify the impacts of land uses on the LAI, along with the impacts from climate, CO₂ concentrations, stand age, etc. Our machine learning approach, combined with hypothetical scenarios that simulate vegetation greenness without LUCC processes, enables us to isolate the impacts of different land-use changes on the LAI by estimating the difference between scenario-based LAI values.

2 Materials and methods

2.1 Study areas

Depending on the availability of various land-use data, we focused our study on a region including Peninsular Malaysia, Sumatra, and Borneo (Fig. 1) in Southeast Asia, which experience rapid LUCC (Geist and F, 2001; Mao et al., 2023). Since the 1980s, insular Southeast Asia has lost at least 1.0 % of its forests annually (Felbab-Brown, 2013; Miettinen et al., 2011), primarily due to cropland and plantation expansion (Wang et al., 2023; Wicke et al., 2011; Xu et al., 2020). Particularly notable is the expansion of oil palm plantations in Indonesia and Malaysia (Chen et al., 2024; Euler et al., 2016). From the 1990s to the 2010s, the extent of oil palm plantations increased from 1.3×10⁶ to 7.7×10⁶ ha in Indonesia and from 2.1×10⁶ to 5.2×10⁶ ha in Malaysia (Xu et al., 2022). A nation-wide study reported that around 55 %–59 % of oil palm plantation expansion in Malaysia and at least 56 % in Indonesia occurred on lands previously covered by forests during the period from 1990 to 2005 (Koh and Wilcove, 2008; Vijay et al., 2016). In parallel, cropland expansion also drives deforestation in our study area, with approximately 15 % of forest loss in Indonesia attributed to this cause (Austin et al., 2019). Rubber, timber, and other plantations have also resulted in deforestation. A recent study revealed that, between 2001 and 2016, approximately 20 % of rubber plantations in Indonesia and 33 % in Malaysia were established on land previously covered by forests, resulting in a loss of about 1×10⁶ and 0.32×10⁶ ha of forest in these countries, respectively (Wang et al., 2023).

Figure 1 Workflow of the study: steps 1 to 4 outline the processes for harmonizing multiple land-cover datasets; steps 5 to 6 show the establishment and interpretation of the LAI prediction and the scenario simulation process.

2.2 Identification of the greening trend

The LAI indicates the total amount of one-sided leaf area per unit of ground surface area (Chen and Black, 1992; Watson, 1947) and often serves as a measure of vegetation greenness (Zhu et al., 2016). In this study, we use the GLOBMAP LAI dataset, which provides a global record of vegetation cover with a 500 m resolution and is one of the main LAI datasets used for global greenness studies (Piao et al., 2020b; Zhu et al., 2016). We used the GLOBMAP LAI (version 3) dataset post-2001, which was generated based on the Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance (Liu et al., 2012), with an advanced algorithm to consider canopy clumping, making it particularly suitable for dense canopies in the tropics (Fang et al., 2019). To assess the trend in the LAI for individual pixels, we utilized the Mann–Kendall test, a non-parametric statistical method that can effectively identify consistent upward or downward trends over time (Mann, 1945). This test provides the pixel-by-pixel magnitude (β) and statistical significance (p values) of the greening trends.

2.3 Mapping different land-use transitions

In our study area, we considered natural forests, oil palms, and croplands as the major land-use types, as they together accounted for over 90 % of the land cover. To delineate the annual changes in all of these land-use types from 2001 to 2016, we harmonized various land-use datasets, which were provided at different spatial resolutions (Table S1), into a unified dataset on a 500 m grid. This harmonization was essential to make land-use changes across various datasets compatible with each other and to match the spatial resolution of the land-use data to the LAI dataset. As shown in Figs. 1 and S1, the workflow for harmonizing multiple land-cover datasets involved the following steps:

1. We first determined the annual percentages of forested (A %) and non-forested areas (B %) within each 500 m grid cell by aggregating the mean of the annual 30 m resolution Global Forest Change (GFC; v1.11) maps (Hansen et al., 2013), using the “reduceResolution()” function in Google Earth Engine (https://developers.google.com/earth-engine/guides/resample, last access: 28 November 2024).

2. Within the forested fraction of each grid cell, we estimated the proportion of oil palm (OP) plantations (A1 %) using an openly available dataset that covers the OP distribution across Malaysia and Indonesia from 2001 to 2016 (Xu et al., 2020). To estimate the proportion of OP, we calculated the frequency of oil palm pixels in each 500 m × 500 m window.

3. After accounting for the OP area, the remaining forested area in each grid was further categorized into the evergreen broadleaf forest (EBF) (A2 %) and other forest types (i.e. deciduous broadleaf forest, coniferous forest, mixed forest, etc.), based on the ratio of EBF to the total forested area provided by the MODIS Land Cover Type Product (MCD12Q1 v6.1) (Sulla-Menashe and Friedl, 2018).

4. Within the non-forested fraction of each grid cell, we used the latest version of the Land-Use Harmonization (LUH2) dataset (Hurtt et al., 2020) to estimate the percentage of cropland (CRO) (B1 %) and other non-forest land uses (i.e. pasture, grass, etc.). In this analysis, we assumed that the fraction of each land-use type in the LUH2 dataset on a 0.25° grid is applicable to the 500 m grid cells within each 0.25° grid cell. To avoid biases when downscaling land-cover data from a coarse to fine resolution, we also tested our results by upscaling all datasets to 0.25° (Figs. S2, S3).

In the end, we obtained detailed information for the EBF, OP, CRO, and Other (including other forests and non-forest vegetated areas) land-use types, at a 500 m spatial resolution. We grouped other forests and other non-forest vegetated areas together, as they represented a minor proportion (less than 5 %) of the land surface (Table S2) and exhibited minimal changes during the study period.

2.4 Extreme gradient boosting model

Tree-based machine learning models, such as extreme gradient boosting (XGBoost), have been widely used in predicting and analysing ecosystem dynamics (Green et al., 2022; Yuan et al., 2019). Compared with neural networks, which often function like “black boxes”, tree-based models offer greater interpretability and are particularly effective on tabular-style datasets (Lundberg et al., 2020). XGBoost is an ensemble learning algorithm that iteratively constructs multiple decision trees and has proven to be effective for both classification and regression tasks (Chen and Guestrin, 2016; Yan et al., 2020). This algorithm employs shrinkage techniques and performs multithreaded calculations to minimize overfitting (Meng et al., 2021).

In our study, we applied the XGBoost algorithm (Chen and Guestrin, 2016) to model the spatiotemporal variations in the mean annual LAI using climatic and LUCC factors as inputs. These factors include the fractions of EBF, OP, CRO, and Other land uses; EBF-stand and OP-stand ages (Besnard et al., 2021; Danylo et al., 2021); CO₂ concentrations (https://gml.noaa.gov/ccgg/trends/, last access: 28 November 2024); and climatic variables (Table S3). The climatic variables in our study include the mean annual temperature (MAT), mean annual precipitation (MAP), wind speed (WS), shortwave radiation (RAD), and relative humidity (RH). These gridded climatic variables were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis product v5 – Land (https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land?tab=overview, last access: 28 November 2024). We aggregated the original hourly data to the annual time step using an annual average.

To fine-tune the parameters of our XGBoost model for LAI prediction, we randomly split the data into training and testing sets using an 80 % : 20 % ratio. We then utilized the GridSearchCV method to test different parameter combinations (i.e. varying numbers of trees from 150 to 400, tree depths from 5 to 15, and learning rates between 0.01 and 0.1) and determined the best parameter combinations through cross-validation (Pedregosa et al., 2011).

2.5 Shapley additive explanations

We utilized the TreeExplainer-based SHapley Additive exPlanations (SHAP) framework to interpret the individual and interactive contributions of LUCC and other factors (including climate variables, CO₂ concentration, and stand ages) to the LAI variations in our XGBoost model. The SHAP methodology, which is based on the concept of Shapley values in cooperative game theory, offers an insightful interpretation of factor importance (Lundberg and Lee, 2017). It provides detailed, instance-specific explanations, termed SHAP values, to quantify the impact of each factor on the model predictions (Lundberg et al., 2020; Yang et al., 2021).

In the SHAP framework, the value for a given factor i in a particular sample x is computed as the average marginal contribution of that factor across all possible combinations. This is mathematically represented as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\varphi }}_i\left(x\right)=\sum _{S\subseteq N\mathit{\left\{}i\mathit{\right\}}}{\displaystyle \frac{\left|S\right|\mathrm{!}\left(\left|N\right|-\left|S\right|-\mathrm{1}\right)\mathrm{!}}{\left|N\right|\mathrm{!}}}\left[f\left(S\cup \left\{i\right\}\right)-f\left(S\right)\right],\end{array}$$

where N is the set of all factors, S is a subset of factors, and f is the prediction model. This formula quantifies the contribution of factor i by comparing the prediction with and without the factor, averaging over all possible subsets of factors.

The SHAP value indicates the magnitude and direction of the impact of a factor on prediction in the specific sample. To be specific, the magnitude (absolute value) of a SHAP value indicates the importance of a factor. Larger absolute SHAP values mean that the factor has a greater impact on the model's output. The sign of a SHAP value (positive or negative) shows the direction of the impact. A positive SHAP value indicates that the factor positively affects the model's output (e.g. increases LAI), whereas a negative SHAP value suggests a negative impact (e.g. decreases LAI). By aggregating the mean SHAP values of all samples, we can also derive global factor importance, which offers a holistic view of variables affecting annual LAI variations.

Furthermore, SHAP values aid in the interpretation of the interactive effects of two or more factors in machine learning models (Lundberg and Lee, 2017). The interactive effect is defined as the change in prediction when the joint contribution of two or more factors is considered, by subtracting the individual contributions made by each factor. The interactive effects of ith and jth factors are expressed as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\varphi }}_{i}j\left(x\right)& =\sum _{S\subseteq N\mathit{\left\{}i\mathit{\right\}}}{\displaystyle \frac{\left|S\right|\mathrm{!}\left(\left|N\right|-\left|S\right|-\mathrm{2}\right)\mathrm{!}}{\left|N\right|\mathrm{!}}}[f\left(S\cup \left\{i,j\right\}\right)\\ & -f\left(S\cup \left\{i\right\}\right)-f\left(S\cup \left\{j\right\}\right)+f\left(S\right)].\end{array}\end{array}$$

We utilized the XGBoost and scikit-learn packages in Python 3.11.0 to develop and train the XGBoost model (Chen and Guestrin, 2016; Pedregosa et al., 2011). Then, we employed the TreeExplainer function from the shap package (Lundberg and Lee, 2017) to interpret the impact of factors on LAI predictions.

2.6 Simulation scenarios

To quantify and compare the impacts of specific LUCC processes, climate change, and elevated CO₂ concentrations on vegetation greenness changes, we adopted the scenario simulation framework from several factorial attribution analyses (Sitch et al., 2015). Specifically, we first estimated the LAI trend under five hypothetical scenarios (S1–S5) using the established XGBoost model. The equations utilized are as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathrm{LAI}}_{i,t,\mathrm{Sn}}& =\mathrm{XGBoost}({{\mathrm{CO}}_{\mathrm{2}}}_{i,t,\mathrm{Sn}},{\mathrm{CLI}}_{i,t,\mathrm{Sn}},\\ & f\mathit{\_}{\mathrm{EBF}}_{i,t,\mathrm{Sn}},\mathrm{f}\mathrm{\_}{\mathrm{CRO}}_{i,t,\mathrm{Sn}},\\ & \mathrm{f}\mathrm{\_}{\mathrm{OP}}_{i,t,\mathrm{Sn}},\mathrm{f}\mathrm{\_}{\mathrm{Other}}_{i,t,\mathrm{Sn}}),\end{array}\text{(4)}& {\displaystyle }\mathit{\beta }{\mathrm{LAI}}_{i,\mathrm{Sn}}=\mathrm{slope}\left({\mathrm{LAI}}_{i,\mathrm{Sn}}\right).\end{array}$$

Here, ${\mathrm{LAI}}_{i,t,\mathrm{Sn}}$ represents the simulated LAI for the ith grid at year of t under scenario Sn; βLAI_i,Sn indicates the LAI trend over the study period; and XGBoost stands for the established model for LAI prediction using the CO₂ concentration (CO₂), climate variable (CLI), and land-cover types, such as the fraction of evergreen broadleaf forest (f_EBF), cropland (f_CRO), oil palm (f_OP), other land uses (f_Other) (see Sect. 2.4).

For different scenarios, we adjusted the input variables according to specific assumptions to progressively incorporate different factors. For S1, we assumed that only the CO₂ concentration varied from 2001 to 2016, while climate and land-use variables (i.e. CLI, f_EBF, f_CRO, f_OP, and f_Other) remained constant at their 2001 values. For S2, CO₂ and climate change over time, with land uses remaining unchanged at their 2001 values. S3 to S5 sequentially considered different land-use processes: S3 involved changes from EBF to CRO using time-varying CO₂, climate, and CRO fractions, while keeping OP and other land-use types constant post-2001; S4 included conversions from EBF to both CRO and OP using time-varying CO₂, climate, CRO and OP fractions, while other land uses remained unchanged post-2001; and S5 encompassed all LUCC changes, with all variables including CO₂, climate, and all LUCC types varying over time.

We then quantified the impacts of each factor on vegetation greening based on differences in LAI trends between scenarios as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{Driver}}_n=\mathit{\delta }\text{LAI trend}=\mathit{\beta }{\mathrm{LAI}}_{i,\mathrm{Sn}}-\mathit{\beta }{\mathrm{LAI}}_{i,S(n-\mathrm{1})}.\end{array}$$

Here, Driver_n measures the impact of the nth driver on LAI trends. The drivers include CO₂, climate change, CRO expansion, OP expansion, and Other LUCC. Notably, Driver₁ quantifies the impact of CO₂, equal to βLAI_i,S1. In addition, our estimation of CRO or OP expansion assumed that the increased areas of CRO or OP since 2001 came from EBF, given that CRO and OP expansion mostly resulted from deforestation in Indonesia and Malaysia (Numata et al., 2022; Wagner et al., 2022).

3 Results

3.1 Land-use changes and greening trends

From 2001 to 2016, the extent of forest in our study area decreased annually by 1.29 %, leading to an EBF reduction from 73.41 % to 53.09 % of the study area. Notably, approximately 25.54 % of the region experienced rapid deforestation, with a forest loss rate exceeding 2 % of the land surface per year. This deforestation rate was especially pronounced in the eastern parts of Sumatra and along the western and southern edges of Borneo (Fig. 2b).

In the deforested areas, we observed the widespread expansion of CRO and OP plantations. The CRO area increased at a rate of 0.63 % yr⁻¹, resulting in an increase in CRO area from 14.45 % to 24.56 % of the region (i.e. 20.76 × 10⁴ to 35.28 × 10⁴ km²) between 2000 and 2016. Meanwhile, the expansion of OP proceeded at a pace of 0.48 % yr⁻¹, resulting in a nearly tripled extent of oil plantations over the past decade (i.e. from 3.91 % in 2001 to 12.05 % in 2016; Fig. 2a). In central Sumatra, the southern edge of Borneo, and the southern part of Peninsula Malaysia, OP showed the largest increase, partly at the expense of decreases in CRO (Fig. 2c).

Figure 2 Land-use composition and its changes in the study area from 2001 to 2016: (a) the changes in the fractions of main land uses, including evergreen broadleaf forest (EBF), oil palm (OP), cropland (CRO), and Other land uses over the study period; (b–d) spatial distribution of the trend in EBF, CRO, and OP as a fraction of land area.

The regional average LAI exhibited a non-significant upward trend over the study period, with a slope of 0.04 × 10⁻² m² m⁻² yr⁻¹ (Fig. 3a). Using another LAI dataset (MODIS LAI), we also captured an increasing trend (Fig. S4). According to Galán-Acedo et al. (2021), areas with over 70 % forest loss are categorized as having high to severe deforestation, whereas areas with less than 70 % forest loss are classified as undergoing low to intermediate deforestation. We found that, for our study area, the non-significant upward trend in the LAI was due to the net effect of a rapid LAI increase ($\mathit{\beta }=+\mathrm{1.07}$ × 10⁻² m² m⁻² yr⁻¹, p<0.05) in areas with low to intermediate deforestation and a significant LAI decline ($\mathit{\beta }=-\mathrm{3.08}$ × 10⁻² m² m⁻² yr⁻¹, p<0.001) in areas with high to severe deforestation (Fig. 3a). Across the study area, 58.50 % of the region showed a significant decreasing LAI trend; these areas were mostly regions with pronounced forest loss (Figs. 2b, 3b).

Figure 3 Trends in the LAI in the study area from 2001 to 2016. Panel (a) shows the regional average LAI trend, the LAI trends in regions with high to severe deforestation, and the LAI trends in regions with low to intermediate deforestation. The classification of the deforestation level is taken from Galán-Acedo et al. (2021): areas experiencing more than 70 % forest loss are classified as high or severe deforestation, whereas those with less than 70 % loss are classified as low or intermediate deforestation. Panel (b) presents the spatial pattern in the LAI trend, with the inset histogram showing the frequency (%) of the pixel-wise LAI trend in the study area.

3.2 Drivers of changes in the LAI

To understand the variations in vegetation greenness in Southeast Asia, we established an XGBoost model to quantify the relationship between vegetation greenness and land uses, climate variables, CO₂ concentrations, and stand ages. The XGBoost model showed high explanatory power (i.e. 98 % accuracy for calibration and 93 % accuracy for validation), underscoring the model's reliability for analysing the determinants of LAI variability (Fig. 4). Using different splitting ratios for calibration and validation dataset, we obtained similar model performance (Fig. S5).

Figure 4 Comparison of observed and predicted LAI values with the XGBoost model for the (a) calibration and (b) validation datasets.

Based on the built XGBoost model, we evaluated the contributions of each factor to the LAI in the SHAP framework (Fig. 5; see Sect. 2). Our results reveal that the fraction of EBF (f_EBF) in each grid holds the greatest mean absolute SHAP value (1.28), indicating that f_EBF had the largest impact on the LAI (Fig. 5b) and that this impact was positive (Figs. 5a, S6a).

The average impact of the fraction of OP (f_OP) on the LAI ranked the second largest (0.21) and had a similar positive impact on the LAI to f_EBF (Figs. 5a, b; S6b). In contrast, the impact of CRO and other land uses on the LAI was found to be negative (Figs. 5a; S6c, d). A higher fraction of CRO (f_CRO) led to a larger negative SHAP value. We further explored the impacts of the interactions of land-use types (i.e. f_CRO, f_OP, and f_EBF) on the LAI. The results demonstrated that, in areas with low f_EBF, an increase in f_OP enhanced the LAI, whereas f_CRO induced LAI decreases. Meanwhile, in areas with high f_EBF, the impacts of both f_OP and f_CRO on the LAI were markedly reduced, suggesting a dominant role of f_EBF in greenness in the region (Fig. 5c, d).

Apart from the impacts of land uses on the LAI, we also found that an elevated CO₂ concentration substantially increased the LAI, with a mean SHAP value of 0.18. In contrast to an elevated CO₂ concentration, the MAT was negatively related to the LAI, with a smaller mean SHAP value of 0.08 (Figs. 5b, S6e). Other climate variables have limited impacts on the LAI. It is noteworthy that the stand ages of both EBF (Age_EBF) and OP (Age_OP) positively impact the LAI. Specifically, Age_EBF has a greater impact than that of Age_OP, as indicated by a higher absolute mean SHAP value for Age_EBF (0.11) compared with Age_OP (0.08). While we found that Age_EBF continuously contributed to LAI increases (Fig. S6k), Age_OP increased the LAI at a younger age (less than 12 years) and then decreased the LAI afterward (Fig. S6l).

Figure 5 The impact of factors on LAI. In panel (a), beeswarm plots show the SHAP values of each LAI factor for each sample. The SHAP value indicates the magnitude and direction of the impact on LAI (see Sect. 2). Each dot represents an individual sample, with the colour indicating the relative values of the specific factor. Panel (b) presents a bar plot of the mean absolute SHAP values of each LAI factor. Panel (c) shows the interaction of f_OP and f_EBF on the LAI, whereas panel (d) presents the interaction of f_CRO and f_EBF on the LAI. The definitions of the abbreviations used for each factor are available in Table S3.

3.3 Impacts of LUCC on greening

Through scenario-based prediction, we quantified the impacts of LUCC, elevated CO₂ concentration, and climate change on the greening trend (Figs. 6f, S8). Compared with the observed greening trend for the study area (i.e. 0.04 × 10⁻² m² m⁻² yr⁻¹), we found that the greening trend increased to 5.71 × 10⁻² m² m⁻² yr⁻¹ under scenario S1, which simulated the effect of elevated CO₂ alone, with the climate and LUCC values remaining constant. This result suggests that elevated CO₂ had a clear positive impact on the LAI in Southeast Asia. Climate change showed a small negative impact (i.e. mostly due to rising temperature) on the LAI. We found that both CRO expansion and OP expansion decreased the LAI trend, with the trend dramatically dropping by −4.46 × 10⁻² m² m⁻² yr⁻¹ under the impact of CRO expansion and by −0.41 × 10⁻² m² m⁻² yr⁻¹ under the impact of OP expansion. These results highlight that CRO expansion was the primary reason for the decrease in vegetation greenness, counteracting the greening trend caused by elevated CO₂ in our study area. In contrast, OP expansion only contributed to a small decline in greenness.

We further examined the spatial variations in the impacts of each factor on greening by quantifying the differences in greening trends under different scenarios at the pixel level. Across the study area, LUCC imposed a negative impact on LAI trends (Fig. 6c, d, e). Consistent with regional average values, we found that the changes from EBF to CRO had a more pronounced negative impact on the greening trend than the conversion of EBF to OP (Fig. 6c, d). In some regions, such as the southern edges of Sumatra and Borneo, OP enhanced regional greening (Fig. 6d). Meanwhile, elevated CO₂ concentrations consistently had a positive impact on greening across the region (Fig. 6a), and climate change showed an overall negative but highly heterogeneous impact on LAI trends (Fig. 6b). From a spatial perspective, we found that elevated CO₂ dominated the increase in the LAI in most areas, accounting for 62.10 % of the study area, while CRO expansion was the primary driver of LAI decreases in other regions (26.33 %), especially coastal areas (Fig. S7).

Figure 6 The spatial distribution of the pixel-wise impacts of each factor on the greening trends (δLAI trend). Positive values mean that the factors considered increase the LAI trend, whereas negative values mean the inverse. We show the spatial patterns of contribution from (a) elevated CO₂ concentrations, (b) climate change, (c) expansion of cropland (CRO), (d) expansion of oil palm (OP), and (e) Other land-use changes. Panel (f) shows the average impact of each factor, with the error bars indicating 1 standard deviation. The symbol “${}^{***}$” denotes a statistically significant difference in LAI trends at the p<0.001 level.

4 Discussion

In this study, we analysed the typical LUCC processes in Peninsular Malaysia, Sumatra, and Borneo as well as their impacts on greenness over the past 2 decades. We found a significant decline in EBF coverage, from 73.41 % to 53.09 %, predominantly due to CRO and OP plantation expansion. Meanwhile, we did not find a significant trend in the LAI in our study area, as the increases in regional greenness due to elevated CO₂ were offset by the decreases in regional greenness caused by CRO and OP expansion. Notably, the negative effect of CRO expansion (−4.46 × 10⁻² m² m⁻² yr⁻¹) was more pronounced than that of OP expansion (−0.41 × 10⁻² m² m⁻² yr⁻¹), indicating a dominant role of CRO expansion on greenness decline in our study area.

4.1 Strong negative impact of cropland expansion on greenness in Southeast Asia

Our results demonstrate that CRO expansion in Southeast Asia contributed negatively to vegetation greenness. This is the opposite of existing reports on the greening trend in China and India, where CRO expansion is suggested to be the main reason for the net increase in the LAI (Chen et al., 2019; Kuttippurath and Kashyap, 2023). This disparity may be partly attributed to the original land-use types before CRO expansion. In Southeast Asia, CRO expansion mainly occurs at the expense of natural forests (Potapov et al., 2022; Zeng et al., 2018), and crops often have less-dense canopies than natural forests (Asner et al., 2005; Foley et al., 2005; Pocock et al., 2010). In contrast, the increase in the LAI in China and India primarily resulted from the intensification of croplands, rather than their expansion; moreover, where expansion did occur, it predominantly took place on lands that were previously bare or sparsely vegetated (Chen et al., 2019).

We also note that agricultural practices in Indonesia and Malaysia are generally less intensive compared with those in India and China, partly due to the less-advanced agricultural technologies deployed in the region (X. Liu et al., 2021). For example, in China and India, intensive agricultural practices, including precision fertilization and advanced irrigation systems, such as drip irrigation in China and spray irrigation in India, are widely adopted to enhance crop growth (Cui et al., 2022; Wang et al., 2013). In addition, the development of specialized crop varieties, such as hybrid rice in China and climate-smart drought-tolerant rice varieties in India (Panda et al., 2021; Zhang et al., 2022), also facilitated plant growth and, consequently, regional greening (Zhao et al., 2021; J. Zhao et al., 2023). In contrast, Indonesia and Malaysia predominantly depend on rainfed irrigation and traditional farming methods, such as the “subak” terraced rice fields in Indonesia, leading to lower crop intensity. Less-intensive practices are less likely to create a dense canopy and high biomass in croplands (Takeshima, 2019).

4.2 Minor contribution of oil palm to regional greenness

Our study observed a small negative impact of OP expansion on greenness. This aligns with studies reporting that OP has LAI values comparable to, or only slightly lower than, native forests (Propastin, 2009; Rusli and Majid, 2014; Vernimmen et al., 2007). Therefore, OP expansion, although a major driver of land-use change in Southeast Asia, was not a main driver of greenness decline in the study region.

Meanwhile, we also found that the OP stand age influences the LAI, as the LAI of OP increases with stand age when OP is young but then decreases after a threshold stand age is reached. This observation agrees with previous studies that reported a plateau and decrease in the oil palm LAI after the age of 10 (Xu et al., 2021) and a plateau and decrease in the oil palm yield after the age of 8–9 (Park et al., 2023). The stand age of most OP plantations in Southeast Asia is approaching maturity (over half were planted before 2009) (Danylo et al., 2021). In response, industrial and smallholder plantations have undergone or are starting to undergo the process of replanting (Danylo et al., 2021; Numata et al., 2022), although replanting in smallholder plantations is often delayed as farmers face more financial constraints (J. Zhao et al., 2023). Hence, we expect to see a more complex influence of OP on the LAI trend depending on the management practices of different types of plantations.

In comparison with CRO, our findings indicate that the impact of OP expansion on vegetation greenness decline is relatively minor (Fig. 6). This can be attributed to the generally higher biomass of OP compared with typical crops, like rice and maize, in Indonesia and Malaysia. In addition, we assumed that the CRO or OP expansion since 2001 came from EBF in our study. While CRO and OP expansion indeed mostly resulted from deforestation in Indonesia and Malaysia (Numata et al., 2022; Wagner et al., 2022), the expansion of CRO and OP can also come from other land uses (i.e. grasslands or pastures). We suspect that transitions from these land uses to CRO and OP might result in a smaller negative impact on vegetation greenness, considering that grasslands and shrubs have a lower LAI than EBF.

4.3 Other LUCC impacts on greenness in Southeast Asia

While our study examined the two most prominent processes of LUCC in Southeast Asia (EBF to OP and EBF to CRO), there are other types of land-use change that we analysed together under the category of “Other” land-use types. The changes in these land-use types are also relevant to deforestation, and they include other plantations, such as rubber (Wang et al., 2023); agroforestry plantations, such as cocoa and coffee (Panda et al., 2021); and grasslands or pastures (Austin et al., 2019). In total, the Other category accounted for 4.70 % of the study area in 2016, much less than the three main land-use types that we studied (51.06 % for EBF, 25.01 % for CRO, and 12.09 % for OP in 2016), and has experienced minor changes since 2001.

Meanwhile, we found that the overall impact of these other LUCC types on greenness was likely small (Fig. 6f). The small impact is contingent on their smaller extent compared with EBF, CRO, and OP (see above). It may also result from the offset of positive and negative impacts from individual other LUCC processes on greenness. For example, rubber plantations exhibit a higher LAI than natural forests (Wang et al., 2022), whereas agroforestry and other plantations generally have a lower LAI than natural forests, therefore leading to different trends in the LAI after deforestation.

4.4 Impact of climate change and CO₂ concentrations on regional greenness

CO₂ fertilization effects appear to be the primary drivers of greening trends observed in global studies (Chen et al., 2022; Ewert, 2004; Zhu et al., 2016). Our research also confirmed the substantial contribution of rising CO₂ levels to the greening of vegetation in Southeast Asia. The impact of temperature on greenness in Southeast Asia was negative, in contrast to the positive effects noted in cold climate zones, such as the Qinghai–Tibet Plateau (Yang et al., 2024; R. Zhao et al., 2023; Zhong et al., 2019) and Arctic areas (Forbes et al., 2010). We suspect that the negative effect of temperature implies that temperature may have exceeded the optimal point for plant growth in parts of Southeast Asia. This aligns with several studies suggesting that ecosystem functions in the tropics are approaching a temperature tipping point (Doughty et al., 2023; Meir et al., 2015; Wu et al., 2019). Additionally, temperature rise might exacerbate the incidence of pests and diseases in tropical forests, negatively impacting plant health and productivity (Ghini et al., 2011). These factors might jointly contribute to the observed decline in greenness with increasing temperatures in Southeast Asia.

5 Conclusions

Our study closely examined the impacts of LUCC on vegetation greenness in part of Southeast Asia. We found that there was no significant trend in vegetation greenness in the study area, which is attributed to the net effect of negative impacts of LUCC on LAI and positive influence of elevated CO₂ on the LAI. Among various LUCC processes, we found that cropland expansion was the primary reason for LAI decrease, whereas oil palm expansion had a small impact on the LAI trends. These results shed light on the interplay between greenness and land-use changes and provide valuable insight into our future studies on terrestrial carbon, water, and energy budgets in the land-use-change-intensive region of Southeast Asia.