Projecting global oil palm expansion under zero-deforestation commitments: Direct and indirect land use change impacts

Summary In the last three decades, global production of oil palm has boomed, which has partly come at the expense of tropical rainforests. Recognizing this, many companies operating in the palm oil industry have committed to eliminate deforestation from their operations, often referred to as zero-deforestation commitments (ZDCs). Here, we estimate that if ZDCs are fully adopted and enforced across all sectors and geographies, the global extent of oil palm plantations may be 11 M ha or 40% smaller in 2030 than in a business-as-usual (BAU) scenario that assumes no compliance with ZDCs. As a result of such land-sparing effects, we estimate that 96 M ha of forests are saved from conversion, of which, 17% would otherwise have been converted (directly or indirectly) due to expanding oil palm plantations. Overall, these figures suggest that ZDCs have the potential to deliver major environmental benefits if they are fully adopted and enforced.


INTRODUCTION
In the last three decades, global production of oil palm (Elaeis guineensis) fresh fruit bunches (FFBs) has boomed, increasing by nearly 600%. 1 Oil palm provides two types of vegetable oil: palm oil and palm kernel oil. These oils are used in a variety of applications including foods, soaps, detergents, cosmetics, pharmaceuticals, and biofuels. 2 Although the yields vary substantially depending on the age of the plantation, 3 the time-averaged yields of palm oil per hectare are by far the highest of all oil crops. 4 Over the last decades, palm oil has become the most consumed and traded vegetable oil in the world and demand will very likely continue to grow in the near future. [5][6][7] However, in contrast to other major agricultural commodities such as rice or wheat, increased demand for oil palm has mostly been met by expanding the total area under production rather than intensification (i.e. increasing yields on existing production areas). 8 This has resulted in an increased area of oil palm plantations of more than 350%-from 6 to 28 M ha-since 1990 (see Figure S1). 1 As oil palm only grows in the tropics, this expansion has partly come at the expense of tropical rainforests. 9 This is especially true for Southeast Asia, which is the largest palm oil-producing region of the world, accounting for 84% of production in 2018. 1 It is estimated that around 45% of all new oil palm plantations established in Southeast Asia since 1989 have replaced forests. 10 At the same time, there are large opportunities for the expansion of oil palm plantations outside tropical forests 11,12 and recent evidence suggests that recent oil palm expansion in Latin America has been largely deforestation-free. 13,14 In recent years, many companies involved in the palm oil industry have pledged to eliminate or reduce deforestation from their supply chains. These pledges are often referred to as zero-deforestation commitments (ZDCs). 15 ZDCs apply to forests within a company's supply chain, regardless of whether that company owns the land or not. As of April 2020, 83% of the palm oil refining capacity in Indonesia and Malaysia 1 -which, together, accounted for 84% of global production between 2010 and 2019-was covered by ZDCs. 16,17 While various studies have considered the past effectiveness of ZDCs, 18-20 less attention has been given to understanding how full implementation and enforcement of ZDCs could potentially play out in terms of future global land use and land use change. The difficulty in predicting the potential impacts of ZDCs arises from the necessity to construct a plausible counterfactual, i.e., a scenario of what would happen in the absence of ZDCs. Such an assessment should account for global marketmediated effects that may occur in response to the implementation of ZDCs as well as for local differences in land availability, land suitability, and baseline changes in the economy. One study 7 provided an ex ante assessment of the overall forest loss that may be avoided due to ZDCs in Indonesia but did not account for economy-wide effects and spillover effects in other countries. Such spillover effects have the potential to seriously undermine the effectiveness of ZDCs. 21 A more recent study 22 offered a comprehensive, global assessment of the amount of deforestation that may be avoided if certain economic policies are adopted by the Indonesian and Malaysian government, but they did not consider the role of ZDCs in curbing deforestation. Furthermore, both studies only assessed aggregate forest loss at the country level and did not consider spatialized deforestation outcomes within countries. As impacts depend on the spatial patterns of forest loss, more spatially explicit assessments are needed to anticipate the potential environmental benefits of the ZDCs. 23 The objective of this paper is to make a spatial assessment of how the worldwide implementation of zero-deforestation commitments across all agricultural commodities, economic sectors, and geographies could alter the spatial configuration of oil palm and other land uses up until 2030 while accounting for expected baseline changes in the economy. The analysis will facilitate an evaluation of the extent to which fully executed ZDCs may moderate conversion of forests and other (semi) natural areas over this time period. We link the results from a computable general equilibrium model with a dynamic land use change model to project how land use may evolve under two scenarios: a business-as-usual (BAU) scenario assuming no implementation of ZDCs and a ZDC scenario assuming full implementation and enforcement of ZDCs across all economic sectors. The spatial distribution of oil palm and other land uses is modeled from 2014 (the most recent reference year for which global economy-wide data are available) up until 2030. The year 2030 coincides with the target year of the 2014 New York Declaration on Forests 24 ; a declaration signed by more than 200 companies, governments, non-governmental organizations, and groups representing indigenous communities to end all global natural forest loss.

Effects on oil palm area and other types of land use
The results of the GTAP experiments show that the global oil palm area in 2030 under the ZDC scenario is nearly 3 M ha or 14% smaller relative to the area in 2014 (see Figure 1). As the global area under oil palm is expected to expand by 42% under the BAU scenario during the same period, the worldwide implementation of ZDCs leads to a 39% smaller oil palm area relative to the BAU scenario up until 2030. The reduction in oil palm area due to ZDCs is largest in Indonesia (6.3 M ha), Malaysia (2.4 M ha), and Nigeria (0.8 M ha; see Figure 2). These results are rather sensitive to how the share of oil palm in the aggregate oil crop sector of GTAP is calculated. If it is assumed that the share of oil palm in 2030 equals the observed share within each GTAP region in 2019, the global oil palm area expands by 19% between 2014 and 2030 under the BAU scenario and reduces by 28% under the ZDC scenario over the same period. Results are shown for two scenarios: a business-as-usual scenario that follows historic trends and an alternative scenario that assumes full implementation of zero-deforestation commitments (ZDCs).  Figure S13).
In addition to boosting yields, the ZDC-induced increase in land rental rates will partially translate into an increase in global commodity prices (see Figure S14) and hence a reduction in the demand for oil palm. Whereas oil palm production increases by 42% under the BAU scenario between 2014 and 2030, it only increases by 1% under the ZDC scenario. Thus, due to ZDCs, the global production of oil palm FFB relative to the BAU scenario is expected to decrease by 29% or 117 M metric tonnes.
The results also indicate that ZDCs will partially displace production of agricultural commodities to new areas and encourage intercrop substitutions. Production of oil crops is expected to increase by 7.4 M metric tonnes in non-oil palm-producing countries (which equates to 0.9% of the global production in 2014). It is likely that some of the production of oil palm will also be displaced to other oil palm-producing regions with more land availability. However, Figure 2 shows that none of these regions experience a net increase in oil palm due to ZDCs, which implies that any such effects are offset by the demand and yield effects. In addition to a major difference in oil palm area, total area devoted to other crops, forestry, and pastureland is also expected to decrease relative to the BAU scenario as a result of the large increases in land rents, although the percentage changes are smaller (ranging from À10% to À20%). Apart from oil palm, major reductions are expected in the area devoted to fruit and vegetables, coarse grains, and rice (see Figure S15).
These results remain largely robust if the alternative estimates of ZDC coverage are used (as described in the STAR Methods section), with the projected reduction in oil palm area relative to BAU conditions varying from 38% to 42%. Results for the other land use types are similar too, with the percentage reductions in area varying in the range of 9%-26%. These effects are to a lesser extent driven by oil palm as shrublands tend to be less dominant in the main oil palmproducing regions (Figure 3). Most of the avoided conversion is in Oceania (notably Australia), where iScience Article shrublands tend to be the dominant natural biome and where many of these areas are covered by ZDCs due to their high levels of above-ground biomass carbon. 25 By comparison, only a small area of grassland is saved from conversion (3 M ha), while a much larger area is converted due to ZDCs (23 M ha, of which less than 1% can be attributed to oil palm). The reason why displacement effects for grassland are larger than for other biomes is that most grasslands fall beyond the scope of ZDCs, which is why cropland expansion tends to be redirected to grasslands. When taken in aggregate, the total extent of natural land that is saved from conversion is 217% larger than the area converted due to ZDCs. This suggests that the ZDCs are likely to deliver substantial environmental benefits within the oil palm-producing world, although there is considerable heterogeneity across regions.

DISCUSSION
This is the first study to provide a modeling experiment of how the worldwide implementation of zerodeforestation commitments could alter the spatial configuration of oil palm and other land uses up until 2030, assuming full adoption and enforcement across industries and regions from 2014 onward. The results suggest that under these assumptions, ZDCs are likely to bring about significant land-sparing effects. Due to increases in yields and an overall decrease in demand for oil crops, ZDCs could hypothetically induce a decrease in global oil palm area of 11 M ha or 39% relative to a BAU scenario that assumes no compliance with ZDCs and that follows historic trends. By way of comparison, FAO statistics for the period 2014-2020 report an increase of 6.8 M ha in the global extent of oil palm plantations, of which 75% occurred in the period 2014-2017. This compares with our estimate of 8.1 M ha over the period 2014-2030 under BAU conditions, which appears to be a reasonable projection in light of the flattening trend observed from 2017 onward. These results remain largely robust when using alternative estimates of the potential spatial coverage of ZDCs. A potential explanation for this is that the spatial differences between these different estimates are relatively minor after masking out areas unavailable for oil palm production, suggesting there is little additional impact when considering more ambitious forest protection scenarios. Apart from oil palm, ZDCs are also expected to induce a decrease in the total area devoted to other crops, forestry, and pastureland in the oil palm-producing world relative to the BAU scenario. Although the total area under forestry and pastureland is projected to decrease regardless of ZDCs (as would be expected given the large contraction in recent years 26,27 ), both types of land use will still be major drivers of local forest loss under the BAU scenario. Overall, ZDCs are estimated to prevent the conversion of around 96 M ha of forests, of which 17% would otherwise have been converted due to expanding oil palm plantations.
Notwithstanding the large land-saving effects, potential environmental benefits associated with reductions in nature loss are likely partially offset by displacement effects to natural areas that fall beyond the scope of ZDCs. Not only does this include displacement effect within oil palm-producing regions but also iScience Article displacement effects to temperate regions. This suggests that temperate oil crops such as rapeseed, sunflower, and soybean are expected to partially substitute for oil palm, which is particularly likely in the case of biofuels, 28 but also for many food applications where hydrogenation of alternative edible oils is a viable option. 29 In addition, although yield increases help to free up crop areas for other uses, anticipated increases in the use of fertilizers, pesticides, and heavy machinery (likely needed to achieve these yield increases) could lead to larger local and global environmental impacts from existing agricultural areas. 30 The degree to which ZDCs are to deliver environmental benefits thus partially hinges on the impact of displacement effects and the technology through which yield increases are achieved. We also find that due to land scarcity and increasing commodity prices, ZDCs are expected to depress consumption of agricultural commodities, posing potential challenges for food accessibility. What is more, imposing ZDCs on individual supply chains without considering the differential capabilities of actors to comply with them could result in the unfair market exclusion of farmers. For example, smallholder farmers may be unable to comply with deforestation-free certification schemes such as the Roundtable of Sustainable Palm Oil by virtue of the high costs involved in certification audits, lack of awareness of ZDCs, or unclear land tenure. 31-33 Such concerns of distributive equity could be mitigated if ZDCs are accompanied by differentiated and locally targeted capacity-building measures. These could include financial support schemes, knowledge dissemination activities, and coordination efforts with other policy-making actors (private and public) to enhance the inclusivity of policies and to exploit synergies between ZDCs and equity. Previous research has shown that such synergies exist as poor communities stand to benefit the most from forest protection measures by virtue of the improved pollination services, forestry, and fishery yields they provide. 34

Limitations of the study
Our findings build on a burgeoning literature that has found that anti-deforestation policies, if implemented at large scales, could induce major environmental benefits outside targeted areas. 22,35 However, previous studies have also found that the risk of leakage effects is much higher if adoption and implementation varies widely across space. [36][37][38] Our results should thus be interpreted carefully as they are projected on the assumption that all industries and regions will be fully covered by ZDCs. It is likely that the anticipated land-sparing effects will be considerably diminished if large consumer markets fall beyond the scope of ZDCs, as is currently the case in China and India. 39 Moreover, the results are rather sensitive to how the share of oil palm in the aggregate oil crop sector of GTAP is calculated, with much less oil palm expansion projected under the BAU scenario between 2014 and 2030 (19% instead of 42%) if the share of oil palm relative to other crops is assumed to remain constant rather than growing at historic rates. However, given that recent national-level statistics show a major increase in the share of oil palm over the period 2014-2019 for most oil palm-producing regions (see Figures S7A-S7C), this appears to be a less plausible assumption. Our analysis is subject to some additional uncertainties that we cannot fully address. First and foremost, it is possible that the region-specific land supply curves underpinning our analysis imply land supply elasticities that are too high for some regions. As a result, we may have overestimated the magnitude of the overall market-mediated effects of ZDCs. The land supply elasticities may be too high because data on current land prices are only available for a limited number of regions in the world. 40 For regions with no data on land rents, it is assumed that the marginal land rents are inversely proportional to marginal yields. While this assumption may be a reasonable approximation of the land rent trajectory in many developed countries, 41 it is less likely to hold true in developing countries with poorly functioning land markets and weakly enforced land tenure rights. 42 A related point is that the analysis is based on a shift in the land supply asymptotes that are assumed to remain unaltered throughout the simulation period. In practice, it is more likely that such asymptotes are dynamic as they depend on annual weather fluctuations, changes in soil quality, the spatial dynamics of legally protected areas, land governance, and the degree of technological progress. A case in point is the Brazilian Cerrado, a savannah ecoregion that was largely considered unsuitable for agricultural production in the 1970s until advances in soil amendment technologies and crossbreeding transformed it into the world's ''soy basket''. 8 It should be noted, though, that given the relatively short time period for our simulation (2014-2030), it is unlikely that such unforeseen dynamics will have a major impact on the amount of available land until 2030. A final major source of uncertainty is that our estimates of the overall impact of ZDCs are dependent on the validity of our counterfactual projections which describe what might happen in the absence of ZDCs (i.e., the BAU scenario). Although the BAU scenario largely follows the SSP2 scenario, which is considered a benchmark for baseline projections, 43 it is likely that there will be deviations between our projections and the actual changes in oil palm area. These deviations may occur due to political economy factors that are not well predicted by a computable general equilibrium model iScience Article (e.g., biofuel mandates, election cycles, and imperfect law enforcement). In addition, we may have underor overestimated future crop yields as they are currently assumed to grow at a linear rate. Although this assumption is largely consistent with historic crop yield trajectories, 44 future crop yields remain a key source of uncertainty as it is unknown what capital investments would be required and how this would influence the adoption rate. 45 On a similar note, our results depend to a large extent on the factor mobility assumptions underpinning the GTAP model, especially with regard to land use type-specific conversions. While the disaggregation of each GTAP region into AEZs region in the GTAP region is an attempt to account for spatial differences in mobility of land across uses, there is still considerable heterogeneity within these units. Future research should thus explore the sensitivity of our results to alternative yield and mobility assumptions. Despite these uncertainties, our study provides a first estimate of how full implementation and enforcement of ZDCs globally may moderate the expansion of oil palm and its encroachment into forest areas. It provides a benchmark against which future estimates can be compared, for example assessments of ZDC implementation with incomplete coverage across sectors and geographies. Although it is unlikely that the coverage of ZDCs will be anywhere near 100% in the coming decade, our study provides strong quantitative evidence of the major reduction in deforestation and agricultural expansion they could deliver if they are adopted and enforced at large scale, thus resulting in potentially large environmental benefits. Given the goals of the New York Declaration on Forests and the Convention on Biological Diversity, this should motivate the international community to increase the uptake, but also the enforcement of ZDCs. Imposing mandatory due diligence rules for importers into major destination markets of deforestationrisk commodities, as currently discussed by the European Union, United Kingdom, and United States, 46 could help accelerate progress toward making ZDCs more effective. In addition, care should be taken to engage participation of stakeholders from tropical forest countries in zero-deforestation policymaking to ensure not only the effectiveness of such policies but also the equitable conservation of tropical forests.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Materials availability
The gridded land supply asymptotes at 1 3 1 km for all region-specific Agro-Ecological Zones in the GTAP-AEZ database under the Business-As-Usual (BAU) scenario can be downloaded from the Environmental Geography website: https://www.environmentalgeography.nl/. The gridded land supply asymptotes under the Zero-Deforestation Commitments (ZDCs) scenarios are partially based on IUCN Red List Data and therefore, cannot be redistributed without the prior written permission of IUCN: https://www. iucnredlist.org/resources/spatial-data-download.

Data and code availability
All the source data used for this study are publicly available (see Tables S1 and S3). All code used to analyze the data in this paper is available from the lead contact upon reasonable request.

Overall approach
We took a stepwise approach to project changes in oil palm area from 2014 up to 2030, as a result of fully implemented and enforced ZDCs. We first employed a multi-commodity, multiregional comparative static computable general equilibrium model to make baseline projections of economic activity up to 2030, assuming no compliance with ZDCs. This model operates at a spatial resolution of 37 world regions, stratified by 18 global Agro-Ecological Zones (AEZs), resulting in 337 region-specific AEZs (see Figures S3 and  S4), and accounts for legally protected areas using data from the World Database on Protected Areas. 48 To account for heterogeneity in the area available for cropland expansion across the different regions, we calibrated the model with new spatially explicit estimates of region-specific supply of available land, also referred to as land supply asymptotes. We then constructed an alternative scenario that assumes full implementation of and compliance with ZDCs but otherwise equal to the BAU scenario. This scenario implies Spatial predictor variables used in the CLUMondo simulations This paper Table S3 Land conversion matrix used in the CLUMondo simulations This paper Table S4 Conversion resistance elasticities for each land system used in the CLUMondo simulations This paper Table S5 List of neighborhood weight factors by land system This paper iScience Article a reduction in the supply of land available for expansion equal to the areas covered by ZDCs, and hence, simulates a potential change in both the supply and demand of oil palm, and other crops, relative to the baseline scenario. The outcomes of these two scenarios were then used to model the spatial distribution of oil palm plantations within each oil palm producing region up to 2030, using a spatially explicit and dynamic land use change model, operating at a spatial resolution of 10 3 10 km. This model accounts for competition between all land uses for space and can therefore simulate displacement of land uses. Comparison of the land use configurations resulting from these two scenarios provides an indication of the potential influence of ZDCs on oil palm and other land uses up to 2030. A stylized flowchart of the methodology is shown in Figure S2.

Computable general equilibrium modeling
Computable general equilibrium (CGE) models are systems of mathematical equations that describe economies as a whole and the interaction among their parts. 49 They are based on the premise of market-clearing, which occurs when the aggregate supply of goods and services equals aggregate demand. Producers are assumed to choose levels of input and output that minimize costs and consumers are assumed to maximize their utility subject to budget constraints. Policy changes can be simulated by changing variables of the model, which leads to a reallocation of labor, capital and land across sectors and geographies, until the system reaches equilibrium again. Responsiveness of producers and consumers to changes in relative prices and income is simulated by empirically calibrated elasticities of supply and demand.
In order to model the global market-mediated effects of ZDCs from 2014 up to 2030, we employed the most recent version of the Global Trade Analysis Project with differentiated Ago-Ecological Zones (GTAP-AEZ) CGE model 34 using the most recent GTAP database. 50 The database describes the world economy in 2014, disaggregated into 65 different sectors and 141 different regions, which we aggregated in 37 distinct regions for computational reasons (see Figure S2). Given that coverage of ZDCs across industries and regions was far from 100% in 2014, our model experiment is hypothetical as it does not inform on the actual effectiveness of ZDCs. Instead, the model experiment is designed to capture the potential effects that could be induced by ZDCs if they were fully implemented across industries and regions. To account for likely macroeconomic trends between 2014 -2030 that are unrelated to the adoption of ZDCs, we constructed a baseline or BAU scenario for the year 2030 that is largely consistent with the second Shared Socioeconomic Pathway (SSP), a scenario where socioeconomic trends broadly follow their historical patterns. 51,52 Following Johnson et al., 34 we used growth rates of real GDP, capital stock, population, unskilled and skilled labor from ECONMAP (v2.4) 53 which are calibrated against the SSP2 scenario. Sector specific productivity growth for ruminants and non-ruminants are taken from Ludena et al. (2007). 54 Due to lack of estimates for forest sector productivity growth, agricultural productivity growth figures are also taken from Ludena. Following Chateau et al. (2020), 55 a 2% productivity growth gap between manufacture and service sectors is also imposed. Finally, to account for baseline growth rates in crop yields, we used data from FAOSTAT (2021) and extrapolated yields up until 2030 based on the observed linear growth rate between 1998 -2014 (given the reference year of the GTAP database) for each agricultural commodity group within each region in the GTAP database. To avoid implausibly high yield increases, we scaled these projections using the global yield projections from the OECD-FAO Agricultural Outlook. 56,57 A special feature of GTAP-AEZ is that, in contrast to the standard GTAP model, 58,59 each region is further disaggregated into spatially heterogeneous land endowments or AEZs. These AEZs are determined based on climatic zones (tropical, temporal and boreal) and the length of the growing period (6 x 60-day intervals, see Lee et al., 2005 60 ), resulting in 18 different AEZs that may intersect multiple countries (see Figure S4). Within each region-specific AEZ (337 in total), the supply of available land is modeled using an empirically calibrated asymptotic curve specifying the relationship between land supply and the real land rental rate. 40,41 The land supply curves are predicated on the assumptions that expansion occurs within the most productive areas that are still available and that the land rental rate monotonically increases whenever the supply of available land decreases (see Figure S6). The maximum area available for expansion within each region-specific AEZ is constrained by a prespecified land supply asymptote. As land supply is critical in our scenarios, spatial differences in land availability need to be accounted for in the most accurate way. Therefore, we updated GTAP-AEZ by including new estimates of the land supply asymptotes for each region-specific AEZ, thereby better accounting for the constraints to land availability defined by earlier studies. 61 iScience Article A challenge for our model experiment is that the most recent GTAP database aggregates all oil crops (including oil palm) into one sector. Hence, to distinguish between oil palm and other oil crops in our analysis of GTAP outcomes, we used data from FAOSTAT (2021)

Land supply asymptotes
To estimate the extent of area that is potentially available for cropland expansion within each region-specific Agro-Ecological Zone (AEZ), we took three steps (see Figure S8 for a flowchart of the methodology). We first identified areas where cropland could possibly expand given certain biophysical, socio-economic, and institutional constraints through a residual approach. This means, we computed the total area within each region-specific AEZ after excluding areas already under cultivation, areas biophysically unsuitable for cropland cultivation, legally protected areas, rough terrains, and urban areas (see Table S1). To harmonize the different input data, all data were resampled to a 1 km 2 grid.
We used data from the European Space Agency Climate Change Initiative 63 ; for the year 2014 to identify existing cropland areas and urban areas. Areas classified as ''Cropland, rainfed'', ''Cropland, irrigated or post-flooding'' and ''Urban areas'' were assumed to be unavailable for expansion. The first two of these three classes do not cover all cropland areas as some cropland areas in the ESA-CCI database are subsumed under two mosaic classes (''Mosaic cropland (>50%) / natural vegetation (<50%)'' and ''Mosaic natural vegetation (>50%) / cropland (<50%)''). Following Liu et al. (2018), 64 we therefore assumed a cropland fraction of 58% and 38% for these two mosaic classes, respectively. To identify areas biophysically unsuitable for overall cropland expansion, we used data from Zabel et al. (2014). 65 These data incorporate suitability estimates for the 16 most important food and energy crops and account for the potential impact of climate change under SRES A1B conditions. Under this scenario, economies will rapidly grow, population growth will be small and there will be a rapid introduction of new and more efficient technologies. Areas classified as unsuitable were assumed to be unavailable for expansion. Furthermore, we excluded areas that are legally protected. 66 Although it has been recognized that encroachment may still occur within these areas, 67 it is unlikely that these areas represent future hotspots of agricultural expansion. 68 Finally, we excluded rough terrains (here defined as areas with steep slopes) as these areas are less likely to be developed. 69 The extent to which rough terrains can be cultivated depends, however, on the degree of agricultural mechanization, 70 which varies across space. Therefore, in an attempt to control for the spatially varying levels of agricultural mechanization, we took a region-specific approach to identify areas that are too steep to convert into cropland. Within each of the 37 regions in the GTAP database, we applied a slope threshold, based on the top 5% slope values within existing cropland areas. Slope data were sourced from Lloyd (2016). 71 As this resulted in unreasonably pessimistic slope thresholds in areas that are dominated by flat terrains, we imposed a minimum threshold of 10 degrees. To account for tiny landscape features that inhibit cropland expansion but cannot be detected at a 1 km 2 resolution (e.g., roads, rocky outcrops, water infrastructure, hedgerows, buildings etc.) we multiplied the extent of the remaining area of each regionspecific AEZ by 0. 85 iScience Article either extremely high or low cropping intensities, see Figure S9), we constrained the cropping intensities to lie in the range 0.5-2.0. Formally, this means we computed the area available for cropland expansion within each region-specific AEZ as follows: where Available area i denotes the available area within region-specific AEZ i in hectares; TA i denotes the total terrestrial area in hectares; UA i denotes the total area unsuitable for cropland expansion in hectares; and CI i denotes the cropping intensity. Cropping intensities were estimated by computing for each regionspecific AEZ the ratio of the harvested area of all crops in the GTAP database 50 to the total extent of cultivated area, based on the 10 arc-seconds resolution European Space Agency -Climate Change Initiative (ESA-CCI) satellite-based land cover map for the year 2014 (see Table S1 of for further information). Finally, we summed the total harvestable area available for expansion and the total area harvested in 2014 to arrive at an estimate of the land supply asymptote for each region-specific AEZ (see Figure S10A).

Implementation of zero-deforestation commitments
Implementing ZDCs involves reducing the (forest) area available for cropland expansion and hence, shifting the land supply asymptote to the left (see Figure S6). There are large uncertainties regarding the spatial coverage of ZDCs as the uptake and specificity of ZDCs varies across individual firms, commodities, and regions. 75 Nevertheless, common criteria outlined in many ZDCs are the protection of High Conservation Value Forests (HCVF) and High Carbon Stock Forests (HCSF). At least 78% of the palm oil refining capacity in Indonesia and Malaysia is covered by commitments which include these criteria 17 and the uptake of similar commitments in other industries and regions is on the rise. 76,77 We used data from Leijten et al.
(2020) 25 on the likely spatial distribution of HCVFs and HCSFs to delineate areas covered by ZDCs (see Figure S11). We focused on a middle-of-the-road estimate, which includes all tropical peatland forests and all forests with at least two overlapping HCVF categories or with at least 75 tonnes of carbon per hectare (t C ha À1 ) if located in the tropics (see Leijten et al. for further details). Updated estimates of the land supply asymptotes that account for the coverage of ZDCs are presented in Figure S10B. To explore the robustness of our results, we present alternative results using the top and bottom range estimates of the likely spatial extent of HCVFs and HCSFs. The low estimate (i.e., least protective) is based on three overlapping HCVF categories with a 75 t C ha À1 threshold while the upper estimate only requires one HCVF category with a 35 t C ha À1 threshold.

Spatial land use modeling
We used the GTAP-AEZ output as input for a dynamic land use model to simulate the spatial dynamics of oil palm and other land systems within all oil palm producing regions.

CLUMondo
CLUMondo is a spatially explicit and recursive dynamic land change model that can be used to simulate the spatial dynamics of land systems. 47, 78 It is an open-source model written in C++ that can be downloaded from: https://www.environmentalgeography.nl/site/data-models/models/. One of the unique features of the model is that it incorporates multifunctional land systems, i.e., land systems that provide multiple types of goods and services. To simulate future land use configurations, the model uses an iterative procedure where grid cells are allocated to the land system with the highest transition potential constrained by the overall regional production levels as determined by CGE models. Transition potentials represent the potential of a grid cell to be converted into a certain land system and are calculated for each individual land system though the following formula: where Ptrans t;i;LS represents the overall transition potential to land system LS in grid cell i at time t, Ploc i;t;LS represents the local suitability of land system LS, Pres LS represents the conversion resistance of land system LS, proxying conversion and investment costs as well as cultural attachments to the current land use, Pcomp t;LS represents the relative competitive advantage of land system LS, determined by the demand for the different products/services delivered by the land system and Pneigh i;t;LS denotes the neighborhood influence, representing agglomeration processes, for example due to economies of scale (2.5.5 To simulate the likely future spatial distribution of oil palm (given the GTAP-AEZ projections), we first created a 10 3 10 km land use map for the year 2014, the reference year of the GTAP-AEZ database. The land use map was primarily based on the 2014 ESA-CCI land cover map, 63 which we reclassified into 8 distinct land cover classes (see Table S2). To distinguish oil palm producing systems from other systems, we overlaid the resulting land cover map with a remotely sensed oil palm map from Descals et al. (2021), 79 aggregated to a 10 3 10 km resolution. We classified all grid cells with at least 50% oil palm cover as oil palm areas. Grid cells with 10-50% oil palm cover were classified as oil palm mosaics, which in addition to oil palm, produce a variety of other agricultural commodities (see below). GTAP regions with no oil palm areas or oil palm mosaic (according to the Descals map) were excluded from the analysis.
Demand for all other agricultural commodities is supplied by either cropland areas or cropland mosaics (which, like oil palm mosaics, represent a mixture of different land uses). As these two cropland categories vary in their agricultural intensity across space, we distinguished between high yielding systems, medium yielding systems, and low yielding (or extensive) systems. To classify cropland areas and cropland mosaics into these three different yielding systems, we overlaid the land use map with spatially explicit data (10 3 10 km) on the yields of all agricultural crops from the International Food Policy Research Institute (2019) and divided the data into quantiles. More specifically, we computed the sum of all yields within each pixel and selected the three highest quintiles (i.e., the three highest cut-off values that divide the distribution into five equal parts) to differentiate between high yielding, medium yielding and low yielding systems, respectively.
Finally, to distinguish grasslands that are intensively grazed from other grasslands, we used spatially explicit data (10 3 10 km) on the distribution of cattle in 2010 from Gilbert et al. (2018). 80 Grasslands within the top octile (12.5%) of the cattle density distribution were assumed to be intensively grazed. Figure S12 shows the resulting land use map. GTAP-regions with no or barely any oil palm area were excluded from the CLUMondo simulations. This includes China and Madagascar, which both accounted for less than 0.23% of global oil palm area in 2014. 1

Local suitability
Local suitability (Ploc i;t;LS ; see Equation 2) of the different land systems was assessed based on empirically quantified relations between land use patterns and several explanatory variables, using a logistic regression analysis for each land system within each oil palm producing region. Table S3 lists all 44 predictor variables used in the regression analysis. Model selection was done through a backward model selection procedure based on Akaike's Information Criterion. To avoid multicollinearity, predictor variables with intercorrelations exceeding 0.8 were removed from the analysis. To allow for dynamic changes in local suitability, we used spatiotemporal predictions of total and rural population density up to 2030 from CIESIN (2018). 81 The data were linearly interpolated for years for which no predictions were available.

Conversion resistance
For each land system in each region, a parameter was specified that captures the degree to which it is resistant to any type of land system conversion (Pres LS , in Equation 2. Land systems that involve large capital investments are typically more resistant to change. 82 Oil palm plantations typically require a large amount of capital investment, especially since they only start bearing fruit in the third or fourth year after establishment. 83 For that reason, we assumed that oil palm plantations tend to be relatively resistant to change, whereas land systems requiring little capital investment (e.g., extensive croplands) were assumed to be more reversible (less resistant to change once established at a location) and thus assigned a lower conversion resistance factor (see Table S5). In addition to conversion resistance, the spatiotemporal dynamics of the simulation are also influenced by a prespecified land conversion matrix that stipulates the type of land system conversions that are allowed and how long it takes before such conversions may take place. Given the time lag between the establishment of an oil palm plantation and the yielding period, the land conversion matrix was specified such that oil palm plantations can only be established after three years with stable land use. Finally, we also assumed, based on the estimates of Liebsch et al. (2008), 84  iScience Article each land system, we assigned a random number to all pixels indicating the number of years that the land system is already in place. All other land conversion rules developed for this analysis can be found in Tables S4 and S5.

Relative competitive advantage
The relative competitive advantage (Pcomp t;LS in Equation 2) of a land system in a certain region is determined through an iterative procedure in which the provided goods and services are compared to the total demand for those goods and services in the same region. Drawing on the GTAP-AEZ projections, we distinguished between five different types of demand: forestry, oil palm, other agricultural commodities, pastureland and urban areas; each of which is supplied by at least one land system (see Table S6). As our GTAP-AEZ results only provide predictions for the year 2030, we interpolated the annual changes in demand between 2014 -2030 based on the compounded growth rate of each type of demand. To ensure consistency with the GTAP projections, both in terms of area allocated to crops as well as to production volumes, the allocation was constrained by both the hectarage and the production volume. Given these constraints, the land allocation procedure determines the area and distribution of production in terms of low, medium and high intensity land systems. Since GTAP-AEZ does not solve for the amount of urban land in a region, we used global projections of future urban land expansion for the years 2020 and 2030 that are consistent with the SSP2 scenario from Chen et al. (2020). 85 Future expansion patterns of urban land were assumed to be independent of the implementation of ZDCs. To determine the quantities of the different goods and services that each land system provides per unit area, we took two steps. First, we used the aggregate demand values for the year 2014 to compute the average yield per grid cell across all the yielding land systems. Second, to calibrate the relative yields of the different land systems, we used several spatially explicit datasets on yields or production area ratios from the International Food Policy Research Institute (2019) 86 and computed the relative yields of the different systems. Spatially explicit data on production forests and grazing grasslands were sourced from Schulze et al. (2019) 87 and (Gilbert et al. (2018). 80 Finally, to attribute changes in land use to changes in oil palm demand, we compared the outcomes of our BAU and ZDC CLUMondo runs with two counterfactual scenarios (one for each scenario) that keep demand for oil palm constant until 2030 but are otherwise equal to the BAU scenario. This approach ensures that the oil palm-driven changes in land use can be distinguished from the other nonoil-palm driven changes in land use.

Neighborhood influence
Land use changes are often influenced by the spatial configuration of land use in neighboring areas. For example, it has been found that the strongest determinant of oil palm expansion in Malaysia is proximity to previously existing plantations. 88 This is especially important for oil palm as several plantations often deliver to one single mill, which has to be in close proximity. To account for such neighborhood influences (Pneigh i;t;LS ; in Equation 2), we used a 3 x 3 kernel function to adjust the transition potential in each grid cell for each land system, depending on the land system configuration in neighboring grid cells. The magnitude of the neighborhood effect was determined by a set of weight factors (see Table S7) and the fraction of the neighborhood that is occupied by each land system. Given that some land systems produce multiple goods (e.g., oil palm mosaics produce both oil palm and other crops), we scaled the weight factors depending on the average area composition of each land system. Finally, in light of the empirically supported theory that intensification is more likely to occur if land availability is scarce, 89-91 we followed Van Asselen and Verburg (2013) 47 by implementing a function that promotes cropland intensification under limited land availability and extensification under high land availability. Land availability within each 3 x 3 kernel was measured based on the land supply asymptotes described above.