Climate risks to Brazilian coffee production

Coffee is highly valuable agricultural export for the global South and accounted for US$846 billion of global agricultural export value in 2017 (FAO 2018). Although it is produced in a small area globally (11 million hectares) (FAO 2013), coffee production contributes substantially to foreign exchange in producing countries (Akiyama 2001). In the last two decades, Arabica coffee farmers in Latin America have endured numerous threats to their livelihoods due to plummeting coffee prices, a shift in the power of the supply chain structure away from producers toward consumer-facing companies (i.e. roasters and retailers), and changing rural demographics (i.e. out-migration and an aging population (Bacon 2005)). Coffee farmers in many regions now face the risk of declining coffee yields and quality due to global changes in temperature and precipitation (Gay et al 2006, Schroth et al 2009, Zullo et al 2011). Given the already low profit margins (and low farm-gate prices) associated with coffee production (Bacon 2005), declines in coffee yields could make existing coffee producing regions both economically and biophysically unsuitable for production in the future (Zullo et al 2006, de Sousa et al 2019).

The most suitable growing regions for Arabica coffee are where the annual temperature average is between 18 °C and 22 °C (Camargo 1985b) and in mountain regions where the altitude is above 1000 or 1200 m (Schroth et al 2015). Both very high and very low temperatures can result in yield and quality losses (Zullo et al 2011c, Camargo 1985b). Climate change is expected to alter temperatures and precipitation, and increase extreme events in many coffee producing regions throughout the world (Läderach et al 2010, Ovalle-Rivera et al 2015). These impacts will undoubtedly compound the economic and demographic challenges currently faced by small-scale coffee farmers.

Throughout the world, there is increasing evidence that coffee farmers are already suffering from the impacts of climate change (Baker and Haggar 2007, Haggar 2016). In their study of Mesoamerican coffee producing regions, Läderach and co-authors concluded that by 2050, changes in temperature and rainfall will decrease the area suitable for coffee production by at least 40% (Läderach et al 2010). In Brazil, the world's largest Arabica and Robusta coffee producing region, climate change is predicted to substantially reduce the amount of suitable area for production (Zullo et al 2006). In the Brazilian state of São Paulo, for example, the proportion of climatically low-risk areas for coffee production may decrease by ∼20% by 2050 under a 1 °C increase in temperature relative to the current climate and as much as 75% under a 5.8 °C increase (Coltri et al 2012). A majority of the 308 000 coffee farms in the country are operated by 'family' (or smallholder) farmers,⁶ who are potentially less capable of adapting to changes in temperatures and precipitation (IBGE 2017a). In the absence of adaptation, these shifts in climatic suitability are likely to substantially reduce coffee farmers' incomes and potentially force them to abandon production.

In this study, we assess the relative climate risk of coffee producing regions in Brazil through the concepts of hazard (potential harms to humans, infrastructure, and ecosystem services), exposure (the presence of humans, infrastructure, and ecosystem services in places that experience hazards), and vulnerability (the propensity of a system to be negatively affected by hazard) as defined by the IPCC (IPCC 2014). The intersection of climate hazards, exposure, and vulnerability influences the level of risk that a community faces. Where risk is high, it is very likely that losses in well-being will occur if steps are not taken to reduce risk, either by reducing exposure to the hazard or by increasing communities' ability to adapt to the changes induced by the hazard.

The impacts of climate hazard on crop yields, particularly annual crops, are becoming increasingly well understood (Knox et al 2012). Yet, the impacts of climate hazard on farmers' livelihoods are more difficult to assess due to a lack of longitudinal data about changes in rural livelihoods. Despite these data limitations, snapshots of the material measures of farmer livelihoods and broader socio-economic conditions do exist and can be used as a proxy to assess farmers' potential adaptive capacity in each region. Once potential adaptive capacity is assessed, it is possible to understand the degree to which high climate hazards and exposure are likely to overlap with high vulnerability. In regions where there is high overlap between hazards, exposure, and vulnerability, additional institutional support for climate adaptation will be most urgently needed.

To date, few studies have combined analysis of hazards, exposure, and risk at broad spatial scales. Nor have any studies in Brazil focused on identifying the relative vulnerability of different coffee communities and how it intersects with climate hazards and exposure. To contribute to this knowledge gap, our paper has two aims (i) to assess spatial and temporal variations in climate hazards, exposure, and vulnerability and (ii) to identify which regions have the highest overall climate risk (overlap between high hazards, exposure, and vulnerability).

To identify regions where coffee communities are most at risk to climate change, our methodological approach combines: (i) econometric modelling of the historical relationship between climate and coffee yield (climate hazard), (ii) mapping of coffee prevalence (climate exposure), and (iii) newly developed municipal-level estimates of relative climate vulnerability, based on underlying socio-economic characteristics. To draw inferences for the future, our analysis focuses on changes in the overlap of hazards, exposure, and vulnerability over time (2006 and 2017) and space (coffee producing municipalities in Brazil) using all available historical data.

Many studies have analyzed climate hazards through the lens of potential impacts on crop yields. These studies draw on panel statistical models of historical weather data to explain the relationship between climate change and crop production at the global, national, and regional level, particularly for individual annual crops (Lobell and Field 2007, Schlenker and Roberts 2009, Schlenker and Lobell 2010, Butler and Huybers 2015). To date, the focus of most climate and agriculture studies has been annual crops, such as maize, wheat, and rice, which are critical to global food security (Knox et al 2012). There are fewer studies of perennial species, which occupy less area globally, yet represent an integral component of many rural livelihoods (Samberg et al 2016, Hong et al 2020).

The existing literature on coffee-related climate change impacts and risks focuses heavily on ecological niche modelling, often coupled with machine learning techniques, to identify transitions in the suitability of regions for coffee production (Davis et al 2012, Ovalle-Rivera et al 2015, Schroth et al 2015, Bunn et al 2015, Pham et al 2019). Relatedly, several studies determine changes in climate risk by assessing to what degree future climate projections fall out of the optimal temperature or water deficit range, and exceed the frost probability threshold (Zullo et al 2006, 2011c, da Silva Tavares et al 2018). Statistical assessments of the impacts of historical changing temperatures and precipitation on coffee yields are limited and thus far constrained to India, Mexico, and Tanzania (Gay et al 2006, Craparo et al 2015, Jayakumar et al 2016). These studies in India, Mexico, and Tanzania all found significant negative yield impacts from increasing temperatures and decreasing rainfall (Gay et al 2006, Craparo et al 2015, Jayakumar et al 2016). To date no paper (to the best of our knowledge) has empirically examined the relationship between historical changes in coffee yield and historical climate change in Brazil. Instead, existing research on the impacts of climate change on coffee production in Brazil has focused on simulating potential changes in yields (Verhage et al 2017) and climate zoning (Zullo et al 2006, 2011, da Silva Tavares et al 2018).

The concepts of vulnerability and risk encompass a wide variety of potential components, many of which are difficult to measure, and few of which have been analyzed for causal impacts. Given these challenges, one leading empirical approach to both vulnerability and risk research has been to develop quantitative indices that vary over space and time (Cutter et al 2003, Cutter and Finch 2008, Hahn et al 2009, Cutter et al 2010, Gbetibouo et al 2010, Pandey and Jha 2012, Shah et al 2013, Ahsan and Warner 2014). These indices are based on existing theory and available data and provide estimates based on a range of variables that are theoretically likely to influence societies' ability to prepare for or adjust and respond to stress and mediate risk. By quantifying the conditions and assessing the variations in a single metric, indices are particularly useful for comparing the relative levels of social vulnerability over time and space (rather than absolute vulnerability). For work that spans larger geographical areas, rather than individual communities, vulnerability and risk research has highlighted the importance of the attributes of place, and thus many indices focus on the aggregate social aspects of a geographical location (Cutter 1996, Cutter et al 2003).

Given that overall risk encompasses both physical stress and human's capability to adapt, vulnerability indices tend to include both physical or socioeconomic data (Harlan et al 2006, Johnson and Wilson 2009). Whereas specifically 'social' vulnerability indices (e.g. Cutter et al 2003, Cutter and Finch 2008) encompass only socioeconomic and geographical characteristics. Since climate hazard and exposure are measured seperately, our vulnerability index focuses on the socio-economic dimensions of vulnerability. The interplay of climate hazard, exposure, and vulnerability then determine overall climate risk—the likelihood that climate hazards will negatively impact human well-being.

3.1. Study region

Our analysis centers on the South, Southeastern, and Center-west coffee producing states of Brazil (Goiás, Mato Grosso do Sul, Paraná, Minas Gerais, São Paulo, Rio de Janeiro, Bahia, and Distrito Federal), which produce approximately 90% of the country's Arabica coffee (IBGE 2017b) (figure 1). This region is expected to encounter substantial changes in climate; mean annual temperature is expected to increase by 4 °C in the summer and 2 °C to 5 °C in the winter by 2100 in the RCP 8.5 emissions scenario (Pachauri et al 2014). Higher temperatures are expected to reduce coffee bean quality and generate more favorable conditions for pests and diseases. Under all scenarios, the ideal climatic conditions for coffee production are expected to shift to the south of Brazil (Zullo et al 2011). Our unit of analysis is all of the municipalities in this region that produce coffee and have data for all variables of interest (n = 935). This is the smallest spatial scale at which socioeconomic data are available across the study region.

Zoom In Zoom Out Reset image size

3.2. Assessing climate risk

Our methodological approach for determining overall risk follows the conceptual approach of the IPCC (IPCC 2014)—a score for overall risk defined as the sum of individual estimates of climate hazards, exposure, and vulnerability:

$$\begin{equation}{\text{Risk}} = {\text{Hazard}} + {\text{Exposure}} + {\text{Vulnerability}}.\end{equation} \tag{ 1 }$$

To generate this overall risk score, the estimates of hazard, exposure, and vulnerability, explained below, are converted into indices along a uniform scale (1–5) using Jenks natural breaks and then summed (resulting in a total value ranging from 0 to 15). We also present results using breaks defined by Equal Weights Intervals and Quantiles.

3.2.1. Climate hazard

To measure climate hazard, we examine how historical changes in precipitation and temperature have influenced coffee yields. We used global gridded (0.25 degree) monthly average air temperature and total precipitation data from 1974 to 2017 (the period for which we also have data on municipality-level coffee yields) by Willmott and Matsuura.⁷ This dataset interpolates weather station data and is drawn from recent versions of the Global Historical Climatology Network (GHCN version 2) and the Global Surface Summary of Day archive. The version of the temperature and precipitation data are 1900–2017 Gridded Monthly Time Series V 5.01. From the dataset, the gridded data are averaged across each coffee growing municipality.

We average the gridded climate data at the municipality level for each coffee phenological season. Given the differential exposure of coffee to weather fluctuations across different phenological stages (Camargo 2010), we separate out the blooming (September to November), ripening (December to May) and harvesting periods (June to August) for coffee plants in the study region. The blooming period is the time where coffee flower bud blooms, initiated by the first rains after the dry season (Barros et al 1999). During the blooming period, high temperatures coupled with a lack of rainfall can impact coffee flower buds (Camargo 1985a). High temperatures (above 23 °C) are also thought to be detrimental during the ripening period (when the coffee beans develop and mature) and harvesting period, by influencing both yield and quality (Camargo 1985a). Finally, excess rain during the harvest period can inhibit optimal harvesting.

The historical relationship between coffee yield, temperature, and precipitation is examined using the following panel econometric model:

$$\begin{align} Log{(Yield)_{it}} & = Tmp\_bloomin{g_{it}} + Tmp\_ripenin{g_{it}} \nonumber\\ & \quad + Tmp\_harvestin{g_{it}} \nonumber\\ & \quad+ Tmp\_blooming_{it}^2 \; \nonumber\\ & \quad + Tmp\_ripening_{it}^2 \nonumber\\ & \quad + Tmp\_harvesting_{it}^2 \nonumber\\ & \quad+ Log{(Pcp\_blooming)_{it}} \; \nonumber\\ & \quad + Log{(Pcp\_ripening)_{it}} \nonumber\\ & \quad+ Log{(Pcp\_harvesting)_{it}} \nonumber\\ & \quad + Log{(Pcp\_bloomin{g^2})_{it}} \; \nonumber\\ & \quad + Log{(Pcp\_ripenin{g^2})_{it}}\nonumber\\ & \quad + Log{(Pcp\_harvestin{g^2})_{it}} \nonumber\\ & \quad + municipalit{y_i} + yea{r_t} + state\_tren{d_j} \end{align} \tag{ 2 }$$

where Log(Yield) is the log of coffee yield (kg per hectare) in municipality i and year t, which spans from 1974 to 2017. Tmp_ is the average monthly temperature and Log(Pcp_) is the log of the average monthly precipitation during a particular season. Both the yield and precipitation variables are highly non-normally distributed. This non-normal distribution leads to heteroscedasticity in the relationship between precipitation and yields. To avoid this problem, we log-transform both variables, including the squared precipitation terms. The temperature variables are normally distributed so they are not log-transformed. We included municipality and year fixed effects, as well as a state (j) time trend, to control for unobserved, non-climatic factors influencing yields.

Predicting yield on climate factors may be too complex to describe using only linear relationships (Watson 1963). Using the Ramsey RESET test, which tests whether the model is missing important non-linearities, we confirmed that the results supported the use of quadratic terms for average temperature and precipitation to specify a quadratic functional form. In addition, a quadratic functional form generates optimum values by calculating the highest point in the curve and enables better modeling of crop yield responses to climate variables (Gay et al 2006).

We conducted a Hausman test to test whether fixed effects or random effects are more suitable for the panel data analysis. The p-value is significant at the 99.9% level, supporting the use of fixed effects. We include municipality fixed effects to account for unobservable spatial variation and time-invariant effects such as soil type and elevation (Welch et al 2010, Blanc and Schlenker 2017). This makes the model less prone to omitted variable bias because the unobserved factors that influence yield are controlled. Additionally, the model estimators of the fixed effects enable consistent estimation of the effect of independent variables. Likewise, a year fixed effect is added to control for time-varying shocks, such as macroeconomic factors. After conducting a joint test to find whether the dummies for all years are equal to 0, the F test was significant at 99.9% confidence level, suggesting that year fixed effects are needed. Lastly, a state- trend is used to control the broad regional productivity improvements that may be occurring due to technological change or changes in economic policies. The state-specific technological trends control the long-term trends in yields unrelated to weather variability.

The climate hazard of each municipality is then defined by combining our coffee yield model with historical trends in climate, following existing approaches to calculate the impacts of historical climate trends on yield (Lobell et al 2011, Butler et al 2018). First, for each municipality, we calculated the temporal trends in temperature and precipitation within each phenological stage using a linear regression. We then identify the net yield impacts of these climate trends by using Eq. 2 to model yields at the beginning and end of our 44-year period, given state-year yield trends. The difference in predicted yields due to climate is defined as the climate hazard, with more negative yield impacts indicative of a greater climate hazard. A limitation of this approach is that it relies on a deterministic model of climate-yield impacts, which may not be as suitable as a flexible time series approach for capturing climate-yield relationships, given the heterogeneity in coffee yield trends and variability across Brazilian states (figure S1 (available online at stacks.iop.org/ERL/15/104015/mmedia)) (Agnolucci and De Lipsis 2020).

3.2.2. Climate exposure

Our definition of climate exposure focuses on the presence of rural livelihoods in places that are likely to incur climate hazards. In case of rural coffee producing communities in Brazil, their exposure to climate hazards is proxied by the prevalence of coffee production in that region. As a whole, municipalities with a large coffee cultivation area as a proportion of their total area are likely to be more exposed to changes in climate that impact coffee yields than municipalities where coffee production occupies only a small proportion of the agricultural area. To map exposure, we generated a variable called coffee prevalence, defined by the percentage of agricultural area that is in coffee production using agricultural census data from the Brazilian Institute for Geography and Statistics for both 2006 and 2017.

3.2.3. Climate vulnerability

Here we develop a Vulnerability Index (VI) to characterize spatial variability in factors likely to influence communities' propensity to be negatively affected by climate change, following several prior studies (Cutter et al 2003a, 2010, Cutter and Finch 2008, Hahn et al 2009, Gbetibouo et al 2010, Pandey and Jha 2012, Shah et al 2013, Ahsan and Warner 2014). Prior social VIs have included age, race, health, income, type of dwelling unit, employment of people living within a region, but are not specific to rural or agricultural vulnerability. In our study, we develop a VI that focuses on variables for which there is a clear causal pathway linking the condition to the ability of people in rural communities to prepare for, respond to, or adapt to climate stresses. Due to data availability, our analysis focuses on average levels of vulnerability at the municipality level (aggregate statistics of individual households within each municipality) in regions where coffee is grown.

The factors considered in our VI include average household and property conditions (age, knowledge and social capital, technology, and household and farm economy), as well as regional conditions (infrastructure and yields) in 2006 and 2017. The data source for the household and property conditions is the Brazilian Agricultural Census (table S2), while infrastructure comes from OpenStreeMap and yields from the Brazilian Agricultural Municipal Surveys. Variables included in the VI were selected based on their theoretical and empirical impacts on vulnerability or its sub-components, from previous studies (Kellerman 1983, Hahn et al 2009, Pandey and Jha 2012, Shah et al 2013, Garrett et al 2013). Specific justifications for each variable in the index are included in the SI.

To generate the VI from these six groups of variables, we use a composite index approach, whereby sub-components are first normalized and then averaged at the group level (i.e. major component) and then each group is averaged into a single index value. The composite approach is a standard way of achieving a single numerical value when the amount of data per group is unbalanced. It also allows the user to weight each group differently based on their theoretical importance. This approach was chosen over a PCA approach due to the small number of candidate variables and the low correlation between them. The SI discusses the justifications for the weightings used in the main text results, as well as the results of a sensitivity analysis using different weightings. The sensitivity analysis indicates that both weightings produce the same distribution of relative vulnerability over space, except that municipalities in Rio de Janeiro have an even higher relative vulnerability using an equal weights method over an approach that more heavily weights baseline yields and household assets (figure S3).

3.3. VI validation

4.1. Climate hazard

Since 1974, temperatures in Brazilian coffee growing municipalities have been increasing by ∼0.25 °C per decade (figures 2 (a)–(c)). Annual precipitation has been decreasing during the blooming and ripening period, and until 2002 during the harvesting period as well (figures 2 (d)–(f)). However, since 2002 precipitation during the harvesting period has been increasing (figure 2(f)). This recent trend may help offset the increasing temperature, sustaining the development of coffee beans and preventing cherries from ripening too soon during the dry season.

Zoom In Zoom Out Reset image size

Municipalities in the north of the study region (Bahia, northern Goiás, and Minas Gerais) have the highest mean temperature. Since 2010 mean temperatures in this region have frequently exceeded the optimal range for Arabica coffee (>23 °C). During the flowering period temperatures in all states increased by more than 1.2 °C. In Bahia, Minas Gerais, and São Paulo these large increases in temperature were coupled with large decreases in rainfall (>10% decrease).

Like past studies of annuals, we found that temperature and precipitation increase yields up to a point, but then detrimentally influence yields (table S1). This concave relationship is constrained to the blooming and ripening period. This supports the hypothesis that excessively high temperatures coupled with a lack of rainfall can inhibit both the initial flower budding and the development and maturing of ripening beans. During the harvest period, temperatures had a weakly convex relationship with yields (the linear term was significant, large and negative, while the squared term was significant and positive, but very small). This suggests non-linearly increasing benefits to cooler temperatures during the coffee harvest period, which contradicts established agronomic understanding of temperature-yield relationships and may reflect a limitation of the model's structural form. Precipitation had a monotonically positive relationship with yields.

The net impact of climate trends since 1974 have been negative overall, with the biggest impacts concentrated in Minas Gerais. The average yield loss ranged from 9% to 29% across the study region. The differences in climate hazard across individual municipalities within each state are separated into quantiles and mapped in figure 3 for climate trends from 1974 to 2006 and 1974 to 2017.

Zoom In Zoom Out Reset image size

4.2. Climate exposure

4.3. Climate vulnerability

4.4. Overall risk (overlap of climate hazard, exposure, and vulnerability)

We find that Minas Gerais and Rio de Janeiro have the highest mean risk under all classification systems due to the combination of high exposure, hazard, and vulnerability (figures 3 and 4). Paraná and Goaís have the lowest mean climate risk due to a combination of low exposure, hazard, and vulnerability.

Zoom In Zoom Out Reset image size

Cumulative risk decreased in all states except for Goiás and Minas Gerais, largely due to decreases in exposure. This decline in risk is accentuated when equal length intervals are used (shown in figure 5). Since each component of risk is a relative score, this means that the number of municipalities that fell into moderate to very high vulnerability and hazard classes in most states decreased over time. Decreases in mean climate exposure were partially offset by increases in vulnerability in several states. Risk decreased the most in Mato Grosso do Sul and Paraná through reductions in exposure.

Zoom In Zoom Out Reset image size

While prior climate zoning studies forecast major reductions in the area suitable for Arabica coffee production in Brazil by 2050 (Zullo et al 2006, 2011), our study shows that historical climate change is already having a substantial negative impact on yields. However, this climate hazard is not evenly distributed. Worryingly, it is concentrated in the Southeast of Brazil, including Minas Gerais, the region with the highest Arabica coffee production. Little adaptation has occurred in the form of reducing exposure by diversifying away from coffee production. In contrast, the states in Brazil that had the lowest climate hazard did reduce their exposure to climate hazards via agricultural diversification.

Differences in all aspects of climate risk mirror biophysical conditions and rural development levels. Many of the mountainous regions in Brazil, which have the highest climate risk, rely heavily on coffee production as a farming activitiy due to the high slope, lack infrastructure and rural services, and lack of other economic development opportunities (Watson and Achinelli 2008). Mountainous regions also tend to experience more dramatic shifts in climate than flatter areas, explaining higher climate hazards (Diffenbaugh and Giorgi 2012).

Yet, our municipality-level estimates of climate hazards and vulnerability may not accurately represent the conditions for all farms in each region. Across all data sources, the mismatch between sub-municipal areas and municipal averages is likely to be larger in regions with greater variation in altitude, orientation, slope, road density, and distance to municipal centers. Additionally, measurement error, particularly with respect to the agricultural censuses, may be larger in more remote regions where it may not be possible to interview all farmers. Another limitation of this work is that it relies heavily on existing understanding of the correlates of climate vulnerability in rural regions to construct the vulnerability index. In much of Brazil, and for coffee specifically, the causes of climate vulnerability are poorly understood. Future studies should supplement this research by analyzing the mechanisms underlying coffee farmers' vulnerability, including through in-depth field interviews.

Since livelihood risks are not constrained to a single crop, future work should examine climate risks across the whole livelihood portfolio of rural households, including off-farm activities. Climate hazards may damage people's ability to maintain or secure material assets and resources, as well as their ability to live a 'good life' in other ways, including non-material goals (Chambers 2013). However, it is also possible that climate change could create new opportunities for regions that become more favorable for coffee farming or new crops, such as cocoa (de Sousa et al 2019). Further in-depth fieldwork is urgently needed to investigate these issues in the Brazilian context.

In this study we sought to quantitatively assess the spatial and temporal variation in the climate risk of coffee communities in Brazil. Unlike past work, we measure climate hazard, exposure, and vulnerability independently and then identify where they overlap to increase overall risk. This approach allows us to assess which regions have the highest overall risk, as well as the major sources of that risk.

Our study finds that since 1974, temperatures in Brazilian coffee growing municipalities have been increasing by ∼0.25 °C per decade and annual precipitation has been decreasing during the blooming and ripening periods. This historical climate change has already resulted in reductions in coffee yield by more than 20% in the Southeast of Brazil. The Southeast, particularly Minas Gerais, is the largest coffee producing region in Brazil, so this high climate hazard translates into high overall climate risk for much of the country's core coffee producing regions. In the mountainous Brazilian highlands, where farmers are generally poorer and more disconnected from markets, the risk posed by high climate hazards and exposure is exacerbated by high vulnerability.

Our results provide useful information for the targeting of agricultural policies and climate planning. They indicate that federal and civil society efforts to prepare for climate change in coffee production regions should focus on Minas Gerais, where a majority of production occurs, rural economic dependence on coffee is very high, and climate hazards and vulnerability are the highest. Minas Gerais has already experienced high climate hazards from rising temperatures and declining precipitation in the blooming and ripening period. Improved coffee varieties, agricultural loans for irrigation and agroforestry might enable coffee farmers to maintain or improve their yield under climate hazard, while infrastructure development and capacity building within existing cooperatives could help increase access to higher value marketing opportunities to offset lower yields. Given its high climate vulnerability, but lower hazard, climate adaptation strategies for the more mountainous Eastern Minas Gerais could focus on broader development interventions, such as improved services, infrastructure, and market access. In contrast, climate adaptation strategies in Southwestern Minas Gerais would be better off focusing on diversifying production, off-farm income diversification, increasing access to irrigation, or expanding climate-related crop insurance.

We thank the three anonymous reviewers for their careful comments, which substantially improved the manuscript. This work would not have been possible without the generous contributions of our co-author Anthony Janetos, who sadly passed away earlier this year. We further thank Mark Friedl, Anne Short Gianotti, Robert Kaufmann, and Ian Sue Wing for their suggestions and comments on earlier versions of this research. This research was supported by the Fulbright Foundation NEXUS Program for the Western Hemisphere, the Summer Fellows program at the Boston University Fredrick S. Pardee Center for the Study of the Longer-Range Future and the Boston University Global Development Policy Center Land Use and Livelihoods Initiative.

The data that support the findings of this study are available from the corresponding author upon reasonable request.