Which Socio-economic Conditions Drive the Selection of Agroforestry at the Forest Frontier?

Models are essential to assess the socio-economic credentials of new agroforestry systems. In this study, we showcase robust optimisation as a tool to evaluate agroforestry’s potential to meet farmers’ multiple goals. Our modelling approach has three parts. First, we use a discrete land-use model to evaluate two agroforestry systems (alley cropping and silvopasture) and conventional land uses against five socio-economic objectives, focusing on the forest frontier in eastern Panama. Next, we couple the land-use model with robust optimisation, to determine the mix of land uses (farm portfolio) that minimises trade-offs between the five objectives. Here we consider uncertainty to simulate the land-use decisions of a risk-averse farmer. Finally, we assess how the type and amount of agroforestry included in the optimal land-use portfolio changes under different environmental, socio-economic and political scenarios, to explore the conditions that may make agroforestry more attractive for farmers. We identify silvopasture as a promising land use for meeting farmers’ goals, especially for farms with less productive soils. The additional labour demand compared to conventional pasture, however, may prove an important barrier to adoption for farms facing acute labour shortages. The selection of agroforestry responded strongly to changes in investment costs and timber prices, suggesting that cost-sharing arrangements and tax incentives could be effective strategies to enhance adoption. We found alley cropping to be less compatible with farmers’ risk aversion, but this agroforestry system may still be a desirable complement to the land-use portfolio, especially for farmers who are more profit-oriented and tolerant of risk.

The production costs of rice and maize are based on information from the Ministry of Agricultural Development of Panama Ministerio de Desarrollo Agropecuario de Panamá, MIDA: (MIDA 2019a(MIDA , 2019b. The annual cultivation costs remain constant over the 20-year period, except more labour is needed in the first year of cultivation to clear the site of secondary vegetation (Table S1).

Table S1
Production costs for rice and maize monocultures, based on information for traditional (non-mechanised) systems with use of pesticides and fertilisers (MIDA 2019a(MIDA , 2019b Maize is not cultivated every year in the alley cropping system, and therefore (unlike the monoculture system) high fertiliser inputs are not needed to sustain constant annual yields over the 20-year period. When growing maize as a monoculture MIDA (2019b) recommends applying 272 kg of both complex NPK fertiliser and urea per hectare. For maize cultivation under alley cropping, we apply this amount of NPK fertiliser but no urea. Furthermore, to conform with the requirements of the Forest Stewardship Council, no pesticides are applied in the alley cropping system.
We reduced the labour days and input costs needed for sowing, fertilising and harvesting maize in the alley cropping system by a factor of 0.83, to account for the reduced area on which these activities occur (compared to the monoculture system). However, we also included two additional labour days per hectare for harvesting maize in the agroforestry system, to reflect the extra time and care needed to protect trees from mechanical damage (Paul et al. 2017). When increased shading reduces expected maize yields by 50%, we apply half the fertiliser and reduce harvest, threshing and transport costs by 50%. Table S2 outlines the establishment and ongoing management costs of pasture, teak plantation and the two agroforestry systems. The timber systems incur an annual fixed cost of $60 per hectare to account for the higher complexity of these systems compared to pure agriculture, and therefore the need for more technical assistance and monitoring. In the cattle systems labour includes fence maintenance, monitoring cows, applying vaccinations, buying and selling cows (Paul et al. 2015).

Fig. S1
Projected standing timber volume (m 3 /ha) of each timber system over the 20-year rotation.

Table S3
Projected tree dimensions at age 20.

Crown radius (m) Comment
Teak 22.6 29.9 4.5 In line with growth measurements from a teak plantation in a comparable site in western Panama (Las Lajas, Chiriquí; Paul et al. 2015), as well as biophysical modelling for the study area using WaNuLCAS (Paul et al. 2017).
Cedar 24 28 2.6 Cedar growth rate within range reported for plantations in Latin America (Cintron 1990) and Costa Rica (Bellow and Nair 2003). Compact crown plausible given the prolonged pruning regime from years four to seven.

Table S4
Categories for yield reductions of maize in the alley cropping system based on canopy shading and tree heights -adopted from Paul et al. (2015).

Revenue from crops and cattle
Expected yields and producer prices for each agricultural crop and cattle are shown in Table S5. Rice and maize yields are compiled at the national level for traditional (non-mechanised) planting systems with use of pesticides and fertilisers. For maize grown in the alley cropping system, we reduce the initial expected hectare yield by a factor of 0.83, to account for the smaller area covered by the crop. In the cattle systems young cows are bought with an initial weight of 250 kg and fattened to 454 kg by the end of the 12-month ceba period. Based on knowledge of local farmers and agricultural experts we selected a stocking rate of 2 cows per hectare for conventional pasture. Beef prices are taken from the local cattle market. The initial stocking rate of the silvopastoral system is reduced by a factor of 0.95 to account for the reduced pasture area: the stocking rate then declines linearly with the canopy growth of cedar trees (see Figure S2 and Equation 1 in the main text).

Timber revenues
Following Griess and Knoke (2011) and Paul et al. (2017), only 64% of the standing wood volume was considered to be marketable timber: this accounts for 20% harvests losses and assumes that 80% of the harvested volume is stem wood. Prices for teak and cedar logs were obtained from the National Forest Office (ONF) in Costa Rica (ONF 2019) -see Table S6. Prices for the very small logs (<15 cm dbh) are taken from Paul et al. (2017). We applied a 0.5 reduction factor to the cedar price (ONF 2019), because timber quality is likely to be much lower in a silvopastoral system compared to plantation (Love et al. 2009).

Net cash flow
Using the expected costs and revenues for each land-use, l, we determined the net cash flow (NCF) for each land-use for each year, t, of the 20-year period (Table S8):   We account for yield and price variability through Monte Carlo simulations, in which expected prices and yields are adjusted based on historical data. For each simulation run, a, the bootstrapping process (random sampling with replacement) selects a year, h, from the historic dataset for each year, t, of the land-use model. The randomly chosen year is used to adjust the expected price (Equation S2) and expected yield (Equation S3) of each commodity, c, for that year of the land-use model to represent market and yield fluctuations respectively. The quotient of a given commodity price or yield for the randomly selected year, h, and the average price or yield for that commodity Reliable data on historic timber yields and prices were not available for Panama. Therefore, following Paul et al. (2017) and Castro et al. (2015) we assume a 10% coefficient of variation for timber yields. For prices we assume a 19% coefficient of variation for teak, and 8% coefficient of variation for cedar, based on the variation of timber prices in Costa Rica between 2007 and 2020 (ONF 2020). Within the land-use model, timber price and yields were adjusted accordingly for each simulation run.

Adding elements of uncertainty
Simulating variation in yield and prices drives variation in three of the socio-economic indicators: NPV, payback periods and food production. To simulate variation in labour demand and investment costs we randomly selected values from the normal distribution, assuming a 10% coefficient of variation in investment costs and yearly labour demand.
By repeating the Monte Carlo simulation 10,000 times we generated a frequency distribution for each indicator. From these 10,000 simulation runs we then calculated the mean (predicted) value G H ),! and standard deviation IJ ),! of each indicator for each land-use:

Formulation of multi-criteria optimisation model
The starting point for the multi-criteria optimisation model are the predicted values G H ),! , which form our estimate of the ability of each land-use, l, to achieve each indicator, i. However, recognising the uncertainty associated with these estimates, the model also considers potential deviations based on the standard deviation, IJ ),! . To do so we compute uncertainty adjusted values, G ),!,7 that span from a best-to worst-case estimate. For our best-case we take the mean score, while for the worst-case estimate we add or subtract a multiple, m, of the SD, depending on the direction of the indicator: for best case for worst case, if more is considered better for worst case, if less is considered better The model generates all possible combinations of best-and worst-case estimates across the considered land-uses in discrete uncertainty scenarios. This results in 2 L scenarios per indicator, where L is the number of land-uses considered in the optimisation. These uncertainty scenarios describe the surface (i.e. provide the corner points) of multi-dimensional boxes that represent our uncertainty spaces Ui for each indicator. The optimisation considers all corner points simultaneously. Therefore the optimal land allocation offers a feasible solution for all input values contained within the uncertainty spaces . In our study the optimisation considers 640 uncertainty scenarios (2 7 scenarios ´ 5 indicators).
Within each uncertainty scenario the model computes the performance of a hypothetical farm portfolio for achieving a given objective. The hypothetical farm portfolio comprises various shares, M ! , of each land-use. We compute the farm-level performance, " !,# , of this portfolio as the sum of the uncertainty adjusted values G ),!,7 (i.e. best and worst case estimates) within a given scenario, u, weighted by the area share, M ! , of each land-use in the hypothetical portfolio: The unit of this farm-level performance depends on the indicator (e.g. for NPV it is measured in $/ha, for food production Mcal/ha/yr). Therefore, to compare performance across different indicators, we normalise " !,# between 0 and 100%. To do so we set the best performing uncertainty adjusted value within each uncertainty scenario as our reference point, i.e. the 100% or target level. For "more is better" indicators the reference point is the highest uncertainty adjusted value within an uncertainty scenario, G ),7(809) * = max ! OG ),!,7(809) P, while for less is better indicators it is the lowest, G ),7(809) * = min ! OG ),!,7(809) P. To ensure robust results, reference points are computed across all uncertainty levels using m = 3, denoted as 0(R = 3). For each uncertainty scenario we divide the difference between the reference point and farm-level performance by the difference between the highest and lowest uncertainty-adjusted values (∆ ),7(809) ), to produce a normalised distance J ),7 to the 100% level: • 100 if more is better The variable J ),7 measures the shortfall between farm-level performance for a given indicator and the target (best possible) level, which can be interpreted as underperformance. We define V as the maximum J ),7 across all uncertainty scenarios: The variable V serves as our objective function, which we want to minimise. To solve the allocation problem we set area shares, M ! , allocated to each land-use as the decision variables. Equations S11 to S14 formulate the optimisation problem.
Minimise V (S11) subject to: By minimising V we aim to minimise the worst underperformance (highest J ),7 ) across all indicators. Inequation 12 is needed to linearise the objective function; the inequation summarises individual constraints (representing 2 L uncertainty scenarios per indicator) as V on the left side and the level of underperformance on the right side. This formulation of the problem can be solved exactly by the Simplex algorithm.

Prioritising individual objectives
In the scenario Prioritising individual objectives we test the effect of weighing individual objectives above the others. We derived weights, ^), by assuming an individual indicator is twice as important as the others (Table S11). For each uncertainty scenario, u, and indicator, i, we multiplied the distance between the achieved and target indicator level, J ),7 , by the weight, ^), of each indicator. The model then determined the land-use composition that minimised the maximum of these weighted distances.

Farmer preferences
In the Farmer preferences scenario we include the general land-use preferences obtained through farmer interviews , see Table S12) as an additional indicator in the multi-criteria optimisation model. These stated land-use preferences may serve as a proxy for more intangible or cultural values that are not easily captured by the other socio-economic indicators (Knoke et al. 2014;Temesgen and Wu 2018). Therefore, including them as an additional (equally weighted) indicator in the multi-criteria model can help to test how farmers' cultural values may influence the optimal land-use composition.

Fig. S3
Composition of the optimised farm portfolio (share of land area allocated to each land-use) for three levels of risk aversion: risk neutral (m = 0), moderately risk-averse (m = 1.5), and strongly risk-averse (m = 3.0), when farmers' land-use preferences (see Table S12) are included as an additional (equally weighted) indicator in the multi-criteria model.  changing the assumptions and coefficients of the land-use model. Input variables of the land-use model are progressively increased or decreased under three scenarios: changes to expected crop yields relate to the Lower crop yields scenario, changes in investment costs to Agroforestry subsidy and changes in teak and cedar price to Higher timber prices. These scenarios are described in Table 5 in the main text. Optimisation carried out from the perspective of a moderately risk-averse decisionmaker (m = 1.5).