Tipping Points in the Supply of Ecosystem Services of a Mountainous Watershed in Southeast Asia

Rubber plantations have expanded at an unprecedented rate in Southeast Asia in recent decades. This has led to a substantial decline in the supply of ecosystem services (ESS) and has reduced livelihood options and socioeconomic well-being in rural areas. We assessed the impact of two land use scenarios on the supply of ESS in a mountainous watershed in Xishuangbanna Prefecture, People’s Republic of China. We combined time-series data derived from spatially explicit ESS models (InVEST) with a sequential, data-driven algorithm (R-method) to identify potential tipping points (TPs) in the supply of ESS under two rubber plantation expansion scenarios. TPs were defined as any situation in which the state of a system is changed through positive feedback as a result of accelerating changes. The TP analysis included hydrological, agronomical, and climate-regulation ESS, as well as multiple facets of biodiversity (habitat quality for vertebrate, invertebrate, and plant species). We identified regime shifts indicating potential tipping points, which were linked to abrupt changes in rubber yields, in both scenarios at varying spatial scales. With this study, we provide an easily applicable method for regional policy making and land use planning in data-scarce environments to reduce the risk of traversing future TPs in ESS supply for rubber producing land use systems.


Introduction
Ecosystem services (ESS) are defined as the goods and benefits people obtain from ecosystems [1].The capacity of ecosystems to sustainably deliver ESS can change and is often threatened by increasing anthropogenic pressure.This pressure includes increased climatic variability, land use change, and the overexploitation of resources [2][3][4].The intensity of external pressure that a system can absorb without changing its characteristic structure and functions is defined as resilience [5,6].The faster the system returns to its former equilibrium state after disturbance, the more stable it is considered to be [5].Socioecological systems are characterized by complex interactions between ecosystem properties and social dynamics.Their resilience is defined as the ability of the socioecological system to sustain a set of ESS under changing environmental and management conditions [7].This ability to generate ESS may shift from desired to less desired states as a result of increasing external pressures, which often act synergistically [8].Resilience may then be compromised and the system reaches a tipping point (TP).Shifts from one system state to another may not be desirable, costly to revise, or irreversible [6,9].
Using the concept of TPs to analyze socioecological systems has gained considerable momentum in academia in recent years [10].Milkoreit et al. [10] deduced the concepts most frequently used for defining TPs across scientific disciplines (natural science, social science, and socioecological systems).In their review, the majority of literature on land use change and socioecological systems involve the concepts of "Multiple Stable States" and "Abruptness".Concepts involving "Feedbacks" are frequently used in socioecological system research, whereas land use change research often involves the concept of "Irreversibility".The term "Tipping Point" is increasingly used in research on socioecological systems to describe phenomena comparable to the better known term "Regime Shift" in ecological systems [10].We refer to a general definition for TPs from van Nes et al. [11]: "Any situation where accelerating change caused a positive feedback [that] drives the system to a new state".To date, only few studies have looked at TPs in ESS supply.
Zhang et al. [12] analyzed time series of 55 social, economic, and ecological indicators, ranging from soil stability, air quality, and water quality to livestock and various crop yields, in order to assess the sustainability of ESS in Eastern China.They found that the analyzed regional socioecological system already passed a TP in the 1970s, coinciding with China's shift from a planned economy to a market economy.After China's accession to the World Trade Organization in 2001, the system is now heading towards a new steady state.Such TPs are of concern, as passing critical thresholds and the associated changes can affect the ability of a system to provide ESS and lead to disastrous consequences for human and ecological societies that depend on them [13].
The concept of a "safe operating space for humanity" was proposed by Rockström et al. [14,15].They interlinked nine critical boundaries of biophysical processes at a planetary scale that frame the safe operating space.Within this space, the risk for unpredictable or catastrophic changes is low.Raworth [16] extended the planetary boundary concept and integrated human well-being into the framework, which has become known as the Oxfam Doughnut.In the Oxfam Doughnut, the borders of the safe operating space are defined by the social foundation and the environmental ceiling.Falling below the social foundation or breaching the environmental ceiling reduces human well-being [16].On a planetary scale, three of the environmental boundaries have already been transgressed (climate change, rate of biodiversity loss, and the nitrogen cycle) and millions of people's living standards fall below the social foundation (e.g., food, energy, income) [16].
In contrast to the planetary scale of these approaches, management of resources and decision making mainly takes place at a regional or local scale.Dearing et al. [17] propose a method to transfer the Doughnut concept into a framework of regional boundaries.They describe four types of time-series analyses to define the boundaries of a safe operating space at a regional scale: (1) linear trends with environmental limits set by scientific, public, or political instances (e.g., urban air quality regulation); (2) the envelope of variability defined by long-term, normal fluctuations of a system (e.g., climate system of planet earth); (3) analysis of systems that crossed a critical threshold in the past (e.g., eutrophication of a lake); and (4) the identification and analysis of early warning signals (e.g., "critical slowing down").At regional levels, processes of deforestation and agricultural expansion are among the most important drivers for regime shifts that impact the supply of ESS in rural areas [18].These processes and their impact on ESS can be simulated and assessed with spatially explicit models [19].Modeling results can then be used to delineate the boundaries of local safe operating spaces in order to predict and reduce the risk of traversing TPs in the future.This is of particular importance in areas where people's current and future livelihoods are directly linked to a sustainable supply of ESS, such as rural agricultural areas.
In Montane Mainland Southeast Asia (MMSEA), the recent intensification in the cultivation of cash crops has led to the demise of traditional swidden farming systems [20].Major reductions in livelihood options, socioeconomic well-being, and the supply of ESS have been the result of such large-scale land use conversions [21][22][23].A prime example of these conversions is the introduction and extensive expansion of rubber tree (Hevea brasiliensis) plantations in many parts of MMSEA.The Prefecture of Xishuangbanna, located in Southwestern China, has seen increases in area covered by rubber plantations, from 4.5% in 1988 to 22.2% in 2010 [24].The rapid expansion of rubber plantations has been linked to increased soil erosion, loss of biodiversity, changes in carbon dynamics, and changes in the hydrological cycle [25][26][27][28][29][30].However, rural areas in MMSEA often lack the amount of long-term and high-quality data needed to advance the scope of multidisciplinary modeling.To date, no research has been done on deriving potential TPs in ESS in rubber-dominated land use systems.
As rubber expansion has proceeded at an unprecedented rate in MMSEA in recent decades, there is an urgent need for methods to assess the future supply of ESS from rubber-producing land use systems, even in data-scarce environments.We propose a method for identifying TPs in the supply of ESS.For this, we use results from spatially explicit ESS models such as InVEST (Integrated Valuation of Ecosystem Services and Trade-Offs) in combination with a data-driven, sequential algorithm, originally developed to detect regime shifts in climate variations.The analysis includes hydrological, agronomical, and climate-regulation ESS, as well as multiple facets of biodiversity (habitat quality for vertebrate, invertebrate, and plant species) for two rubber-related land use scenarios, simulated at the watershed scale over a 25-year period.We aim at providing a simple tool for regional policy making and land use planning to reduce the risk of traversing future TPs in ESS for rubber producing land use systems.This method can also be transferred to other comparable land use systems outside of MMSEA.

Research Area
The research area is the Naban River Watershed National Nature Reserve (hereafter referred to as Nabanhe Reserve).It spans an area of roughly 271 km 2 and is located in Xishuangbanna Prefecture, Yunnan Province, Southern China (22 • 08 N 100 • 41 E).Xishuangbanna is characterized by a subtropical climate, strongly influenced by monsoon cycles and a distinct wet season from May to October [31].The mean annual temperature lies between 18 and 22 • C and annual average precipitation amounts to 1100-1600 mm [27].The elevation ranges from 500 to 2300 m.a.s.l., with rubber plantations being located in lowland valley bottoms and, more recently, expanding into higher altitudes [32].The Nabanhe Reserve is part of the Indo-Burma biodiversity hotspot [33] and is not only rich in floral and faunal species diversity but also in terms of the multiple ethnicities living in the area [34,35].The largest share of area in the Nabanhe Reserve is covered by natural or seminatural tropical mountainous rainforest (~60%).Other land uses include tea plantations and bushland (~11%), rubber plantations (~9%), bamboo forests (~6%), annual crops such as maize and sugar cane (~6%), paddy rice (~4%), and small patches of perennial crops (~1%) (Figure 1a).

Scenario Definition
The scenarios used for this study were developed as part of the Chinese-German project SURUMER (Sustainable Rubber Cultivation in the Mekong Region) [36].The project had the objective to develop sustainable land use strategies for rubber cultivation through a close interaction of science and practice [37].A team of relevant stakeholders in collaboration with a multidisciplinary team of scientists developed several future land use scenarios [38].We used two of these scenarios to assess potential TPs in the supply of ESS in the research area: (1) The Business-As-Usual scenario (BAU) and (2) the Balanced-Trade-Offs scenario (BTO).In BAU, past rubber expansion in Xishuangbanna [39] and the Nabanhe Reserve [28] was linearly continued at an annual rate of 2%, independent of the suitability of the land.In BTO, rubber expansion was restricted to areas suitable for rubber cultivation.These were defined as below 900 m.a.s.l. and with slopes below 23 • .These restrictions were chosen based on reduced rubber yields at high altitudes [40] and the aim to reduce soil erosion on sloping areas [27].BTO additionally included the establishment of riverine buffer zones around the main rivers in Nabanhe Reserve.These consist of secondary forest vegetation and are 30 m wide.Similar reforestation measures were introduced in the initial years of BTO as water protection zones around water sources for domestic use.The land use map of 2015 (Figure 1a) served as the initial condition for both scenarios.The scenarios have been described in more detail in Thellmann et al. [41].Land use changes were simulated for 25 years, from 2015 to 2040.The focus in Thellmann et al. [41] was on deriving realistic land use change scenarios and information on ESS preferences and evaluation through a close interaction with stakeholders.Here, we refrained from social or economic valuation for ESS and focused on analyzing the biophysical changes in ESS supply for two contrasting land use scenarios across varying spatial scales and initial land use conditions.Figure 1c-f shows the spatial extent of land use changes introduced in both scenarios for the final year of the simulation (2040).
Sustainability 2018, 10, x FOR PEER REVIEW 4 of 15 [41] was on deriving realistic land use change scenarios and information on ESS preferences and evaluation through a close interaction with stakeholders.Here, we refrained from social or economic valuation for ESS and focused on analyzing the biophysical changes in ESS supply for two contrasting land use scenarios across varying spatial scales and initial land use conditions.Figure 1c-f shows the spatial extent of land use changes introduced in both scenarios for the final year of the simulation (2040).

Ecosystem Service Assessment
Ostrom proposed a framework for analyzing the sustainability of a socioecological system [42].We used this framework as a basic guideline to identify relevant variables: (1) resource systems-the Nabanhe Reserve; (2) resource units-quantification of ESS supply from spatial models; (3)

Ecosystem Service Assessment
Ostrom proposed a framework for analyzing the sustainability of a socioecological system [42].We used this framework as a basic guideline to identify relevant variables: (1) resource systems-the Nabanhe Reserve; (2) resource units-quantification of ESS supply from spatial models; (3) governance system-future land use plans based on past land use change and governmental land use plans; and (4) users-feedback from local stakeholders to identify locally relevant ESS.
We assessed a set of four ESS using InVEST (Integrated Valuation of Ecosystem Services and Trade-Offs, Version 3.3.3., The Natural Capital Project: Stanford, CA, USA) [43] and used a self-developed model for assessing rubber yields, as rubber is the most important provisioning service in Nabanhe Reserve.These ESS were chosen based on several discussions with stakeholders, including village heads, farmers, prefecture administration, and politicians at the provincial level [38,41,44].The InVEST submodels we used were: (1) habitat quality (as a proxy for biodiversity), (2) carbon storage and sequestration, (3) water yield, and (4) sediment retention (as proxies for regulating ESS).The code of InVEST was modified to iterate model runs in 25 steps for the 25 years of land use changes outlined in the scenario description.Otherwise, we implemented InVEST according to the developers' user manual [43] and published literature [30,[45][46][47] and used a substantial pool of field studies for model parameterization [41].Environmental input variables such as long-term mean annual precipitation and soil properties were kept constant throughout the simulation period [48,49].This was done in order not to bias the effect of land use changes on the model results with changes in other parameters.Note that we considered a decrease in sediment export as beneficial.On the other hand, due to water scarcity in the dry season, decreases in annual water yield were considered to be unfavorable.The rubber yield model was based on local survey data considering average rubber yields with regard to plantation altitude and plantation age [50,51] and was performed with ArcGIS (Version 10.3.1).
We normalized the biophysical model results so that the values of each ESS are set to 1 in the initial year of the simulation.ESS values for the following years in the simulation were then calculated in relation to the initial year.We calculated the arithmetic mean of the five normalized ESS values for every year.This we refer to as the ESS z-score.The time series of the ESS z-score was then used as an input for the TP identification algorithm described in Section 2.4.A detailed description of the conceptualization, parameterization, and results of the InVEST submodels and the rubber yield model can be found in Thellmann et al. [41] for the entire Nabanhe Reserve.Here, we focused on three subwatersheds in the Nabanhe Reserve in more detail: SW1, SW2, and SW3 (Figure 1).We selected these subwatersheds as (a) data was available for calibrating the submodels for water yield and sediment export (in SW1) [27], and (b) they feature varied conditions regarding their spatial extent, initial land use conditions, and land use change trajectories in the simulated scenarios.Therefore, they allowed us to analyze the provisioning of ESS in relation to different rubber expansion rates as well as different initial land use conditions.The initial land use conditions are listed in Table 1.SW1 is the smallest subwatershed, at 6.92 km2 in area.SW2 spans 27.97 km 2 and SW3 has an extent of 21.34 km 2 .

Identification of Tipping Points
In order to identify potential TPs in the time series of the ESS z-score in our scenarios, we adapted the method of Zhang et al. [12], who used long-term time series of economic and environmental data to analyze the past supply of multiple ESS in eastern China.The algorithm we used to statistically identify TPs was developed to detect regime shifts in time series of climate variations by Rodionov [52].We chose the Rodionov algorithm (henceforth referred to as the R-method) as it is capable of identifying TPs even in the presence of a background trend and outliers in the data.The R-method does not require prior knowledge of the timing of a potential TP as it is an exploratory data-driven analysis.We refer to the periods separated by TPs as regimes.The input for the R-method was the annual time series of the ESS z-score for both scenarios in the simulation period between 2015 and 2040.Years with potential TPs are indicated by the regime shift index (RSI).A positive RSI signifies a shift to a new regime, where the mean ESS z-score for the entire regime is greater than the mean of the previous regime.The opposite is true for a negative RSI.The higher the absolute value for RSI, the more distinct is the subsequent regime from the previous one.To prevent the misidentification of very short regimes, the R-method was set to a moving average window of half the length of the simulation period (cutoff = 12.5), with significance of p = 0.01 and Huber's weight parameter w = 2. Huber's weight parameter defines to which extent outliers are taken into account in the algorithm of the R-method.

Biophysical Model Results of InVEST and Rubber Yield
Table 2 shows the results of the InVEST submodels and the rubber yield model for the initial condition (2015) and the final years (2040) of BAU and BTO for the entire Nabanhe Reserve as well as each subwatershed.In both scenarios, the simulated land use changes show the same general trend in affecting the supply of ESS, regardless of the spatial scale of the analysis.These trends are: (1) a decrease of the regulating ESS and habitat quality in BAU and (2) increases of habitat quality and the regulating ESS in BTO in the Nabanhe Reserve and the subwatersheds.Rubber yields increase in both scenarios but are highly variable throughout the simulation period as plantations shift in and out of their economic life cycle in different parts of the Nabanhe Reserve.The final year of BAU shows an increase of 1.13 × 10 6 kg rubber yield compared to the initial condition.A slightly lower increase can be seen in BTO (1.07 × 10 6 kg).Even though there is a decrease of exported sediments in BTO, there are only minor changes in water yield.The land use changes introduced in both scenarios have only minor effects on water yield (±1%, which is masked by rounded numbers in Table 1).

Land Use Change and Tipping Points in ESS Supply
Figure 2 shows the temporal changes in the normalized supply of ESS as well as rubber-and forest-related land use changes for both scenarios.Additionally, the ESS z-score, the arithmetic mean of the five normalized ESS values (light blue line in Figure 2) calculated on an annual basis, is shown.Years with potential TPs are depicted as purple columns in Figure 2 (RSI: regime shift index).In both scenarios, we identified potential TPs in all case study subwatersheds, all of which are significant with p > 0.01.
In the BAU scenario, for each subwatershed, as well as the entire Nabanhe Reserve, habitat quality and the three regulating ESS decrease in response to the rubber expansion and the loss of forest areas.Rubber yields are increasing in the initial years of the simulation, then remain relatively constant until 2036, and see a steep decline at the end of the simulation.The year of a potential TP in total ESS supply is 2036 for NR, SW1, and SW2.For SW3, a potential TP is found in 2034.The trajectory for sediment retention in SW3 shows a more pronounced decline in comparison to the sediment retention trajectory in the other subwatersheds.SW2 has the smallest share of rubber areas in the initial year of the simulation (7.6%), and with that, also the lowest total annual rubber yield in comparison to the other subwatersheds.As the rubber yield trajectories are all relative to the initial condition (set to 1), a lower initial rubber yield results in a steeper increase of the trajectory, as can be seen for the rubber yield in SW2 (Figure 2e).SW2 is the only subwatershed in the BAU scenario where a positive TP can be observed (2024).
In the BTO scenario, the reforestation measures and restricted expansion of rubber plantations led to increases in the regulating ESS and habitat quality.Rubber yields follow trajectories that are comparable to those seen in the BAU scenario, although with lower total values.In the BTO scenario, a positive RSI indicates a potential TP for the entire Nabanhe Reserve in the year 2023.In contrast to the BAU scenario, no further TPs occurred in the remaining years of the simulation for the Nabanhe Reserve.In both SW1 and SW2, positive RSIs indicate TPs in the year 2024, although less distinct than for the entire Nabanhe Reserve.No positive TPs are observed for SW3.Similar to the BAU scenario, negative RSI values in the year 2036 indicate TPs for SW1, SW2, and SW3.In SW1, the steep gain in forest areas in the initial years of BTO are a direct result of the establishment of water protection zones around water sources.In SW3, from 2018 onwards, the changes in the forest categories were only about ±1%.Changes in the rubber category were even less pronounced, ranging between 11.1% and 11.7% of areal coverage of SW3 throughout the simulation period.

Tipping Points in the Supply of ESS
The land use changes introduced in the BTO scenario generally improved regulating ESS and habitat quality.We observed the opposite in the BAU scenario over the 25-year simulation.These results confirmed the expectations developed in the scenario design process with local stakeholders [41].Similar results have also been reported in other ESS scenario simulations.Bai et al. [53] applied InVEST in Northern China to assess agriculture, forestry, and urban expansion scenarios and found the establishment of riparian buffer zones to be the optimal land management strategy, as it balanced agricultural production and hydrological ESS supply.In contrast to the distinct differences of regulating ESS and biodiversity between BAU and BTO, the time series of total rubber yields in the Nabanhe Reserve and the analyzed subwatersheds follow comparable trajectories for both scenarios but vary in the magnitude of rubber gains and losses.The time series of rubber yields turned out to be a decisive factor in dictating whether a TP has been reached or not.This can be seen in Figure 2a,c-f,h, as the models all agree on a negative TP for the year 2036, which coincides with the steepest declines in rubber yields.
As perennial systems, rubber plantations have an economic lifespan of about 30 years [54].The amount of latex which can be tapped during this period is highly variable.Productivity is low when tapping begins, peaks during the middle years of a plantation's lifespan, and declines again when the plantation approaches the end of its economic lifespan [40].In the study area, a broad range of plantation ages is present at the start of the simulation [51], which results in high variations of total rubber yield throughout the simulation period as existing plantations shift in and out of their economic lifespan in addition to the variations introduced by newly established plantations.The drop in total rubber yield leading to the TPs in 2034/2036 is a result of multiple plantations in the Nabanhe Reserve reaching the end of their economic lifespan simultaneously and having to be renewed.Rubber plantations are not the only source of income for the population in the Nabanhe Reserve [55].Paddy rice, maize, and fruit plantations are also relevant [55], but we assumed no changes of those categories in our scenarios.Since our analysis focused on deviations from initial conditions, we argue that the additional provisioning services provided by these land use categories can be neglected within the aims of this study.In addition, we argue that the strong dependency on rubber for the inhabitants of the Nabanhe Reserve justifies rubber yield as the defining provisioning service (like agricultural GDP in Hossain et al. [56]).
The analysis of subwatersheds, varying in their initial land use conditions and land use change trajectories, proved to be a crucial point for the interpretation of the results.For the entire Nabanhe Reserve in BTO, the land use changes led to a positive TP in 2023 (Figure 2b).This means that at a large scale, the increases in the regulating ESS and habitat quality are able to buffer against the decline in rubber yield in the later years of the simulation.We expected to find similar results for the three analyzed subwatersheds.These expectations have only partly been met, as we identified positive as well as negative TPs for smaller spatial extents (SW1, SW2, SW3) in BTO.For SW1, SW2, and SW3 in BTO, the results agree on negative TPs in 2036, comparable to the negative TPs identified in the BAU scenario.However, for SW1 and SW2, positive TPs indicate that the system shifted to a comparatively more beneficial state of ESS supply in 2024.This means that the negative TP found for SW1 and SW2 in 2036 in BTO represents a shift back to a comparable system state as that of the initial regime (2015-2023).In SW3 in BTO, this is not the case, as only a negative TP was identified in 2036.This indicates that the increases in other ESS are not sufficient to buffer against the steep drop in rubber yields in 2036 for SW3 in the BTO scenario.
In the BAU scenario, on the other hand, as linear growth of rubber areas and linear decrease of forest areas are set to occur with constant annual rates, the ESS z-scores in the Nabanhe Reserve and analyzed subwatersheds show very similar trajectories, regardless of the differing scales.The exception to this is SW2 in BAU.The high gain in rubber yields shifts the system to a comparatively more beneficial state of ESS supply in 2024.This is very similar to what was observed in SW1 and SW2 in the BTO scenario, the difference being that we observed increased habitat quality and regulatory ESS in BTO, as opposed to increased rubber yields alone in the BAU scenario.
The strengths of our method include the focus on smaller subwatersheds, since this revealed results, which would have gone unnoticed due to aggregation, if the focus had been on the Nabanhe Reserve as a whole.Furthermore, ESS research often focusses on comparing present land use conditions with altered land use "snapshots" of one specific target year in the future, lacking any analyses of the timeframe in between the initial condition and the final year of a scenario (see [30,45,57] as examples).As rubber is a perennial crop and the expected yield can vary substantially with plantation age and location, the annual resolution of our method is more robust in comparison to snapshot approaches.In addition, our method can be transferred to other regions of land use change situations, as a both the InVEST models and the R-method are freely available through online platforms [43,52].
As suggested by Rodionov [52], the applicability of the R-method should be tailored to the topic of study at hand.In Rodionov [52], the R-method was exemplified using the January PDO (Pacific decadal oscillation), which is well known to experience regime shifts.In our case study, no prior knowledge existed for potential future or past shifts in regimes of ESS supply.Analyzing past shifts is very data demanding, as exemplified in Zhang et al. [12].They analyzed time series (1900-2006) of 55 social, economic, and ecological indicators to assess the sustainability of ESS supply in Eastern China.In Zhang et al. [12], normalized trajectories of regulating services and provisioning services were grouped and depicted as separate time series.This revealed a clear trade-off between the two ESS categories, with declining regulating services and increasing provisioning services.In both scenarios analyzed in this study, the relationship of regulating services (and habitat quality) with rubber yields (the only provisioning service in our case study) is not as clear because of the variable amount of rubber yield throughout the simulation period.
Depending on the focus of research, definitions for TPs in socioecological systems and land use change can, but do not necessarily have to, include different concepts [10].Here, we did not focus on specific concepts to define TPs.Instead, we describe TPs as deviations from normal fluctuations in the presence of underlying trends [17].Nevertheless, we address some of the concepts (indicated by quotation marks in the following paragraph) commonly used in research on socioecological systems and land use change according to Milkoreit et al. [10].
In our results, "Multiple Stable States" can be seen as the regimes between TPs in time series of the ESS z-score.With the chosen models and their structure, we were not able to cover the concept of "Irreversibility".For some submodels we used (e.g., sediment retention and water yield), a reversal of land use change rates in BAU would simply result in an approximately linear increase for these ESS, as opposed to the approximately linear decrease shown in Figure 2a,c,e,g.The same can be said for the ESS trajectories of habitat quality and carbon storage.Improvements for both submodels could be made by including an analysis of forest patch connectivity or forest edge effects [58,59].The habitat quality submodel of InVEST does include edge effects at agricultural areas, as well as for roads and villages, but does not include any landscape connectivity indicators.We were able to include "Feedbacks" in a qualitative manner through the rules and outlines of the scenarios, which were codesigned by stakeholders.However, quantitative feedback is not included in any of modeling structures of InVEST or the rubber yield model.This means that changes in one ESS do not affect other ESS in the simulation.Therefore, the only dynamically changing input parameter throughout the simulation period is land use.The concept of "Abruptness" is crucial, as a majority of identified TPs are related to abruptly declining rubber yields in the case study areas.Hughes et al. [60] highlighted how slow regime shifts provide additional time for taking management actions in order to prevent a change in system states.Our analysis revealed habitat quality and the three regulating ESS to be variables that slowly react to changes in land use at the scale of the entire Nabanhe Reserve.Depicted as the ESS z-score, the decrease in habitat quality and the regulating ESS is masked by increases in rubber yield.This can be seen in the first regime in BAU until the years 2036 (Figure 2a,c,e) and 2034 (Figure 2g).With an abrupt decrease in rubber yield, the algorithm of the R-method revealed significant differences between the mean ESS z-score of the new regime (2036-2040) in comparison to the previous regime (2015-2035).

Methodological Limits and Future Research
The trajectories of habitat quality and most of the regulating ESS depicted in Figure 2 can be approximated by linear functions.These would qualify as type 1 of time-series categories described in Dearing et al. [17]-linear trends with environmental limits set by scientific, public, or political instances.We refrained from an analysis based on environmental limits, as no known thresholds for, e.g., habitat quality, exist for our study area.At the moment, without the quantification of environmental limits, our method can only provide the general direction in which the state of the socioecological system is moving, towards or away from the boundaries of the safe operating space.This, in turn, raises the question of if the system had already been in a state of unsustainable management before the start of the simulation, as all further analysis is related to the initial system state.The next step in improving our method would be to integrate environmental limits for the supply of the analyzed ESS.Our results revealed further weaknesses of using the R-method for the purpose of identifying TPs in modeled time series of ESS supply.For the BTO scenario in SW1, we expected the rapid gain of forest areas to result in a positive TP after the initial two years (Figure 2d).However, no TP was identified.This was due to the following reason: The R-method allows for the reduction of influence of outliers on the regime mean in the time series.The ESS z-score of the initial two years of BTO in SW1 are counted as outliers and, therefore, these values enter the calculation of the mean of the respective regime with reduced weights.Furthermore, the methodology of this study does not take any stakeholder preferences for specific ESS into account, as all ESS enter the calculation of the arithmetic mean (ESS z-score) with equal weights.With the inclusion of stakeholder preferences for ESS, the identified years of potential TPs may be subject to change [41].The potential impacts of stakeholder preferences and weighing methods have been discussed, e.g., in Thellmann et al. [41].Other points of improving the methodology would be to include natural disturbances (e.g., extreme weather events or pests) in addition to the human drivers we analyzed in this study.
In parts, our results are an example of the modifiable areal unit problem, meaning that results may change depending on where the boundaries are set [61].To circumvent this problem in future research, agent-based models of ESS provision might be used, which add explicit ESS "use regions" as a link between ESS "source" and "sink regions" [62].Recent advancements in simulating socioecological systems have shown promising results with the use of system dynamics models for the delineation of regional safe operating spaces [63].This approach of system dynamics is well suited to capture potential feedbacks for longer simulation periods but is also more dependent on past empirical studies in comparison to the approach we used here.In our case, a simulation period of 25 years was chosen as a realistic timeframe to predict potential future land use changes for rubber cultivation in the Nabanhe Reserve [41].Future research with the method proposed here could be conducted with longer simulation periods in combination with an earlier year for the initial conditions to: (1) improve validation of model runs, (2) capture long term effects of rubber plantation turnover times, and (3) reveal the impact of past rubber expansion on ESS in the Nabanhe Reserve.

Conclusions
In an era of global change, decisions may have to be made without full knowledge of potential consequences while making the best possible use of what is known at the time [64].Within this background, we developed a method of combining spatially explicit modeling with a data-driven, sequential algorithm, both of which are freely available, as an easy-to-adapt concept for land use planning in data-scarce environments.We used this method to identify potential TPs in ESS for rubber expansion scenarios in MMSEA, but the method can be adapted to other areas facing comparable land use change situations, such as deforestation driven by timber or palm oil production in other parts of Southeast Asia.In the study area, we discovered unexpected differences in the results of the TP analysis related to questions of spatial scale.The application of the same TP identification methodology for the same scenario (BTO) resulted in positive TPs for the entire Nabanhe Reserve (large scale), whereas for the analyzed subwatersheds (small scale), conditions remained unchanged or were comparatively worse at the end of the scenario when compared to the initial state.From this, we conclude that sophisticated land use planning is able to provide benefits in the supply of ESS at watershed scale, but that potential trade-offs at subwatershed scales should not be neglected.Even if land use plans aim at a more sustainable manner of production, specific local conditions may prevent them from being adapted.Developing management plans for socioecological systems on multiple spatial scales without exceeding the limits of regional safe operating spaces is a critical challenge that remains for future research.

Figure 1 .
Figure 1.Depiction of land cover in the Nabanhe Reserve in (a) the baseline year of 2015, (c) the final year (2040) of the Business-As-Usual scenario, and (e) the final year (2040) of the Balanced-Trade-Offs scenario.The 2015 land cover map features a resolution of 30 × 30 m and was derived from Rapid Eye satellite data.Three subwatersheds (SW1, SW2, SW3) are delineated as focal points for further analysis.Detailed excerpts of the three subwatersheds are depicted in (b) for the initial condition of 2015, (d) for the final year (2040) of the Business-As-Usual scenario, and (f) for the final year (2040) of the Balanced-Trade-Offs scenario.

Figure 1 .
Figure 1.Depiction of land cover in the Nabanhe Reserve in (a) the baseline year of 2015, (c) the final year (2040) of the Business-As-Usual scenario, and (e) the final year (2040) of the Balanced-Trade-Offs scenario.The 2015 land cover map features a resolution of 30 × 30 m and was derived from Rapid Eye satellite data.Three subwatersheds (SW1, SW2, SW3) are delineated as focal points for further analysis.Detailed excerpts of the three subwatersheds are depicted in (b) for the initial condition of 2015, (d) for the final year (2040) of the Business-As-Usual scenario, and (f) for the final year (2040) of the Balanced-Trade-Offs scenario.

Figure 2 .
Figure 2. Depiction of the modeled and normalized ecosystem services (ESS) indices for Habitat Quality, Rubber Yield, Sediment Retention, Water Yield, and Carbon Storage, as well as their annual arithmetic mean value (z-score) for the Business-As-Usual (BAU) scenario (left column; (a,c,e,g)) and the Balanced-Trade-Offs (BTO) scenario (right column; (b,d,f,h)) in relation to rubber-related land cover changes throughout the simulation period between 2015 and 2040.The rows of graphs relate to the whole Nabanhe Reserve (a,b), SW1 (c,d), SW2 (e,f), and SW3 (g,h).The land use changes are given in percent of the extent of the Nabanhe Reserve or respective subwatershed for each graph.Years with potential regime shifts are indicated by purple columns (RSI: regime shift index).Note the different scaling on the y-axis for the z-score of the different subwatersheds and the differences in scaling on the axes for rubber and forest coverage.The coloring of the letters (a-h) indicates if the regime of ESS supply at the end of the scenario is comparatively better (green), worse (red), or approximately the same (black) as the initial condition.

Figure 2 .
Figure 2. Depiction of the modeled and normalized ecosystem services (ESS) indices for Habitat Quality, Rubber Yield, Sediment Retention, Water Yield, and Carbon Storage, as well as their annual arithmetic mean value (z-score) for the Business-As-Usual (BAU) scenario (left column; (a,c,e,g)) and the Balanced-Trade-Offs (BTO) scenario (right column; (b,d,f,h)) in relation to rubber-related land cover changes throughout the simulation period between 2015 and 2040.The rows of graphs relate to the whole Nabanhe Reserve (a,b), SW1 (c,d), SW2 (e,f), and SW3 (g,h).The land use changes are given in percent of the extent of the Nabanhe Reserve or respective subwatershed for each graph.Years with potential regime shifts are indicated by purple columns (RSI: regime shift index).Note the different scaling on the y-axis for the z-score of the different subwatersheds and the differences in scaling on the axes for rubber and forest coverage.The coloring of the letters (a-h) indicates if the regime of ESS supply at the end of the scenario is comparatively better (green), worse (red), or approximately the same (black) as the initial condition.

Table 1 .
Proportions of land cover categories in the Nabanhe Reserve and the analyzed subwatersheds for the initial year of the simulation (2015).

Table 2 .
[41]al output of the Integrated Valuation of Ecosystem Services and Trade-Offs (InVEST) submodels and the rubber yield model for the initial state of 2015 and the final year of each land use scenario (2040) for the whole Nabanhe Reserve (NR)[41]and each subwatershed (SW1, SW2, and SW3).