Cascading Impacts of Payments for Ecosystem Services in Complex Human-Environment Systems

The theory and practice associated with payments for ecosystem services (PES) feature a variety of piecemeal studies related to impacts of socioeconomic, demographic, and environmental variables, lacking e orts in understanding their mutual relationships in a spatially and temporally explicit manner. In addition, PES literature is short of ecological metrics that document the consequences of PES other than land use and land cover and its change. Building on detailed survey data from Fanjingshan National Nature Reserve (FNNR), China, we developed and tested an agent-basedmodel to study the complex interactions among human livelihoods (migration and resource extraction in particular), PES, and the Guizhou golden monkey habitat occupancy over 20 years. We then performed simulation-based experiments testing social and ecological impacts of PES payments as well as human population pressures. The results show that with a steady increase in outmigration, the number of land parcels enrolled in one of China’s major PES programs tends to increase, reach a peak, and then slowly decline, showing a convex trend that converges to a stable number of enrolled parcels regardless of payment levels. Simulated monkey occupancy responds to changes in PES payment levels substantially in edge areas of FNNR. Our model is not only useful for FNNR, but also applicable as a platform to study and further understand human and ecological roles of PES in many other complex human-environment systems, shedding light into key elements, interactions, or relationships in the systems that PES researchers and practitioners should bear inmind. Our research contributes to establishing a scientific basis of PES science that incorporates features in complex systems, o eringmore realistic, spatially and temporally explicit insights related to PES policy or related interventions.


Introduction
. Many grand challenges -global warming, loss and fragmentation of forest areas, biodiversity loss, wildlife extinction, desertification, and the like -are besetting humanity at unprecedented rates. Yet virtually all these grand challenges can be traced back to rapidly growing human population and various human activities. Humans are degrading or destroying ecosystems at an alarming rate, jeopardizing their vital "life-support services of tremendous value" such as food, water, clean air, soil, and forests (Daily & Matson ; Millennium Ecosystem Assessment ). Even so-called protected areas are not exempted from such degradation worldwide (Curran et al. ; Liu et al. ). Facing such crisis, the International Convention of Biological Diversity's Aichi targets (https://www.cbd.int/sp/target/) have called for protecting natural habitats (Target ), Goal and research questions . The complex systems approach o ers a great potential to understand the cascading impacts of many ecological / environmental, socioeconomic, and demographic factors in a spatially and temporally explicit manner when their complex, reciprocal relationships are considered. The goal of this paper is therefore to measure and project the social-ecological impacts of PES over a relatively long term. We will use not only land use and land cover metrics as in previous studies, but also wildlife occupancy indicators.
. Specifically, we aim to answer three questions: ) what are the e ects of PES on human demography, livelihoods, and related activities? ) what are the spatially explicit e ects of PES on wildlife habitat use a er incorporating human demography, livelihoods, and related activity data? Moreover, ) what factors may make PES programs ine ective in the long run?

Study site
. Fanjingshan National Nature Reserve (FNNR), located in the northeast part of Guizhou Province, China, is a flagship reserve of subtropical ecosystems in China ( Figure ). FNNR is part of the global biodiversity hotspots (Myers et al. ), replete with over plant, animal, and bird species (GDF & FNNR ). Totaling km in area , the reserve (roughly . • N âĂŞ . • N and . • E -. • E) encompasses low elevation ( m) evergreen broadleaf ecosystems, mixed deciduous-broadleaf ecosystems at mid-elevations ( -m), and subalpine, meadow, and conifer ecosystems at higher elevations ( -m), manifesting large variability in biophysical conditions (Yang et al. ). The reserve is also home of the last and only population (around animals) of the Guizhou golden monkeys (Rhinopithecus brelichi; GGMs herea er), also known as Guizhou snub-nosed monkeys, an umbrella, engendered species highly sensitive to human presence, activity, and the resultant habitat degradation (Yang et al. ). .
FNNR is the home of over , villagers ( % are Tujia and Miao minorities), who live within or near FNNR boundaries in a subsistence lifestyle. These villagers are allowed to enter non-core habitat areas for resource collection and livestock herding, though illegal production of charcoal, wood collection, and poaching in the core habitat area also occur occasionally year round (Yost et al., in press). Over the last two decades, local villagers have increasingly migrated to cities for higher-pay jobs, and only returned to home villages for short times (e.g., during a spring festival). Migration in this paper is defined as outmigration outside of county boundaries, and the dominating migration destinations are coastal big cities. The total number of out-migrants was ( . %) out of all the households we surveyed in (Sections . -. ). In when GTGP started, there were very few out-migrants ( for all migrants). Once GTGP was in operation, a substantial (near exponential) increase in the number of migrants occurred, suggesting that the related laborers, once released due to GTGP, may have migrated out at a slower pace at earlier times yet at an accelerated speed later. We therefore consider GTGP to be a trigger and catalyst of migration because it not only freed labor that had been bonded to farmland, but also provided a start-up fund (e.g., for travel expenses) for those relatively poor households to migrate out. .
This study focuses on the northern area of FNNR, where two villages (Taohuayuan and Pingsuo) closest to Yangaoping (approximately . • N, . • E) are located. Yangaoping is the breeding area where GGMs gather to mate in April and September (the red circle in Panel D, Figure ). Human settlements typically are located around the borders of the reserve, which are at lower elevations.

Social survey .
We administered household interviews in . Using the roster of all , households used in the FNNR census as our sampling frame, we surveyed households (yellow squares in Panel A, Figure ) based on a stratified random sampling strategy (detail in Yost et al., in press). The survey focused on ) individual-level characteristics: age, gender, education level, etc. of each member in the surveyed household; ) household characteristics: living conditions, household economy including incomes, expenses, and time, number of people, and income related to migration and local o -farm jobs; ) household land use and PES characteristics: , households surveyed in , and sites with cameras mounted in ˜ . In Panel A, the area above the pink, horizontal line represents the northern area we simulated, the black rectangle the area shown in Panel D, and red circle Yangaoping. total land area, amount of time (labor), and income or compensation related to GTGP; and ) household extraction of local resources: time, frequency, amount, and location (recorded in Google Earth maps) related to various resources (Yost et al. ). .

In
, we revisited all these households with similar questions, although focus was more on household land use, participation in GTGP, and detailed household income. We ended up with a survey of households in . In we surveyed all households in five natural villages within two administrative villages of Pingsuo and Taohuayuan (blue triangles in Panel A, Figure ). The aim was to obtain full coverage of all households in these two administrative villages for agent-based modeling, verification, and data analysis purposes. We surveyed all the available households in these two administrative villages, ending up with a total of households. All the surveys were conducted under permits from the San Diego State University's Institutional Review Board (Protocol #: ).

Ecological survey .
To investigate vegetation structure and habitat use of wildlife under human disturbance, we established sampling plots in FNNR. Each plot was m x m (a subset of them is shown as pentagons in Panel D, Figure ). Location of plots was decided based on accessibility, elevation, distance to other plots, and suggestions provided by FNNR sta and local field guides with the goal to spread out plots across FNNR. For each plot, we recorded species of understory, midstory, and overstory vegetation and estimated the percentage of cover for each species. We also deployed a Bushnell Trophy Cam infrared camera at each plot to monitor presence of mammals (> . kg) and pheasants from April to January , relocating cameras to other places to increase site diversity (Chen et al. ). Typically, cameras were mounted on trees . to m above the ground. We set cameras at auto-sensitivity to record three photos upon each detection, with a -second delay between .
The ABM creates an x grid spanning the entire FNNR, with a spatial resolution of approximately m. The spatial resolution is chosen to balance the need to represent seasonal movement of monkey family-groups and human resource extraction activities and simulation speed. Among the cells within the lattice, are within the FNNR boundary with the following attributes assigned to them: elevation and land cover or use types: Bamboo, Coniferous, Broadleaf, Mixed Forest, Lichen, Deciduous, Shrublands, Clouds, Farmland, Household, Farm, PES Forest, and Outside_FNNR. These vegetation categories were generated through a combination of supervised and unsupervised image classification routines from near-anniversary-date summer Landsat satellite imagery (Tsai et al. ). Elevation data were extracted from a m ASTER DEM GeoTIFF released by NASA in (Tsai et al. ).

.
Our representation of the environment and ecology is from two aspects. First, we incorporate land use and land cover data as mentioned above. Second, we use golden monkey habitat occupancy as an indicator of intactness of a certain land area for GGM habitat use (the inverse of this measure is the degree of human intrusion): Occupancy = # of simulated GGM visits on a pixel per year ( )

.
Our choice of this occupancy measure is based on the ecology of GGMs: They avoid direct encounter with humans (Yang et al. ). This metric for the impact of PES programs enables us to address the need to quantify broader, ecological dimensions of PES as mentioned in Section PES and Challenges.
With GGM life traits well-documented (e.g., Bleisch et al. ; Wu et al. ; Yang et al. ; Zhang et al. ) an agent-based model was developed and posted online, which simulates the GGM agents' seasonal movement and their life events such as birth, death, and grouping (Mak ). The literature shows that GGMs live and travel in family groups of -, while larger groups (Bleisch et al. ; Yang et al. ) and all-male groups are observable. Yet in the model larger groups of more than individuals were split into two. For space limit we do not elaborate other GGM features in the model (Mak ). .
Building on (Mak ), the ABM in this paper models GGM habitat use or occupancy primarily with consideration of elevation, vegetation types, and avoidance of human presence in the context of human livelihoods and demographic changes. This modeling choice also hinges upon the fact that FNNR has high primary forest cover, representing excellent habitat that provides food and water for golden monkeys (detail in Yang et al. , p. ). This situation makes it reasonable to not incorporate GGM feeding or water seeking behavior in our model.

.
GGMs prefer to stay in areas from to m (Wu et al. ), and in extreme situations to m (Bleisch et al. ). Accordingly, the model limits GGM family-groups to elevations between m and m, human settlements (i.e., they escape when approaching a radius = m centering each settlement), farms (radius = m), and GTGP enrolled land parcels (radius = m). Finally, whenever away from human settlements, GGM family-group agents either move according to free semi-random walk based on their vegetative surroundings or are traveling to and from Yangaoping in mating and breeding seasons. The monkey agents tend to stay in mixed, deciduous, conifer-broadleaf, and broadleaf forests, but also travel through other types of vegetation (Bleisch et al. ; Wu et al. ; Yang et al. ; Zhang et al. ).
Human demographics and livelihoods .
The ABM generates land parcel objects (one special kind of agents) that belong to a total of household agents. Then the model assigns household attributes (e.g., household ID, village, x and y coordinates, dryland area, paddy land area, GTGP land area, non-GTGP land area) to each of these parcels. These parcels may convert between GTGP and non-GTGP land uses, depending on the PES program implemented within the model as well as a set of socio-demographic variables (Section GTGP participation and land use). A sum of human agents are also created with individual attributes such as person ID, age, gender, education level, marriage status, and working status (work on own farm, o -farm agricultural work, non-farm business, being student, not working, unable to work) assigned to each agent based on the survey data. Each person agent goes through relevant life events: being born, growth, going to college, marriage (simultaneously bringing in a female or male from outside the two villages ), bearing a child, death, migration (outmigration and re-migration listed below separately), and collecting resources. For space limit, we give detail of bearing children (for the rest we refer interesting readers to the Python code at https://www.comses.net/codebases/ 0852f2b8-6517-4b83-b7fa-8304eb538421/releases/1.0.0/): .
Once a woman reaches (the minimum age for a woman to marry), assign her a birth plan (or fertility rate) representing the number of children for her whole life: she may have , , , , , or children in her whole life with the probabilities of .

.
When ˜ years (a parameter that is randomly chosen to be , , or years) elapse since her marriage (i.e., she reaches years old), she may have her st child. Then a er another ˜ years (another parameter that is randomly chosen to be , , or years), she may have her nd child. . . till her birth plan is accomplished (An et al. ).

Resource extraction .
Based on our survey data of the households, each household sends a person agent to designated cells to collect resources such as fuelwood, medicinal herbs, bamboo shoots, mushroom, fodders (for pigs and oxen), and others while the order of collecting these resources is randomly chosen. At each time step, the person agent only gathers one resource at a time. The person agent is usually the household head (if within the range ofyears old); if s/he migrates out, the next capable person is designated as a substitute head of household and collect these resources.

.
When human agents venture out to collect resources at a certain time step, their paths are marked to act as a barrier at the same step for GGM family-group agents that travel around the reserve freely but avoid humans. Along with human avoidance and seasonal migration to breeding habitat (the circle in Panel A, Figure ), Guizhou snub-nosed monkeys also avoid low elevations and choose weighted semi-random paths based on favored vegetation.
. Each person agent gathers this resource at the corresponding frequency determined by the survey. Each time a resource collection action takes o en one day, part of the -day time step that includes the day. During the step, all the cells that the person agent traverses (from the household cell to the collection cell) receive a weight of / for avoidance. This weight is proportional to the collection frequency. The mismatch in temporal resolution may a ect the influence human presence exerts on GGM habitat use, which will be discussed in Section Discussion. Human resource gathering and other activities in the northern area (Panel A, Figure ) may be suppressed with influences such as higher income or higher rates of human out-migration. Below we describe how we represent outmigration.

Outmigration .
Modeling human decision has been a big challenge in understanding and envisioning human-environment systems, and ABM may o er unique power in this regard (An ). Migration happens only if the household under consideration has a person between and . According to our field survey in (Yang ), the annual probability of outmigration can be calculated based on a logistic function: where annual local o -farm income (income_local_o _farm) and remittances from migrants (mig_remittances) are measured in Yuan, number of laborers (num_labor) is the number of people between and , migration network (migration_network) is defined as a binary variable about whether the household has social connections to earlier migrants ( for no and for yes), the area of non-GTGP farmland per laborer (non_gtgp_land_per_labor) is the average amount of farmland divided by the number of people between and , GTGP participation (gtgp_part) is also binary ( for no, for yes), marriage is binary and refers to whether someone is married or not ( if no, if yes), and farm work status (farm_work) is a binary variable ( for no and for yes)). The rest of the variables are self-evident.
. However, the migrants can choose to return to their rural home, and we model the corresponding probability as a function of the age at the time of the most recent migration (age) and the number of years s/he has migrated out (yr_mig) according to a logistic regression by Yang ( ): Therefore, a person might migrate out, return, re-migrate out, and continue this cycle if s/he is between and years old and the randomly generated numbers are less than the probabilities determined by Equations and are greater than the threshold. The coe icients in Equations and are results of the above-mentioned logistic regression based on survey data collected in (Yang ).

GTGP participation and land use
. The ABM reads in parcel level data for each household from a CSV file. The CSV file contains whether a parcel is enrolled in GTGP (named a GTGP parcel) or not (named a non-GTGP parcel), x and y coordinates, whether the land is dryland or paddy land (land_type), and the amount of time to travel to the corresponding land parcel from the household (time_land).
. Local farmers may fallow some parcels even without any payment for various reasons. Therefore we combine fallowing of land with GTGP participation, resulting in an action of returning farmland to other uses, including planting cash crops, fallowing land, or planting ecological trees (Yost et al. ). A logistic function is used to calculate the probability of returning farmland to other uses (including participating in GTGP and leaving the parcel fallow) at parcel level. Yet for terminology consistency and conciseness, we still call it GTGP participation and the corresponding probability is named GTGP_par_prob: This equation was derived through logistic regression analysis based on data from (Yang ), where land_type and time_land are explained earlier, and GTGP_net_cash is the di erence between GTGP payment and income from growing crops on the same land parcel.
. When enrolling non-GTGP parcels, there are three conditions to meet simultaneously. First, the area of non-GTGP land in the household must be greater than or equal to . mu (a parameter based on the data from ; subject to sensitivity analysis) to provide vegetables. Second, a randomly generated number is less than the probability calculated using Equation . Third, the GTGP payment is greater than zero. Note that if GTGP payment is zero (GTGP ends), there is still a small chance that the household may fallow the parcel (as mentioned earlier we still call this action participating in GTGP for terminology conciseness) at a much lower probability. This low probability defaults at GTGP_par_prob * . , where parameter . is based on data from and subject to sensitivity analysis (same for . below). .
When deciding to reconvert a GTGP parcel to non-GTGP parcel (i.e., farmland), a very small probability (GTGP_par_prob * . ) is calculated, where . is a parameter subject to sensitivity analysis (Sections . -. ). Once a parcel is enrolled in GTGP, it will be reconverted to farmland once a randomly generated number is less than this small probability.

Model evaluation .
The model evaluation is comprised of two elements: model verification that warrants that the code does what it intends to do and is free of coding bugs, and model validation that compares model structure, processes, and results with expected ones based on data, theory, or experts' opinions (An et al. ; Manson ). As model verification is progressive (debugging continues), the results below contribute to both verification and validation. .
We first show the dynamics of model projected human population size, births, deaths, and number of migrants over years (on a monthly basis) ( × steps) for the two villages in the northern area ( Figure ). As we do not have empirical data for population size a er , we are not able to compare our projected population sizes with empirical data. Instead, we calculate various measures such as population increase rate, births, and deaths, and then compare them with data from other sources. The overall population dynamics shows a slow increase (i.e., from in to around in ), which can be explained by the accumulated numbers of births, deaths, and migrants in Table : There are births and deaths, implying natural increase of persons. There are marriages, bringing in people from outside. So the total number of increase in population size is + = ; if subtracting the increase in number of out-migrants ( -= ), the total increase through migration is -= persons. So the total population in should be (base) + = , which explains the simulated population size . in Year (or Step ) in Table . These results exactly corroborate the projected population size in Month (Year ) in Figure (under BP = . ). When birth plan (fertility rate) is set to be . and . , the total population in changes to and respectively (see more in Sections . -. and Section ). .
We also validate these demographic results using external data and information. From the simulation results in Table , the annual birth rate is ( -) / ( × ) = . %, or . ‰. According to China's National Bureau of Statistics ( ), the nationwide birth rate is . % or . ‰. Considering the decreasing trend in the most recent years (e.g., from . ‰ in to . ‰ in ), our average annual rate of . ‰ for the next JASSS, ( ) , http://jasss.soc.surrey.ac.uk/ / / .html Doi: . /jasss. years is reasonable. Similarly we calculate the annual death rate to be ( -) / ( × ) = . % or . ‰, while the same national rate is . % in (National Bureau of Statistics ).
Step ) and Year ( ), yet we can take a close look at specific age groups for verification purposes. For instance, the male, -group is % of the total population (le ), which turns to be around % ten years later (Middle) and % years later (right, Figure ). We explain these decreases using the migration, growth (age-out), and death data: The male group at age -had people in Year , while died and migrated out at Years ˜ . Given that the simulated total population sizes at Years ˜ are , we calculate the percentage of this group at Year to be ( --) / , or / = . %, which matches up with the Year bar for males aged -in Figure . Similarly, we explain other percentages in Figure   .
Third, we examine the projected total area of farmland, which decomposes to GTGP and non-GTGP land area for each of the HHs (Table ). For verification purposes, we examine land use areas at time , time , and time . The average number and area of non-GTGP parcels both decline with increasing payment, while at the same time the same numbers for GTGP parcels increase. This finding is consistent with literature regarding the impact of payment on PES enrollment (Wunder ; Yost et al. ).
Fourth, we simulate GGM births, deaths, and overall population size from Month # to Month # ( Figure ). The overall GGM population dynamics (Figure a) is consistent with literature (Yang et al. ). Yet the overall population dynamics alone may not warrant correctness of processes that constitute such dynamics. For instance, a mistakenly high birth date may o set the impact of a high, yet wrong mortality rate (at some age group) and thus generate "correct" population size, yet giving rise to unbalanced, likely incorrect age structure. Below we examine this potential issue and explore the projected GGM age structure. .
We run the ABM and generate GGM age structure (i.e., percent of individuals of -, -, -, -, -age groups in total GGM population) at Year and ( Figure ). It seems that the two age-structures at Years and are similar to that at Year that is based on literature, yet there are some small di erences (e.g., lower % for age -group and higher % for age -group). For detailed information and discussion we refer to Mak ( ). Figure : Simulated GGM age-structure at years , and . .
Fi h, we evaluate habitat occupancy using several methods by: ) presenting the ABM-generated occupancy maps to FNNR managers, seeking their qualitative evaluation of these maps; ) visually comparing the ABMgenerated occupancy maps with a paper habitat map (n o digital file) generated from FNNR sta 's long-term fieldwork (Yang et al. ); and ) by comparing occupancy outcomes generated by the ABM (Figure ) with a small subset of camera-based GGM density data, i.e., only data from the top camera sites (where we have recorded human activity from the households). For each cell, we standardize the GGM occurrence number to a gradient between and by dividing the observed captures by the max captures we have observed there. Then we standardize the ABM predicted occupancy numbers using the same procedure. We calculate the Pearson's correlation between them. Out of the above three methods, the last one ( ) is quantitative, yet subject to several limitations (details in Section ). .
The spatial evaluation of habitat occupancy ( Figure ) from experts' opinions is very positive, suggesting that the overall trend in GGM habitat use has been captured (Y. Yang, personal communication). The qualitative comparison between the FNNR paper map and our ABM-generated map also suggests a good match. The quantitative test gives a moderate result, in which the Pearson's correlation coe icient is . ( . when a pair of data suspicious of equipment error are removed). We ascribe this to the small sample size (n = ), short term of data collection ( year) compared to our simulation length ( years), and equipment bias (e.g., no-detection of GGMs). However, worthy of mention is that the ABM-based occupancy results have been compared to the results from an independent, di erent modeling approach called MaxEnt, and the results largely match (Mak ). . Last, we perform sensitivity analysis using the parameter-sweeping approach (An et al. ). In Sections . -. , the model performance may be sensitive to three parameters, i.e., three radii (centered on a human settlement, farm, and GTGP land parcel) that each define a zone of avoidance by GGMs. We therefore design a GGM-escape test (Table ). As each parameter has values for test, there are a total of × × = experiments. For each experiment, we ) generate and report habitat occupancy data at the end of simulation and ) calculate the kappa index between each experiment map and the baseline map at the end of simulation. .
Similarly, there is uncertainty regarding minimum amount of cropland that a household hopes to keep and probabilities of enrolling cropland in GTGP or reconverting GTGP land back to cropland (Section GTGP participation and land use). We therefore design a second test named land-decision test (Table ). As each parameter has values for test, there are also a total of × × = experiments under this land-decision test. For each experiment, we collect data regarding ) the number of GTGP parcels and number of non-GTGP parcels over years ( months); ) occupancy data at the end time, and ) kappa value between experiment occupancy map and baseline occupancy map at the end time.  The above sensitivity analysis results are positive. The two tests do not give any model crash, nor lead to uninterpretable outcomes. The Kappa values in the GGM-escape test have a largely decreasing trend when the threshold (i.e., threshold count for # of steps a cell is used by GGMs: above this value the cell is counted as occupied) increases. The spatial locations of occupied habitat also agree with our experience and experts' opinion (Appendix B). In land-decision test, the numbers of GTGP-parcels and non-GTGP parcels still converge even when di erent parameter combinations are used in the model, suggesting that the outcome in Figure is robust, not likely being an ad hoc outcome (Appendix B).

Experiment design .
Once we have evaluated the model and deemed to represent processes in a satisfactory manner, the model is then used to perform scenario experiments. Here we are interested in seeking insights into the impact of GTGP policy on long-term land use, migration, and GGM habitat use. We consider three scenarios: ) Scenario , where GTGP policy stops (payment = ); ) Scenario , in which GTGP stays as is (payment = Yuan/Mu); and ) Scenario , in which the payment is doubled (payment = Yuan/Mu). .
We also examine the potential impact of varying population pressures. In this regard, we design two population scenarios with GTGP payment set at default ( Yuan/Mu): Scenario features a lower birth plan or fertility rate ( . children per woman), in which an eligible woman may have , , , , , or children in her whole life with a probability of . , . , . , . , . , and . . In contrast, the probabilities are .
, . , and . in the baseline simulation (corresponding to . children per woman; see Section Human demographics and livelihoods). Scenario features a higher birth plan of . , in which an eligible woman may have , , , , , or children in her whole life equally with a probability of . , . , . , . , . , and . .

.
For each of the above scenarios, we intend to estimate the amount of land enrolled in GTGP land (accumulated GTGP land area of HHs), number of accumulated migrants (note that returnee migrants are subtracted), and the related densities of GGMs at FNNR over time. To quantitatively measure the degree of change in GGM occupancy between scenarios, we adopt the Cohen's kappa statistic that o ers an overall agreement (no change in our application) and disagreement (change in our application) on the cell-by-cell basis (Carletta ).

Experiment Results
. ) and increases to around under Scenario (BP = . ). The three birth plans do not a ect the number of migrants (Figure ), which may arise from the fact that our time span ( years) is not long enough for most new births to grow to migration age. .
Next, we turn to the impacts of GTGP policy and population pressure on land use (Figure ). Payment levels do make a big di erence in increasing GTGP enrollment in earlier times, e.g., prior to Month in Figure a. This finding is consistent with relevant literature, where the positive impact of payment levels on participation has been reported extensively (Wunder ; Yost et al. ). Yet at later times, some interesting patterns emerge: The non-GTGP land decreases to a nadir, then rises slowly (a concave curve), while the GTGP land rises till a peak is reached, followed by slow decrease (a convex curve). The di erence made by the PES payments (or no payment) is the level of increase or decrease and the timing of reaching the peak and nadir. Interestingly, the number of parcels in all three scenarios turns to converge to around two parcels near the end of the -years' time span (we will discuss this finding later). When examining the impact of population pressure, the above convergence pattern is still prevalent at all birth plan values. Yet birth plan values seem to have no impact on land use for the first half of simulation time, and small di erences begin to occur at later times: the increase in non-GTGP parcels and decrease in GTGP parcels tend to become slower, resulting in a later convergence time (Figure b). .
Worthy of mention is that the above convergence pattern occurs in other parameter settings. As shown in the sensitivity analysis section, this convergence trend also occurs in the experiments with very di erent parameter values (Appendix B) and is not likely an ad hoc outcome due to this specific parameter setting.
Figure : Impacts of GTGP payment level (a) and impacts of birth plan (b) on enrollment (i.e., # of GTGP and non-GTGP parcels). .
The changes in spatial patterns of habitat use densities are quite informative. When interpreting these maps, it is preferable to reference to Figure , which portrays that the surrounding areas are less used by GGMs. In this context, we can see increasing payment level does improve habitat quality in edge areas, especially areas surrounding the two villages. This arises from local villagers' resource extraction patterns; mostly they go to nearby areas but do visit remote places occasionally. Therefore, there is no substantial improvement in remote (o en core) habitat areas due to a higher GTGP payment. We discuss why this comes out in Section Discussion.
Higher or lower fertility rates (birth plans) do not generate big di erences in GGM occupancy visually (Scenarios and vs. Baseline; Figure ), yet quantitative analysis of Kappa index does show that a decrease in birth plan (from . to . ) makes bigger changes in GGM occupancy more than an increase (from . to . ) does a er taking into account the randomness in generating the occupancy measures (Table , Appendix C). .
For the northern part of FNNR (our study site) with households, the number of migrants shows a slight, nearlinear increasing trend from around to around individuals from to ( Figure ; also see Table  ), representing an annual rate of . person/household. This trend is lower than the general trend between and for all the households we surveyed in . There were almost no migrants in and in , representing an annual rate of increase of . person per household. This higher rate during and is due to the rapid increase near ˜ due to the implementation of GTGP in (Yang ). On the other hand, the payment levels did not change the total number of migrants over time substantially except that stopping payment gives rise to slightly lower number of migrants from Months to ( Figure ). . When evaluating the di erences in occupancy due to changes in payment, it appears that compared to the status quo scenario (Scenario ), the payment level makes a moderate di erence in occupancy when threshold is relatively high. When using a threshold of visits/cell (i.e., visits/year; we chose this relatively high value to assure that changes are not by chance, but by substantial changes in GGM visit frequency a er consulting experts on GGM ecology and behavior) as a threshold to classify a cell as occupied (a binary classification is required by the Kappa calculation), the Kappa statistic between and payment levels (i.e., Scenarios vs. ) is .
, while that between and payment levels (i.e., Scenarios vs. ) is . . This suggests that changes made by increasing payment are larger than changes (degradation) made by canceling the payment.

Discussion
. As reported above, the average number of parcels from model simulations based on all three GTGP scenarios and two population pressure scenarios converges near the end of the -year time span (Figure ). Specifically, GTGP land parcels increase until a peak is reached and then decrease; correspondingly non-GTGP land parcels decrease until a nadir is reached and then increaseâĂŤthis convergence pattern occurs regardless of payment level (Figure a). We explain this outcome first by food security. Local households tend to keep a small portion of land unenrolled (i.e., non-GTGP land does not decrease to zero) despite great benefits associated with enrolling in GTGP. During our interviews, respondents repeatedly mentioned this choice, and we translate this choice into a lower bound of non-GTGP land (parameter min cropland in Table , default at . mu). This may reflect the great attraction of labor from local businesses and particularly migration (Zhu ).
. Interestingly, the convergence pattern occurs regardless of birth plans (fertility rates) as well, and birth plans have little impact on GTGP enrollment, especially at earlier times of simulation ( Figure b). As household size has a small, positive impact on GTGP enrollment (see Equation ), an increase in household size due to a larger birth plan (fertility rate) may stimulate the household to enroll more farmland in GTGP, yet other factors such as the above lower bound of non-GTGP land and distance between farmland parcel and household (Equation ) may stop or dilute this increase in enrollment. For instance, if the increase in household size happens to be in a household with available farmland quite far away, then the tendency in probability increase (due to bigger household size) may be o set by a decrease due to higher household-land distance. Under this condition, capturing spatial heterogeneity (as we did here in terms of mapping all households and their land parcels and thus being able to calculate household-land distances) is very important. Therefore, this surprising convergence pattern is jointly accounted for by the complexity existing in the system (Liu et al. ), likely associated with feedback loops and interactions among system components. Although with initial conditions specified by users, any ABM so ware program could decompose to a set of discrete-time or discrete-event di erence equations in principle, such analytical representations become very a dauting task when the system becomes highly complex (Tesfatsion ). In this regard, ABM has unique power when addressing questions in complex systems as shown above. .
Why does higher GTGP payment increase habitat in edge areas? This finding may arise from the following paths. When GTGP pay increases, the GTGP net cash increases (Equation ), which will escalate the probability of GTGP participation and simultaneously decrease non-GTGP land. According to Equation , a decrease in non-GTGP land will promote out-migration -once people of -migrate out, local households tend to limit visits to surrounding areas for resource extraction such that GGMs can visit these places at a higher frequency. Again, this finding, though insightful and understandable, is also not one that is intuitive or was expected beforehand. For instance, another potential outcome -GTGP leads to habitat degradation -is possible if this pathway exists: GTGP participation leads to considerable labor released from farming, then local farmers may spend "extra" time extracting local resources for self-use and sale at local markets. However, our data and model results do not support this outcome. .
In summary, the merits of this paper lie in the following aspects. First, our research moves forward PES research, which is in nature embedded in both human and environment systems. Participants of PES programs are humans, who make decisions and act in response to various human-environment contexts including PES regulations. Following a set of empirically-established behavioral rules, GGM agents rover on the landscape and respond to changes in the environment, especially human activities. All such processes take place on a spatially explicit landscape, leading to changes in the complex human-environment system. In this context, we have employed a complex systems framework, characterized by an agent-based model that integrates data of varying types, processes or relationships operating at multiple temporal and spatial scales, and knowledge from both natural and social science disciplines. .
Our research is innovative in using wildlife (GGM) occupancy as a measure of environmental changes (Scullion et al. ), which may be applicable to evaluation of other environmental policies. In addition, movement analysis and modeling is a hot topic in space-time analysis and animal behavioral ecology (An et al. ; Tang & Bennett ). Our semi-random walk approach -in the context of human disturbance -may provide better understanding of animal behavior. Furthermore, our agent-based modeling practice contributes to ABM testing (e.g., using structural measures such as age structure, parameter sweeping test), model transparency and reusability (all code in the CoMSES repository), and modeling of human decision -a daunting challenge that attracts increasing recognition in many simulation and modeling (including ABM) domains. Such contributions may help advance complexity science (Manson ), addressing the YAAWN syndrome in the ABM arena (O'Sullivan et al. ) and many other challenges in the ABM domain (An et al. ). .
Caveats arise concerning spatial and temporal resolutions. In our model, the temporal resolution is days, while resource collection takes place on a daily basis. This mismatch ( -day vs. -day time step) has forced us to spread the -day influence into a -day period (see Section Resource extraction). This may bias the simulated influence of human visits on GGM habitat use, especially the temporal dimension of such influence. In the future, a finer temporal resolution (e.g., daily or even hourly) should be tested if simulation speed is not sacrificed much, or parallel or cloud computing techniques can be employed. A similar issue is for the relatively coarse spatial resolution (near m). In the future, a m spatial resolution may be adopted, which enables us to map GGM and human activities on a more precise basis.
. Also worth of mention is the moderate quality of spatial projections about GGM occupancy. Although able to capture generic trends of GGM habitat use, the model results are not in high agreement with our observed camera-trap data. In addition to the reasons mentioned above (low sample size, short term of data collection, and equipment bias), we also o er one important reason: The spatial locations of cameras, subject to human accessibility, are not very representative of all GGM accessible habitat areas. At such human-accessible locations, collection of human activity data may be more challenging due to reasons such as local people's unwillingness to share their activities or spatial inaccuracy when mapping human activities.

Conclusion
. This paper aims to reveal complex, reciprocal relationships among PES, human livelihoods, and the environment (both land use and land cover and habitat occupancy) at Fanjingshan National Nature Reserve, China.
The agent-based model contributes to integrating data from di erent spatial and temporal scales and disciplines, revealing land use and habitat patterns that are di icult to obtain otherwise. Instead of being used as a predictive tool, we recommend that our model be used as a platform to study and further understand complex human-environment systems, shedding light on key elements, interactions, or relationships in such systems. E orts in this regard will help us establish PES science that incorporates features in complex systems, such as heterogeneity, feedback, and nonlinearity, o ering more realistic, spatially and temporally explicit scenarios related to human policy or intervention. All these contribute to achieving the United Nations' Sustainable Development Goals (United Nations ), especially Goal that aims to protect, restore and promote sustainable use of terrestrial ecosystems.

Appendix A: ODD protocol
The ODD (Overview, Design concepts, & Details) Protocol for an agent-based model is a standardized document which outlines a model's purpose, variables, framework, schedule, and data. The format was conceptualized by a team of twenty-eight authors who had previously published or worked with agent-based models, and serves as a universal set of guidelines for describing a model (Grimm et al. , ).

Purpose
This agent-based model serves a variety of inter-connected purposes: • To simulate the demographic changes of humans living in the FNNR, including births, deaths, marriages, out-migration and re-migration. These changes are modeled on the existing human data gathered from the FNNR. Heads of household are also designated. Various other statistics, such as income level and education level, are projected as well. Finally, human age structures are recorded.
• To simulate GGM movement within the FNNR in a movement sub-model, which follows seasonal patterns of migration to a mating area and avoids human settlements or low elevations. Movement is also weighted according to the nearby vegetation: monkeys are more likely to move to vegetation that they are modeled to favor.
• To simulate human resource collection in the movement sub-model, which may impact the movement of humans upon GGM habitat.
• To simulate Green-to-Grain Program (GTGP) enrollment or dis-enrollment and GTGP land conversion, which is based on factors such as current income and types of land owned.
• To simulate the demographic changes of the Guizhou Golden Monkey (referred to henceforth as GGM) of the Fanjingshan National Nature Reserve (referred to henceforth as FNNR) over time, including births, deaths, formation of new families from a large group or of all-male groups, and mating behavior. Monkey demographic (age and gender) structures are also recorded.

Entities, State Variables, & Scales
Time -Each time-step of the model represents approximately days. Therefore, every steps of the model, the model "advances" one year. Individual processes such as birth, death, and adulthood are continuous, and may occur at any time-step once conditions are met. While not an entity in itself, the passage of time will trigger events, such as birth or the formation of new groups. Spatial: One pixel represents a non-moving household, each of which has zero or one human agents who is a fuelwood collector and a varying number of total human agents within it. Any human agents that travel will always return to their household pixel.
Temporal: Like monkeys, humans age, sometimes reproduce, and die. Their lifespan and other variables are set in the model.

Environment (grid cell)
Potential Subtype: Household, Farm, PES, (Managed) Forest Elevation Vegetation Type Spatial: The extent of the FNNR. Elevation is adapted from m raster DEMs, scaled down by to create a manageable grid.
The sq. area of the rectangular grid containing the FNNR extent -that is, not the FNNR itself -is approximately sq. km. One environmental grid cell represents approximately m in diameter, or .
sq. km in area.
Temporal: None, though the agents interact with the environment by avoiding areas of low or high elevations. Temporal: Every step, there is a small chance GTGP conversion will be checked. If it is checked, then there is a chance that a land parcel may convert to a GTGP or non-GTGP land parcel, which will affect household-level variables such as income or land type. This chance is based on a regression formula that includes factors such as land type, current household income, and unit compensation.

Process, Overview, & Scheduling
All agents move in a random order, which does not a ect the processes within each step. The processes are carried out in the same order for every agent; for example, for individual monkey agents, the death check always occurs at the end of every step. Possible agent actions, which may be restricted by certain agent attributes such as gender or age, are as follows: Model (Step only): Create GGM family agents, monkey agents belonging to each family, human agents, land agents that refer to household-level lists that humans access, resource agents, and environmental grids Land Parcels (each step, land submodel only): simulate GTGP conversion, update household income, update GTGP area changes from conversion Humans (each step, movement/visualization submodel only): Head to resources, gather resources, head back to house, and randomly select another resource to gather GGM Family (each step, movement submodel only): Avoid humans if paths cross, avoid areas of low or high elevation (seasonally), move to Yangaoping (breeding site), move away from Yangaoping, and move to neighboring cells according to a correlated random-walk (path determined by vegetation and elevation) as described by Ahearn et al. ( ).
GGM Individual (each step, population submodel only): Age, possibly give birth (if female and of age), possibly out-migrate to another group (depending on age, gender, and current group sizeâĂŤmale monkeys may defect to an all-male group, and large families can split into smaller groups), and die Human Individual (each step, population submodel only): Age, possibly give birth (if female and of age), possibly marry, possibly out-migrate (sometimes in a special case due to college, where they typically do not re-migrate from) and re-migrate, gain education levels, change work status, possibly change head-ofhousehold status, and die

Design concepts
Discussed here are the eleven design concepts of the ODD protocol: basic principles, emergence, adaptation, objectives, learning, prediction, sensing, interaction, stochasticity, collectives, and observation (Grimm & Railsback ).

Basic principles
Visualization submodel âĂŞ-This assumes spatial patterns among all GGM families, such as yearly migration to Yangaoping; movement patterns are also calibrated according to movement needs as determined by known travel speeds, vegetation preferences, and behavior around humans from the literature (Grimm & Railsback ; Yang et al. ). Because the model input is directly based on observations from a field study, the default output -a showcase of movement over ten years -is expected to be fairly predictable. What is new about this model is the comparison of its output -a point density map of its movements -to a Maxent model through Cohen's Kappa, and a discussion of how di erent versions of this model -based on configuration settings such as GTGP unit compensation -may di er through a similar comparison.
Population submodel -This confirms GGM population structure observations from Yang et al.'s field study ( ) by modeling changes to the population over ten years based on birth, mortality, and birth-interval rates per age category. It does not consider intermediate factors that are not currently well-understood by the literature, such as low genetic diversity and the impact of this phenomenon on birth defects or miscarriages, or whether or not closely-related monkeys can breed. However, it considers the observed patterns of male departure from groups and refreshed fertility a er a recent loss of an infant. The population submodel also assumes a stable human age demographic, stable migration, and consistent birth, death, and marriage rates in relation to China's national averages.
GTGP conversion submodel âĂŞ-This assumes that households will enroll in GTGP given that they meet a certain threshold as determined by the regression formula, and that they are likely to revert a er a number of years without compensation a er the program ends.

Emergence
Visualization submodel âĂŞ-Avoidance of humans and human settlements may have more or less of an impact than expected; vegetation weights may more or less than an impact than expected.
Population submodel -âĂŞ Fertility "recovery" a er loss of an infant may have more or less of an impact than expected.
GTGP conversion submodel -âĂŞ GTGP unit compensation may have more or less of an impact than expected; the model also does not account for any land area changes as a result of GTGP enrollment in the visualization submodel.

Adaptation
Visualization submodel -âĂŞ Family agents avoid less-desirable cells with human settlements and lower elevations, and choose their neighboring cell by weights determined by vegetation type.
Population submodel âĂŞ-Fertility "recovery" or the lessening of between-birth intervals, a er the loss of an infant is an adaptive trait that female individuals may have.
GTGP conversion submodel âĂŞ-GTGP enrollment a ects income and GTGP land area, which plays a role in future GTGP enrollment.

Objectives
Visualization submodel -GGM Family agents primarily avoid humans, move to ideal vegetation, and visit Yangaoping for breeding and giving birth. The human agents' only objective is to gather resources and return to their homes.
Population submodel -Adaptations -and therefore objectives -are more chance-based rather than individualobjective based. For example, human and GGM agents have a small chance of dying every step, which cumulates to a yearly mortality rate every steps.
GTGP conversion submodel -Adaptations -and therefore objectives -are more chance-based rather than individual-objective based. For example, land parcels have a chance of conversion and adjusting household income every step.

Learning
This model does not have agents change their traits responsively as part of its process, but can be modified to do so using "self" variables. The closest feature it has is using the outputs of processes, such as increased income from GTGP conversion, to inform other future decisions.

Prediction
Agents do not estimate future consequences of decisions; decisions made are based on information available at the time, and o en impact current decisions immediately therea er, or on a cumulative basis.

Sensing
Visualization submodel -GGM Family agents sense the presence of humans gathering resources and human settlements, and will move so as not to overlap them. They also "sense" the vegetation and elevation around them.
Population submodel -Female monkey agents give birth at the correct ages, and if male, may migrate to other groups.
GTGP conversion submodel -Income and land parcel interaction a ect each other, so higher temporary income may result in lower GTGP enrollment, which creates a negative feedback loop (lower GTGP enrollment may also temporarily lower income).

Interaction
Visualization submodel -Humans gather resources, but otherwise do not interact with GGMs. The model may be changed later, e.g. to implement a poaching behavior.
GGMs avoid humans, and will not move to occupy the same or sometimes even an adjacent pixel (each pixel is a -meter space) as a human agent.
Population submodel -if any monkey infant dies, its mother may give birth again the following year (a "recovery" of fertility), even though the normal birth interval rate is years. Humans will gather resources less e iciently if they are marked as a resource-gathering household, but the head of the household/workers of age are currently migrated or die without a replacement.
GTGP conversion submodel -increased household income from GTGP conversion may result in higher outmigration rates; the end of the PES program will cause some households to revert their GTGP-enrolled land parcels.

Stochasticity
Visualization submodel -GGM agents move to their -cell (diagonal and adjacent) neighbors according to a weighted choice, which in turn is informed by the vegetation type of each neighbor. Males can also travel between groups or form an all-male group when of a su icient age.
Population submodel -For both monkey and human agents, births and deaths are random chances determined by a birth-interval rate or yearly rate.
GTGP conversion submodel -A weighted-probability-based formula decides whether or not a land parcel converts to GTGP, which in turn a ects human household income.

Collectives
In both the movement and population submodels, GGM individuals belong to a family; occasionally, all-male families can break o .
Land parcels belong to a household collective, which humans also belong to; the collective is accessed via lists, so there are no household agents.

Observation
The output of the movement model is a .csv file of all points that individual agents traveled to, from which a heatmap or point density map may be generated.
The output of the population submodel are .csv files tracking changes in the population and age/sex structure over time.
The output of the GTGP conversion submodel is a .csv file of non-GTGP and GTGP land parcel counts, areas, and household income tracking.

Initialization
All initial values for humans and GGMs are either taken from Yang et al.âĂŹs ( ) field study of the GGMs (in the case of population structure rates), or from data which was gathered as part of the greater FNNR project (in the case of environmental, resource, or human household data).
Visualization submodel -Initial settings (number of years the model runs for, number of monkey families, GTGP compensation structure, PES program span, etc.) are set by the user in fnnr_config_file.py before running the model.
Population submodel -Each family group contains -agents, and in total, there are between -monkeys in the reserve (very likely -).
GTGP conversion submodel -Income is determined at the start by current land income and o -farm income combined. Some land parcels are already enrolled in GTGP.

Submodels
The processes for the three submodels are outlined here.
Step âĂŞ all agent and environment types (families, individual monkeys, humans, resources, land parcels, environment) are parametrized and created.

Input Data Data Source
Guizhou Golden Monkey Seasonal movement behavior (qualitative), birth rates, death rates, birth-interval rates, groupmigration or infant-loss behavior (qualitative) Field study (Yang et  . Avoid humans and human settlement bu ers (cannot enter certain occupied cells or face lower random odds of entering weighted cell).
. Head to or from Yangaoping if it is directly before, during, or a er mating season or birthing season (e.g. September, or steps -in a year of steps).
. Avoid certain low or high elevations. If traveling to or from Yangaoping, they may temporarily pass through those cells.
. Move to neighboring cells (usually -times in a step to match the distances covered over five days, or one step, as noted in the literature) based on a random choice a ected by weights assigned to each neighbor, which in turn is determined by one of nine vegetation types: mixed, broadleaf, deciduous, conifer, bamboo, shrublands, lichen, clouds (usually random; artifact from classification process), or farmland.

Visualization submodel ordered priorities for human agents:
. If at home, choose a random resource from an imported list of resources their household gathers, and head towards the resource in a shortest-distance path.
. Once at the resource, head back home to deposit the resource, and repeat the process.
Since humans gather resources faster than the time resolution of the model (one time-step represents days), the visualization will show human agents "jumping" back and forth between resources up to each step; however, the coordinates traveled along the paths are recorded.

Population submodel for GGM individuals:
. Face a low-level mortality rate each five-day step (slightly less than the yearly mortality rate divided by , since * = days in one year, because of compounding probability). Mortality rates di er by age category.
. If female and of age, birth interval increases every step; if it exceeds (with no recent infant loss), give birth. If it exceeds (recent infant loss), also give birth. Once a female has given birth, their birth interval resets to , and builds up again over time.
. If male and of age, and if enough are "flagged" by a low-level chance once they reach of age, an all-male group may break o , and males will change families.
. Each step, population dynamics are recorded (in some versions of the model; other versions only record the first and last step), and at the final step (ten years = steps; twenty years = steps), a .csv file is generated. This file includes the starting/ending population and average age/sex structure of the population.

Population submodel for Human individuals:
. Face a low-level mortality rate each five-day step (slightly less than the yearly mortality rate divided by , since * = days in one year, because of compounding probability).
. Every step, there is a low chance a new baby will be born. If that chance is met, a random married female who has not given birth in the last two years will be selected to bear a child, and the household size will increase. Once a female has given birth, their birth interval resets to , and builds up again over time.
. If single and of age, there is a low chance of marriage occurring every step. Divorce is not accounted for in the model.
. Miscellaneous variables such as education and work status update semi-randomly depending on other weighted factors such as the person's age.
. Humans may migrate out or re-migrate back if they have migrated. These are based on regression formulas that consider gender, age, income, migration networks, ratio of land owned to laborers in the household, and education level.
. Migration changes who is designated as the head of the household, which may in turn impact resource gathering (alternative heads of households, especially if not of age, may gather resources more slowly).
. Each step, population dynamics are recorded (in some versions of the model; other versions only record the first and last step), and at the final step (ten years = steps; twenty years = steps), a .csv file is generated. This file includes the starting/ending population and average age/sex structure of the population.

GTGP conversion submodel:
. Each step, a small probability for GTGP enrollment or dis-enrollment is evaluated for each land parcel. This is based o a formula that considers the head of the land parcelâĂŹs householdâĂŹs gender, age, education, and income. It also considers the time taken to travel to the parcel from the household, the land type, and the number of humans in the household (household size).
. Non-GTGP land area, GTGP land area, and household income (for all land parcels of that household) are updated to reflect changes from GTGP enrollment.

Results of GGM-escape test
Below is the GGM habitat occupancy under the GGM-escape test, in which three radii (centered on a human settlement, farm, and GTGP land parcel) are exposed to sensitivity test. Note that each radius defines a zone of avoidance by GGMs. For each scenario or parameter setting, the simulation is run (or depending on the scenario) times with results collected and saved in a certain directory. Once all simulation results are collected, we calculate values of several key variables such as human population size, # of migrants, number of GTGP or non-GTGP parcels, etc.  Table : Kappa results for the GGM-escape test. Notes: * Here for each of the non-baseline simulations (call them S , S , . . . , S ), there are runs. Run of S (the st non-baseline simulation) is compared with Run of the baseline setting ( . , . , and . are for Min cropland, Prob multiplier , and Prob multiplier), and the Kappa is calculated for Run ; we replicate this for Run , Run , . . . up to Run . Finally, an average Kappa is calculated for S : KappaS . The numbers reported in each cell below are the average of Kappa values calculated this way, e.g., .
= (KappaS + KappaS + . . . + KappaS ) / . * * Here Run of baseline simulation ( . -. -. for the three parameters) is compared with Run of the same baseline simulation and Kappa is calculated. We repeat this for Run , . . . , Run . Finally, the average Kappa is calculated at each threshold level. Visually comparing maps for spatial locations of habitat occupancy (Figure ), we do not see big di erences due to the three parameters. Then we turn to examine the Kappa values. It appears that the model is not very sensitive to changes in the three parameters settlement radius, farm radius, and land radius ( Table ) given the magnitudes of Kappa index among baseline results (the right column in Table ). When the threshold increases from to , the di erences in Kappa between a certain GGM-escape test and the baseline tend to increase as well.

Results of land-decision test
From the dynamics of the number of GTGP parcels and that of non-GTGP parcels, we can see that the converging trend still exists (Figure ). Note that we sweep the three parameters min cropland (Mu), probability multiplier , and probability multiplier at values specified in Table of the main text. Below we show the di erences in Kappa due to changes in the three parameters of min cropland (mu), probability multiplier , and probability multiplier (Table ). Comparisons between two Kappa values at the same threshold show that changes in the three parameters give small-to-moderate changes in habitat occupancy.   Table : Kappa di erences as a function of birth plan. Notes: * BP = Birth plan, which is total fertility rate (average # of children born by each human female who reaches childbearing age); see Table for calculation of Kappa; * * Because there exist stochastic processes in the model, we run another runs under the same parameter setting (BP = . ) to show Kappa di erences between scenarios and between threshold counts. See Table for calculation of Kappa.

Notes
According to an online survey, the numbers of articles and authors reporting the development or use of ABMs have risen steadily in an exponential rate since the mid-s, covering a wide range of disciplines (An et al. ) In FNNR has been approved to a World Natural Heritage Site, and its territory will likely be expanded (Shi, personal communication).
Such groups consist of males aged -that have migrated out from original groups.
The lower bound could be m (Bleisch et al. ) m in extreme situations (Niu et al. ).
This and other two radii are based on personal communications with expert in GGM behavior, Mr. Yeqin Yang (February -March , ). They are parameters subject to changes. See Sections . -. for our sensitivity test following An et al. ( ).
There is a small chance that two people within the two villages may marry, thus not bringing in a person from outside. As the probability is very small according to An et al. ( ), we overlook this to simplify the model.
If a return migrant comes back, then s/he is not considered as migrant. I want the total number of people who are outside of FNNR and considered as migrants at each month (aggregate to -time steps).