Land-cover and land-use change trajectory hopping facilitates estate-crop expansion into protected forests in Indonesia

Protected areas (PAs) have been regarded as a critical strategy to protect natural forest (NF) and biodiversity. Estate-crop expansion is an important driver of deforestation in Indonesia. Yet, little is known regarding the temporal dynamics of PA effectiveness in preventing estate-crop expansion into NF. We employ Cox proportional hazard models and their extensions to characterize the dynamics of estate-crop expansion into NF in Indonesia during 1996–2015. The results show that PA effectiveness in Sumatra decreased over time and became insignificant in 2012–2015. A multistate modeling analysis shows that hopping in land-cover and land-use change (LCLUC) trajectories with shrub and/or bare ground as intermediates has decreased PA effectiveness and facilitated the expansion. Pre-venting LCLUC trajectory hopping becomes crucial to biodiversity conservation because it tends to occur at lowland forest, diminishing natural habitat area and increasing NF isolation.


INTRODUCTION
Protected areas (PAs) are regarded as crucial for reducing deforestation and conserving biodiversity, but their effectiveness remains disputed. At the national level, PAs have been shown to experience lower deforestation rates than unprotected areas (Andam et al., 2008;Geldmann et al., 2013). However, PA effectiveness varies widely over space depending on anthropogenic pressures and resource availability for protection (Geldmann et al., 2019;Jones et al., 2018). Land-use change is the dominant factor for biodiver- Wilcove et al., 2013). As global demand for estate-crop commodities, especially oil palm products, grows (Purnomo et al., 2020), estate-crop expansion is expected to continue, driving forest loss (Austin et al., 2019;Taheripour et al., 2019;Xin et al., 2021) and increasing pressure on PAs in the country (Curran et al., 2004;Poor et al., 2019). Evaluating PA effectiveness in preventing estate-crop expansion into natural forest (NF) is critical for improving conservation efforts. However, existing evaluations are sparse and limited to the province scale (Poor et al., 2019), do not account for variations in estate-crop expansion patterns across different islands (Xin et al., 2021), and have limited generalizability. This study quantifies the PA effectiveness across Sumatra and Kalimantan, where over 90% of Indonesia's estate-crop expansion occurred (Xin et al., 2021).
Partially attributed to rapid estate-crop expansion, Indonesia has experienced the largest increase in human pressure since the 21st century, even within PAs (Geldmann et al., 2019;Graham et al., 2021). Understanding how PA effectiveness changes over time in response to increasing human pressure is vital for cost-effective resource allocation and maximizing conservation benefits. Although the importance of examining the time-variant PA effectiveness is recognized (Geldmann et al., 2015), empirical research investigating how the PA conservation outcomes evolves over time remains scarce in the literature. Existing research focuses on evaluating PA effectiveness in reducing deforestation, which commonly relies on techniques such as matching methods Geldmann et al., 2013Geldmann et al., , 2019Graham et al., 2021), multifactor regressions (Brun et al., 2015), and Before-After Treatment-Intervention (Gaveau et al., 2007;Shah & Baylis, 2015). However, these methods do not account for temporal details, and are thus unable to capture the temporal dynamics of PA effectiveness.
Survival analysis, which is specialized for event data modeling and explicitly deals with the occurrence and timing of events, is able to address the temporal complexities in longitudinal studies (Wang et al., 2013). Extending by including covariates interacted with time, it can address the possible time-dependent effects under rapid land changes (Chen et al., 2016). We employ the extended survival analysis to address the potential changes of PA effectiveness in preventing estate-crop expansion into NF over time. Standard survival analysis concentrates on the timing of a single event of interest; however, there are many examples in which a subject may experience a variety of intermediate events during the study period. This is especially true for dynamic landscapes in Sumatra and Kalimantan (Jamaludin et al., 2022), where landcover and land-use change (LCLUC) is typically featured by a sequence of changes along the temporal trajectory . In several cases, deforestation is not immediately followed up with cultivation, but with degraded land of relatively low biomass as intermediate states (Jamaludin et al., 2022). Conversion of NF to estate-crop plantation could either occur directly (rapid conversion) or with trajectory hopping (NF → intermediate states → estate-crop plantation). Therefore, we adopt the multistate model, which is capable of jointly considering all initial, intermediate, and final states, to address the complexities along with LCLUC trajectories.
In this research, we assess the dynamics of PA effectiveness in preventing estate-crop expansion into NF over 1996-2015, with temporal details and varying LCLUC trajectories, in Sumatra and Kalimantan. The estimation results of the Cox proportional hazard models and their extensions indicate that PA effectiveness in Sumatra decreased over time and became insignificant in 2012-2015; the temporal trend was not significant in Kalimantan as Kalimantan is a latecomer with relatively abundant land resources. By applying multistate models, we demonstrate that PA effectiveness is reduced by LCLUC trajectory hopping with shrub and/or bare ground (BG) as intermediates in both Sumatra and Kalimantan.

METHODS
We acquired land-cover and land-use (LCLU) maps of 1996, 2000, 2003, 2006, 2009, 2012, and 2015 from the Indonesian Ministry of Environment and Forestry (MoEF) to analyze the LCLUC trajectories in Indonesia, focusing on Sumatra and Kalimantan. We reclassified LCLU classes to six categories-NF, shrub, dry agriculture (food crops except for paddy), estate crop (dominated by oil palm), BG, and others ( Figure 1)-and quantified the LCLUC among the categories in each 3-or 4-year period. The description of maps and processing methods are presented in Supporting Information (Note S1.1; Table S1). LCLUC occurred within 3(or 4)-year intervals are defined as direct conversions, since it usually takes 2-4 years to have sufficient plant growth to be detected by remote sensing (Austin et al., 2019). Sankey diagrams (Note S1.2) were used to visually display the LCLUC trajectories of deforestation. Deforestation, synonymous with forest loss, is defined as the transformation of forests to a nonforest state at the mapping scale. We use survival analysis to quantify the extent and timing of estate-crop expansion. We run the models on 1 × 1 km grid cells with over 90% NF coverage in 1996 (Note S1.3). We define the estate-crop expansion transition in a grid cell as a single event, with hopping to intermediate states (e.g., BG in NF → BG → estate-crop process) being treated as part of survival time. In this way, whether the estate-crop expansion into NF occurred directly or with intermediates makes no difference and both are included in the analysis. We used the Cox proportional hazard model to calculate the hazard rate-the risk of estatecrop expansion at a time of interest. The Cox proportional hazard model takes the following form: where we use the occurrence of estate-crop expansion transition as the dependent variable ℎ ( ). The term 0 is the baseline hazard. represents the protection status in the starting year of each period, and 1 ⋅ ⋅ ⋅ ik are a set of biophysical and socioeconomic variables selected to proxy the pressure of converting NF to estate-crop plantation (Note S1.4).
The Cox model assumes explanatory variables exert constant effects on hazards over time, which might be violated F I G U R E 2 Graphical representation of the multistate model. Nodes represent possible states and links represent possible conversions between states, where parameters that influence the hazard are indicated. The possible trajectories include (1)  in the process of rapid LCLUC. Therefore, we extend the model with time-dependent effects (Note S1.5), following Chen et al. (2016). We further check whether the timedependent effects have been addressed by the extended models successfully; if not, we apply strata models to these variables to see their effects in each period (Note S1.6).
We used multistate models to reveal the different risks of direct conversion from NF to estate crop and conversions with hopping to different intermediates, conditioned on the characteristics (e.g., , ) of each grid cell. We use a semi-Markov chain with estate-crop plantation as the final states, as shown schematically in Figure 2 (Notes S1.7 and S1.8).

LCLUC trajectories in Indonesia
Estate-crop plantation in Indonesia expanded from 56,900 km 2 in 1996 to 127,300 km 2 in 2015, with approximately half occurring at the expense of NF in 1996.
Conversions from NF to estate crop could either occur directly within one observation period (3-4 years) or take a longer process with hopping to degraded land, especially shrub and BG, as intermediates ( Figure S1; Note S2.1). A total of 55.50% of estate-crop expansion occurred as direct outcome of deforestation and 39.94% as indirect outcome of trajectory hopping via shrub or/and BG (Figures 3 and S1). Four major trajectories that together account for over 95% of all estate-crop expansion into NF are (1) direct transition from NF to estate crop, and three processes with trajectory hopping: Figure 3).

Estate-crop expansion into protected NF
In 1996-2015, appropriately 13,660 km 2 NF was lost in Indonesian PAs, with 6772 km 2 in Sumatra and 2432 km 2 in Kalimantan. Around 922 km 2 estate-crop expansion occurred in Sumatra PAs in 1996-2015, 71% of which were at the expense of NF in 1996. More than 90% of the direct conversion occurred in 1996-2000 ( Figure 4c). Estate-crop expansion into protected NF in Kalimantan was limited (around 87 km 2 ), with around 70% of the expansion occurring in 2012-2015 ( Figure 4d). In both Sumatra and Kalimantan, shrub was the largest sink of NF loss in PAs before 2012, while BG became the largest sink in 2012-2015 (Figure 4c,d). Once hopped to shrub, the risks of further conversion to estate crop increased by about 12 times (∼0.008 to ∼0.1) in Sumatra. The majority of BG established after 1996 was further converted to estate crop in the following period in both islands ( Figure S3).

Effectiveness of PAs
The results of Cox proportional hazard models demonstrate that estate-crop expansion into NF tended to occur at lowland forest with high oil palm attainable yield but low annual average temperature, at gentle slope, as well  as close to old plantation and palm oil processing mills (Table 1, Base model). The extended models reveal that the effects of infrastructure accessibility (e.g., accessibility to old plantations and processing mills) decreased over time and/or even changed directions in 2012-2015, while the effects of biophysical suitability (e.g., elevation, slope, attainable yield) were not time variant (Table 1, Model with time-variant effects; Note S2.2).
As trajectory hopping to shrub and/or BG before estatecrop expansion is prevalently observed, we use shrub and BG as the intermediate states of the multistate models ( Figure 2). With the multistate analysis, we find that 8-year shrub (8 years since deforestation to shrub) had the highest risks of further conversion to estate crop in Sumatra, and 3-year shrub and BG had the highest risks of further conversion to estate crop in Kalimantan ( Figure S2). Although  Note: t refers to time. *, **, ***, and **** stand for the significant level of 10%, 5%, 1%, and 0.1%, respectively.
conversion with intermediates made more share in recent years (Figure 4a,b), the effects of biophysical and socioeconomic drivers (except PA) did not vary much between the total and the direct estate-crop expansion into NF (Note S2.3). The patterns of estate-crop expansion into NF in PAs (Table S8), as well as deforestation to shrub and BG (Tables S9.3 and S9.4), generally mirror those across the whole islands (Note S2.4). The risks of estate-crop expansion into NF in both Sumatra and Kalimantan were significantly lower in PAs than outside PAs (Table 1, Base model); however, the extent of differences in Sumatra decreased after 2000 and became insignificant in 2012-2015 ( Figure 5). In Sumatra, PAs were highly effective in preventing direct deforestation to estate crop as the hazard level of PAs was only 2%−13% of that in non-PAs at the 95% confidence level and had no significant temporal trend ( Figure 5). The decreasing trend of PA effectiveness was present with hopping to shrub as an intermediate state ( Figure 5). After 2000, the hazard levels of estate-crop expansion into shrub (established after 1996) became higher in PAs as time went by and became significantly higher than that in non-PAs in 2012-2015.
Although the hazard levels of deforestation to shrub were significantly lower in PAs (48%−64% of the hazard levels in non-PAs), the differences were much smaller than direct deforestation to estate crop (2%−13% of the hazard levels in non-PAs). In Kalimantan, although PA effectiveness was significant in preventing the total estate-crop expansion into NF and showed no significant temporal trend, the effectiveness in preventing deforestation to shrub and BG was remarkably smaller than in preventing direct deforestation to estate crop ( Figure 5). The high risks of BG establishment inside PAs, especially in 2012-2015 (Figures 4d and 5), substantially contributed to forest loss in Kalimantan PAs. As PAs were not effective in preventing conversion from BG to estate crop, the land may eventually be utilized for estate-crop production.
PA effectiveness is significantly influenced by the edges ( Figure S4). NF at the boundaries of PAs had significantly higher risks of conversion to estate crop than NF completely inside PAs, although the extent of differences became insignificant in 2012-2015 in Sumatra. The decreasing trend in Sumatra was mainly due to the conversions with hopping to shrub as an intermediate state (Note S2.5).

DISCUSSION
Estate-crop expansion into PAs has been observed in many countries, especially in tropical regions critical for biodiversity conservation and climate mitigation (Castellanos-Navarrete, 2021; Riggio et al., 2019). Our study is the first F I G U R E 5 Effectiveness of protected areas (PAs) on estate-crop expansion into natural forest as revealed by the multistate regressions. The vertical axis is the exponent of the PA coefficient, exp( 1 ), indicating the hazard a grid would be subject to because of PAs in each 3(or 4)-year period. The horizontal axis shows different land-cover and land-use change (LCLUC) types. The first six represent the six direct conversions in the multistate process, and the seventh, Estate-Crop Expansion, represents the total estate-crop expansion into natural forest. Vertical bars correspond to 95% confidential intervals. For each transition type on each island, the vertical bars from left to right represent time period 1996-2000, 2000-2003, 2003-2006, 2006-2009, 2009-2012, and 2012-2015, respectively. Only transition types that PAs are significantly effective in at least one period are shown in the figure. The effectiveness of PAs in conversion from shrub to bare ground in Sumatra and conversion from bare ground to estate crop in Kalimantan was insignificant during the whole study period (p > 0.3 and p > 0.6); therefore, they are not shown in the figure. The dashed horizontal line is Hazard = 1. If the vertical bars fall completely under the dashed horizontal line, the grids inside PAs have significantly lower risks of experiencing the LCLUC than those out of PAs; if the vertical bars fall completely above the dashed horizontal line, the grids inside PAs have significantly higher risks of LCLUC than those out of PAs; if the vertical bars intersect with the horizontal line, the effects of PAs are not statistically significant.
to examine the effectiveness of PAs in curbing estate-crop expansion into NF across Sumatra and Kalimantan, the epicenter of palm oil production, contributing over 50% of the global supply (Xin et al., 2021). Unlike previous research that relies on snapshot data and does not account for temporal details (Brun et al., 2015;Graham et al., 2021), we employ a longitudinal study with survival analysis and multistate models to capture temporal aspects of PA effectiveness and LCLUC hopping. Our findings highlight the role of the trajectory-hopping mechanism, which involves degraded land as intermediates, in facilitating estate-crop expansion into NF and decreasing PA effectiveness.
Our research shows that estate-crop expansion risks in Indonesia are significantly lower inside PAs and away from boundaries, confirming that PAs effectively prevented deforestation Graham et al., 2021), especially in core areas (Poor et al., 2019;Watson et al., 2014). However, the time-variant survival analysis reveals decreasing PA effectiveness over time, with this trend more pronounced in Sumatra than in Kalimantan because Kalimantan is a latecomer with relatively abundant land resources (Austin et al., 2017;Xin et al., 2021). As land resources become scarcer, estate-crop expansion encroaches further into PAs and goes deeper into the central area. Our results align with existing literature, demonstrating that estate-crop expansion tends to occur at lowland forest areas, which are biophysically and socioeconomically suitable for estate-crop plantation (Pirker et al., 2016;Xin et al., 2021). Unfortunately, these areas are also known for their high levels of biodiversity richness (Kinnaird et al., 2003;Poor et al., 2019). Accessibility to infrastructure becomes less constraining over time, allowing estate-crop expansion into more remote area. This growing incursion into protected NF not only results in habitat loss but also isolates forest patches (Kinnaird et al., 2003) and introduces anthropogenic disturbances (e.g., infrastructure construction, use of fertilizer and pesticides), adversely affecting biodiversity (Barlow et al., 2016). Therefore, it should be a main concern of conservation activities. To maximize conservation benefits with limited management resources, it is vital to prioritize PAs in regions with scarce land resources but suitable for estate-crop plantation.
The multistate models show that the risk of direct conversion from NF to estate crop remains minimal in PAs; however, the effectiveness of PAs is significantly reduced by the LCLUC trajectory-hopping tactic. Forest loss to shrub and BG inside PAs closely resembles that outside of PAs (Laurance et al., 2012) and is a key driver of deforestation (Austin et al., 2019). PA effectiveness is notably smaller in preventing deforestation to shrub and BG compared to direct deforestation to estate-crop plantations. Once NF is transformed into shrub and/or BG, the risk of conversion to estate-crop multiplies, especially in recent years (Geldmann et al., 2019). Over half of the degraded lands established in Indonesia since 1990 are due to deliberate clearing (Parker, 2022), and this share would be higher if intentional fire were included (Harrison et al., 2009). The rapid conversion of degraded land to estate-crop plantation (risk is highest with 8 years for shrub and 3 years for BG) and limited natural regeneration (e.g., BG to shrub, shrub to NF; Figure S5) suggest a considerable proportion of the degraded land is being actively used as a land banking to facilitate estate-crop expansion into NF, especially in illegal areas (Xin et al., 2021).
Forest loss to degraded land is often seen as a natural process, while estate crop replacing degraded land is frequently overlooked in monitoring and evaluation, particularly in Indonesia's complex policy landscape (Astuti et al., 2022;Parker, 2022). Mobilizing different stakeholders in PA management and identifying trajectory-hopping mechanism are vital for forest conservation. The recent "One-Map Policy" is a promising start, as it clarifies land boundaries, establishes unified geospatial information data, and promotes interinstitutional collaborations (Astuti et al., 2022). However, given the economic focus of government policies (Astuti et al., 2022) and limited local involvement in PA management (Myers et al., 2017), nongovernmen-tal organization-led conservation coalitions that aggregate stakeholders from local to international level can play a key role in managing PAs and mitigating deforestation (Ruysschaert & Hufty, 2020). To limit forest loss to estate crop via degraded land, policies and enforcement measures should aim to prevent clearing and intentional burning of NF and prohibit sales of land and estate-crop products on recently established degraded land. In this connection, a near-real-time and open-access monitoring system can help detect the LCLUC timely, disseminate the knowledge, facilitate broad cooperation among stakeholders, and strengthen public oversight by visualizing the PA effectiveness (Tacconi et al., 2019;Taheripour et al., 2019).
Trajectory hopping not only reduces PA effectiveness but also contributes to estate-crop expansion into NF outside PAs in Indonesia. Although many producers commit to sustainable requirements (e.g., zero deforestation, RSPO, ISPO), deforestation to nonproductive land is sometimes used as a land-banking mechanism prior to certification and creates a "blind spot" in sustainable palm oil production (Astuti et al., 2022;Xin et al., 2021). Recent remote-sensing data indicate a decline in estate-crop plantation and deforestation after 2016, but this may be due to dropping palm oil prices rather than effective policies (Parker, 2022). The 2020 deforestation increase in Kalimantan, alongside palm oil price recovery, underscores the need for effective regulations to prevent forest loss, especially through LCLUC trajectory hopping.

A C K N O W L E D G M E N T S
We thank the Ministry of Environment and Forestry of Indonesia for creating and making available the land-cover and land-use data used for this study; the latest map may be viewed via the following map service: https://sigap. menlhk.go.id/sigap/peta-interaktif. Inputs for this study consisted of historical epochal land-cover and land-use information from 1990 through 2015. The study was made possible through the financial support of the National Aeronautics and Space Administration's (NASA) Land Cover and Land Use Change program (Grant number: 80NSSC18K0723).

D ATA AVA I L A B I L I T Y S TAT E M E N T
The sources of data that support the findings of this study are provided in Table S2.