https://orcid.org/0000-0001-5397-0201
Figures
Abstract
Global coffee production is at risk from synchronous crop failures, characterised by widespread concurrent reductions in yield occurring in multiple countries at the same time. For other crops, previous studies have shown that synchronous failures can be forced by spatially compounding climate anomalies, which in turn may be driven by large-scale climate modes such as the El Niño Southern Oscillation (ENSO). We provide a systematic analysis of spatially compounding climate hazards relevant to global coffee production. We identify 12 climate hazards from the literature, and assess the extent to which these hazards occur and co-occur for the top 12 coffee producing regions globally. We find that the number of climate hazards and compound events has increased in every region between 1980 and 2020. Furthermore, a clear climate change signature is evident, as the type of hazard has shifted from overly cool conditions to overly warm. Spatially compounding hazards have become particularly common in the past decade, with only one of the six most hazardous years occurring before 2010. Our results suggest that ENSO is the primary mode in explaining annual compound event variability, both globally and regionally. El Niño-like sea-surface temperatures in the Pacific Ocean are associated with decreased precipitation and increased temperatures in most coffee regions, and with spatially compounding warm and dry events. This relationship is reversed for La Niña-like signatures. The Madden Julian Oscillation also shows a strong association with climate hazards to coffee, with increased activity in the Maritime Continent related to a global increase in the number of cold or wet hazards and a decrease in the number of warm or dry hazards. With climate change projections showing a continued rise in temperatures in the tropics is likely, we suggest that coffee production can expect ongoing systemic shocks in response to spatially compounding climate hazards.
Citation: Richardson D, Kath J, Byrareddy VM, Monselesan DP, Risbey JS, Squire DT, et al. (2023) Synchronous climate hazards pose an increasing challenge to global coffee production. PLOS Clim 2(3): e0000134. https://doi.org/10.1371/journal.pclm.0000134
Editor: Liqiang Xu, Hefei University of Technology, CHINA
Received: August 21, 2022; Accepted: January 17, 2023; Published: March 8, 2023
Copyright: © 2023 Richardson et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The United States Department of Agriculture publishes annual coffee reports containing yield data at https://usda.library.cornell.edu/concern/publications/m900nt40f. Gridded monthly Berkely Earth Surface Temperature data is available from http://berkeleyearth.org/data/. The Global Precipitation Climatology Centre FD_M_V2020_100 data set can be accessed at https://opendata.dwd.de/climate_environment/GPCC/html/fulldata-monthly_v2020_doi_download.html. Monthly VPD is calculated with 2m temperature and dew-point temperature data from the ERA5 reanalysis product, available at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means. Our MJO indices are calculated from ERA5 hourly velocity potential data and can be retrieved from https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels. Monthly sea-surface temperature data from the Met Office Hadley Centre Sea Ice and Sea Surface Temperature data set is available from https://www.metoffice.gov.uk/hadobs/hadisst/. Maps in the relevant figures are made with Natural Earth. Free vector and raster map data at naturalearthdata.com. Natural Earth Terms of Use are available at https://www.naturalearthdata.com/about/terms-of-use/. The code used for this work is openly available at https://doi.org/10.5281/zenodo.7582943.
Funding: DR, DPM, JSR, DTS and CRT were supported by the Australian Climate Service https://www.acs.gov.au/. JK and VMB were supported by the Federal Ministry for the Environment, Nature Conservation and Nuclear Safety—International Climate Initiative (IKI) https://www.international-climate-initiative.com/en/. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
Coffee production supports the livelihoods of tens of millions of farmers in developing countries [1]. Annual coffee exports were worth an estimated USD 35.6 billion in 2018, a more than fourfold increase from 1991 [2]. The vast majority of the world’s coffee is made up of two species, Coffea arabica (Arabica) and Coffea canephora (robusta). Coffee, especially Arabica, is considered a sensitive crop, vulnerable to climate variability and change [3–5].
Broadly speaking, there are two ways in which the climate influences coffee cultivation. First, the local climatology dictates a region’s suitability to growing coffee. The optimal annual temperature range for growing coffee is commonly cited as between 18°C and 22°C (for Arabica), or between 22°C and 28°C (for robusta) [6–8]. Likewise for annual precipitation, optimal ranges are between 1400 mm and 2000 mm (Arabica), or between 2000 mm and 2500 mm (for robusta) [6]. These climatic values are derived from historical surveys on the distribution of these species and are used to infer the climatic conditions suitable for coffee (i.e. if the species is able to survive and reproduce). Substantial departures from these climatic ranges are taken as an indication that an area is either currently, or in the future under climate change, unsuitable for growing coffee [3].
Anthropogenic climate change is expected to alter the global distribution of coffee suitability. The area of land suitable for coffee cultivation may be reduced by up to 50% [3]. Globally, rising temperatures are the main driver of these projected losses, but regional studies indicate that changes in precipitation totals and seasonality is crucial [5]. Production could in theory be moved to more suitable regions, such as to higher altitudes. However, many farmers, and especially small land holders, often lack the resources and flexibility to relocate [2]. Such a large scale change in land use would result in significant socio-economic and environmental challenges.
The other way that climate affects coffee cultivation is through the impacts of inter-annual variability on the annual production cycle and on yields (as opposed to suitability studies, which are based on the presence and/or absence of coffee farms). During any given year, climate hazards such as heatwaves, droughts, frosts and floods can each affect coffee yield [8–14]. Sub-optimal temperatures and precipitation deficits have negative effects on yield and bean quality, and climate acts as a control on pests and diseases [5]. The timing is also important, as the vulnerability of coffee to climate variables changes depending on the stage of the plant’s life cycle [8, 15].
As with other crops, a systemic risk to the global coffee trade is posed by synchronised crop failures. These failures are characterised by large-scale yield deficits, and can arise as a result of widespread, spatially compounding climate anomalies [16–20]. For coffee productivity, however, the impacts of climate variability are typically analysed on national or regional spatial scales [5, 8, 10, 13, 14]. On a global scale, the historical variability and changes in the frequency of spatially compounding events that affect coffee production is unknown.
It is likely that large-scale climate modes play a role in forcing spatially compounding events [21–24]. However, the extent to which climate modes explain spatially occurring climate hazards important for coffee production has not been investigated. The El Niño Southern Oscillation (ENSO) is the dominant mode of Pacific Ocean sea-surface temperature (SST) variability that influences tropical precipitation and temperature [25, 26]. ENSO has been directly linked to world coffee production and prices [27–30], and has been implicated in forcing synchronous crop failures [16, 17, 19]
Aside from ENSO, previous studies have highlighted the teleconnection of modes of variability in the Atlantic and Indian Oceans to surface conditions in the tropics. The Indian Ocean Dipole (IOD), the Atlantic Niño, and the Tropical North and South Atlantic (TNA and TSA) modes have varying influence across coffee-growing regions [22, 31–34]. The Madden-Julian Oscillation (MJO) is the primary mode of intraseasonal variability in the tropics that modulates global convection [35]. MJO activity plays an important role in driving tropical precipitation and temperature variations [36, 37], including during crop flowering seasons [38]. To what extent the MJO affects climate hazards important for global coffee production is an open research question.
In this article, we present analysis of spatially compounding climate hazards relevant to global coffee growing area suitability and production. We consider a range of coffee-producing regions, each with their respective production cycles. From the literature, we identify the climate processes that influence coffee productivity and a region’s suitability for cultivation, which we term ‘climate hazards’. These broad-ranging factors, incorporating various hazards, regions and seasons, are synthesised in an approach that allows us to concisely explore the changes over time of spatially compounding events.
We also assess the relationship between these compound events and potential drivers of the tropical climate, and highlight which climate modes are most important in driving compound event occurrences. The article is structured as follows. In Section 2, we detail the data and methods used. We provide a description and discussion of the results in Section 3. The article concludes in Section 4 with a summary of the presented work.
2 Data and methods
2.1 Selecting coffee regions
We select the top 12 coffee-producing countries according to 2019 yield data from the United States Department of Agriculture (Fig 1a) [39]. These countries together account for 90% of global production. While some countries grow both Arabica and robusta, we only consider the primary species for the regions here. Due to Brazil’s size and importance to global coffee supply, we split the country into northern and southern Brazil. Southern Brazil is the top global grower of Arabica (53% of total Arabica supply), while northern Brazil is second only to Vietnam in producing robusta (26% of robusta supply; Fig 1b).
(a) Percentage of 2019 coffee yield by primary coffee species per region. Robusta regions are to the right of the red line, plus northern Brazil. Dashed lines indicate regions used to calculate ocean modes indices. (b) Flowering and growing seasons, with percentage of 2019 coffee yield by species shown on the right-hand axis. The map base layer is available from Natural Earth at https://www.naturalearthdata.com/downloads/110m-physical-vectors/110m-coastline/.
The production cycle of coffee varies by region. We consider two stages of production: the flowering season and the growing season. These seasons are shown in (Fig 1b). Colombia and Uganda have two distinct production cycles and we represent both here. We estimated the flowering and growing (i.e. cherry and fruit development) seasons for each country based on the literature [13, 40–43].
2.2 Identifying climate hazards
From the literature, we select 12 climate hazards, six for each coffee species (Table 1). These hazards are defined as when some variable surpasses a biophysical threshold. Growing season average temperature and annual precipitation ranges reflect the suitability of a region to growing coffee [6–8]. Growing season vapour pressure deficit (VPD) and average maximum daily temperature thresholds are chosen as hazards important for Arabica production [44]. For robusta production, too-cold minimum daily temperatures in the flowering season, and too-warm minimum daily temperatures in the growing season, have been identified as important variables [13]. As we will show, the hazards do not necessarily reflect extreme climate conditions. Instead, the thresholds represent sub-optimal conditions, where if surpassed could have a detrimental effect on coffee growing suitability and yield.
2.3 Calculating climate hazard events
We obtain monthly averages of daily maximum, daily minimum, and daily mean temperatures from the Berkeley Earth Surface Temperature data set [45], with a 1°x1° spatial resolution. Monthly precipitation data on the same grid are from the Global Precipitation Climatology Centre FD_M_V2020_100 [46] data set. Monthly mean VPD (with units of kPa) is calculated as
(1)
where c1 = 0.611 kPa, c2 = 17.5, c3 = 240.978°C, T is monthly mean temperature (°C) and Td is monthly mean dew-point temperature (°C) [47]. VPD is calculated using the ERA5 reanalysis product [48], with a spatial resolution of 0.25° × 0.25°. For all data we consider the period 1979 through 2020.
The temperature variables and VPD are aggregated to flowering and growing season means according to the months shown in Fig 1b. Precipitation is aggregated by accumulating the data over a 12 month period ending on the last month of the growing season. Where seasons cross over calendar years, they are labelled according to the year in which the last month falls. This process yields annual data for the flowering and growing seasons for the 41-year period 1980–2020.
Grid cell-level climate hazard events are defined as when a variable surpasses the thresholds listed in Table 1. Only grid cells that correspond to coffee growing regions are considered. We use a mask of coffee production intensity in 2010 provided by the Spatial Product Allocation Model with a spatial resolution of 0.083° × 0.083° [49]. This mask is aggregated to the spatial resolution of the relevant variable to exclude grid cells that do not grow coffee. A lower-resolution grid cell is counted as coffee-producing should any grid cell from the high resolution mask fall within it. Climate hazard events at the grid cell-level will be used to assess the regions’ overall (climatological) susceptibility to the hazards, and in the calculation of region-level events.
As we are concerned with widespread climate events likely to have impacts on national and global scale production, we consider a region-level climate hazard to occur when the proportion of a region’s area that surpasses the threshold is greater than one standard deviation above the 41-year mean. These region-level hazards are therefore indicative of systemic risk to coffee production, rather than local risk.
We use a relative threshold here, rather than an absolute threshold (e.g. exceeding 50% of the region’s area), because we are interested in changes in climate hazard susceptibility. There is substantial variation in the regions’ susceptibility to the climate hazards (Table A in S1 Text). For example, the mean areal percentage of Colombia 2 experiencing too-high precipitation totals is 70%. An absolute threshold of, say, 50% would result in an event every year for Colombia 2, without showing any potential changes (for example from 60% to 80% between 1980–2020). It is this reasoning that led us to choose a relative measure.
A consequence of our choice of threshold is that a region-level event can be triggered by a small number of grid cells. In regions with little year-on-year variation, this standard deviation threshold is small (Table A in S1 Text), with only a few additional grid cells required to exceed it. However, a small number of grid cells difference year-to-year has much bigger implications for small regions than large regions. In the largest region, northern Brazil, a grid cell accounts for at most 0.47% of total coffee-growing land. In the smallest region, Nicaragua, this figure is 12.6%.
We assessed the sensitivity of our results to the choice of threshold by using alternative definitions, of the 41-year mean areal proportion plus 5%, 10% and 20% of the region’s area (as opposed to the additional increment of the standard deviations shown in Table A in S1 Text). While our results did show some sensitivity to the choice of threshold, those from the 5% and 10% increments do not change the overall message. A 20% increment is much stricter than the other thresholds, and did not yield enough region-level events from which to draw any conclusions.
To test whether there has been a temporal shift in the type (e.g. from cold to warm) and frequency of climate hazards, we test for a monotonic trend in the global number of region-level hazards per year. Following [24], we calculate the Mann-Kendall test statistic, S, [50, 51] on the observed data as
(2)
where
(3)
where X = {X1, X2, …, XT} is the time series of annual hazard events. The test statistic is compared with 10,000 test statistics computed on samples generated using a circular block bootstrap procedure [52]. Block bootstrapping is used to try and mitigate the effects of serial correlation, which is a common feature of meteorological time series and which violates the assumption of the Mann-Kendall test that the data are independent and identically distributed. The block size is estimated using
(4)
where n′ is the effective sample size given by
(5) n is the sample size and ρ1 is the lag-1 autocorrelation coefficient.
The procedure is “circular” because the first B values of the time series are appended to the end of the time series. This avoids the problem of under-sampling the tails of the time series when bootstrapping.
The block size for the hazard frequencies is B = 4 years. We deem the observed test statistic as statistically significant if it is smaller the the 5th percentile or greater than the 95th percentile of the distribution of bootstrap test statistics.
2.4 Climate mode indices
We explore the role of six tropical climate modes in driving the occurrence of climate hazards to coffee production. For this part of the analysis, we first remove the climate change signal by linearly detrending variables used to calculate the climate hazards (temperature, precipitation and VPD) and climate mode indices. As the biophysical thresholds can’t be applied to these detrended data, we use relative thresholds of one standard deviation from the mean (above or below the mean as appropriate for the hazard).
For five ocean modes, we use SST observations from the Met Office Hadley Centre Sea Ice and Sea Surface Temperature data set [53] for the period 1979–2020. SST anomalies are computed by removing the mean of the whole period. We represent ENSO with the Niño3.4 index, defined as the SST anomaly averaged over 5°N-5°S, 120°-170°W. The Dipole Mode Index (DMI) quantifies the Indian Ocean Dipole (IOD), computed as the difference between SST anomalies averaged over western (10°N-10°S, 50°-70°E) and south-eastern (0°-10°S, 90°-110°E) regions in the Indian Ocean. We use three indices of the tropical Atlantic Ocean. The Tropical North Atlantic index (TNA) is the SST anomaly averaged over 25°N-5°N, 55°-15°W. The Tropical South Atlantic index (TSA) is the SST anomaly averaged over 0°N-20°S, 30°-10°W. The Atlantic Niño index is the SST anomaly averaged over 5°N-5°S, 20°-0°W.
To represent the MJO, we use eight indices defined as the number of days per month that the MJO is in each of its eight phases, expressed as anomalies relative to the long-term average. We use the method described in [54] to obtain a daily representation of the MJO using ERA5 velocity potential data. For a full description see the text, Figs A and B in S1 Text. The eight MJO phases represent the west-to-east path of the MJO through four general regions (Fig B in S1 Text): the western Hemisphere and Africa (phases 1 and 8), the Indian Ocean (2 and 3), the Maritime Continent (4 and 5), and the western Pacific (6 and 7). Our monthly MJO indices for each phase are denoted as MJOi, for i = 1, …, 8.
The monthly ocean and MJO indices are averaged over the growing seasons for each region. We are therefore assessing the concurrent relationship between climate hazards and potential drivers. Lagged relationships are not considered, as the growing seasons are long enough (minimum four months) to assume that any influence of the drivers on the hazards would manifest over those time scales. We do not consider the Tmin,fl climate hazard as it is relatively less important than growing season climate hazards [13].
2.5 Regression analysis
To explore which climate modes are the most important in explaining compound climate hazards for each coffee growing area, we fit regression models to the number of region-level hazards per year and region. Following the hazard classifications provided in Table 1, we assign ‘warm’ or ‘dry’ events a value of 1, and ‘cold’ or ‘wet’ events a value of −1, and sum for each region and year. For example, if a region experiences one dry and one warm event during a particular year, that year will have a value of 2. A value of 2 could also be obtained by the occurrence of three warm or dry hazards plus one cold or wet hazard. A value of 0 indicates that either no hazards occurred, or there was one cold or wet hazard, plus one warm or dry hazard. We will show that these instances of cancellation are few.
The response variable, y, is therefore an integer in the interval [−2, 4] for Arabica regions, and in the interval [−2, 3] for robusta regions. While y is best described as ordered categorical data, the number of categories (six or seven) suggests that treating the data as continuous is appropriate [55]. We therefore fit a Gaussian generalised linear model (GLM) with the identity link (equivalent to multiple linear regression).
To test the continuous data assumption, we also fit a GLM with a multinomial conditional distribution and the cumulative logit link function (i.e. ordinal regression). We describe results only for the Gaussian GLM, as model performance (see below) is worse when using ordinal regression for all regions except Indonesia and Uganda 1, but with similar sets of explanatory variables in the best-fitting models.
For each model, we consider the five ocean mode indices and the eight MJO indices as explanatory variables. We fit every possible combination of explanatory variables, except if any pair of explanatory variables have an absolute Spearman correlation greater than 0.7 to avoid potential collinearity issues [56]. The strength of the correlation changes depending on the season, and so depends on the growing season of each region. The total number of unique sets of explanatory variables ranges between 1,042 (Mexico) and 4,082 (Peru). Prior to model fitting, all data are standardised by dividing by their standard deviation in time. This ensures the relative importance of the regression coefficients can be assessed [57]. We consider the best models as those that minimise the small-sample corrected Aikake’s Information Criterion (AICc).
3 Results and discussion
3.1 Climatology of climate hazards
As mentioned in Section 2.2, the biophysical hazard thresholds are not always ‘extreme’ values of each variable. The historical frequency of each climate hazard varies greatly by region (Fig 2). In some regions, some hazards have never occurred in the period 1980–2020. The vast majority of coffee regions never experience too-cold growing season temperatures, for example (Fig 2c). Conversely, there are regions in which a hazard occurs every year, such as growing season temperatures over 22°C in southern Brazil (Fig 2d). In these cases, the ‘hazard’ is just a feature of the climate that must be managed to attain the best possible yields given those conditions. For example, regions with too-high temperatures often use shading as a management technique [58], while irrigation is employed to mitigate water stress in regions that do not receive optimal precipitation [59, 60].
Number of years between 1980 and 2020 during which climate variables surpass biophysical coffee thresholds for (a) high growing-season VPD (VPDgr) or low flowering-season minimum temperature (Tmin,fl), (b) high growing-season maximum temperature (Tmax,gr) or minimum temperature (Tmin,gr), (c) and (d) low and high growing-season mean temperature (Tgr), (e) and (f) low and high annual precipitation (Pan). Robusta regions and corresponding hazard definitions are to the right of the red line, plus northern Brazil. The map base layer is available from Natural Earth at https://www.naturalearthdata.com/downloads/110m-physical-vectors/110m-coastline/.
Southern Brazil exhibits the greatest spatial variability, with varying degrees of susceptibility to VPDgr, Tmax,gr and low Pan (Fig 2a, 2b and 2e). As by far the largest grower of Arabica, global coffee production is heavily dependent on favourable climate conditions in southern Brazil, or on effective management of poor conditions.
Due to its higher inter-annual variability, precipitation could be the most important variable that affects large-scale coffee production due to the spatial heterogeneity of Pan hazards. This is apparent for too-dry conditions in northern and southern Brazil (Fig 2e) and for too-wet conditions in Colombia, Peru and Indonesia (Fig 2f). By comparison, the majority of regions either never experience temperature-based hazards, or they are a feature of the climate.
The north of southern Brazil and eastern Ethiopia are the regions most susceptible to climate hazards, experiencing close to the maximum four hazards every year (Fig 3). In general, everywhere experiences at least one hazard per year on average. Robusta regions appear to suffer from fewer hazards per year on average, though this may be due to the differences in hazard types. As we will show later, it is not uncommon for Arabica regions to experience VPDgr and Tmax,gr hazards concurrently, but it is rare for robusta regions to have too-cold minimum temperatures in the flowering season followed by too-warm minimum temperatures in the growing season (i.e. Tmin,fl < 15.8°C followed by Tmin,gr > 18.6°C). Arabica is native to Ethiopia, and so unsurprisingly parts of this country do not usually experience any hazards during the year.
Robusta regions and corresponding hazard definitions are to the right of the red line, plus northern Brazil. The map base layer is available from Natural Earth at https://www.naturalearthdata.com/downloads/110m-physical-vectors/110m-coastline/.
3.2 Changes in climate hazard and compound event occurrences
Changes in temperature-related hazards show a clear climate change signal. Arabica regions are now much more likely to experience widespread high maximum temperatures and VPD (Fig 4a). Until around 2000, robusta regions were susceptible to cold overnight temperatures in the flowering season. Since then, these hazards have become less common, replaced by warm minimum temperatures in the growing season. A similar picture is evident in growing season mean temperatures across all regions, with a clear shift over time from below- to above-optimal temperatures (Fig 4b). In contrast, changes in sub-optimal precipitation totals do not display a clear trend (Fig 4c).
Time series show the occurrence of climate hazards for each region, for (a) growing season VPD (VPDgr) and maximum temperature (Tmax,gr; Arabica), and minimum temperature in the flowering (Tmin,fl) and growing (Tmin,gr) seasons (robusta), (b) growing season mean temperature (Tgr) and (c) annual precipitation (Pan).
Region-level hazard events occur when a relatively large proportion (here, above one standard deviation) of the region’s grid cells experience a hazard. This definition means that there are no region-level events in regions for which the hazard is an ever-present feature of the climate. One example is for southern Brazil, which has never suffered a region-level event for Tgr (Fig 4b). This is because no grid cell has ever had growing season temperatures below 18°C (Fig 2c). In addition, as virtually all grid cells are above 22°C every year (Fig 2d), there is never a situation in which the proportion of the region experiencing the hazard is large enough to satisfy the definition of a region-level event. Changes in these events therefore reflect whether a region is becoming more or less susceptible to a hazard, rather than the region’s general suitability for coffee cultivation.
The annual number of hazards globally has increased since 1980 (Fig 5). Most regions, with Uganda and India being the exceptions, appear to have experienced a greater number of hazards in the past decade. Since 2010, there have been five years during which at least 20 hazards occurred across all regions, compared to just once prior, in 1998 (top bar plot of Fig 5). This implies a greater risk to large-scale coffee supply from spatially compounding climate hazards. In 2016, for example, every region experienced at least one hazard, with most of the top growing regions (northern and southern Brazil, Colombia and Indonesia) experiencing multiple hazards simultaneously.
The main panel shows the number of hazard events per region and year. Shading indicates whether the majority of hazards are ‘warm or dry’ (brown) or ‘cold or wet’ (green) according to Table 1. On four occasions, one hazard from each of these classifications occurred (pink). The right-hand bar plot shows the number of hazards per region over the whole period. The top bar plot shows the number of hazards per year across all regions. These bar plots are shaded according to whether they are ‘warm or dry’ or ‘cold or wet’.
The climate change signal evident for individual hazards is stark for compound events, with a clear shift towards warm or dry hazards. This shift is supported by trend test results. We find that the Mann-Kendall test yields a statistically significant (p = 0.0003) upward trend for the annual number of warm or dry hazards per year, and a significant (p = 0.001) downward trend for the number of cold or wet hazards per year. With the lack of evidence for trends in precipitation hazards, the increasing number and shift in the type of hazard is driven by increasing temperatures (and related increases in VPD).
Southern Brazil, the world’s most productive Arabica producing area, has experienced among the fewest climate hazards, with only seven years featuring hazards between 1980 and 2013 (right-hand bar plot of Fig 5). This implies that, by and large, the majority of coffee producers in southern Brazil have not often had to adapt cultivation practices or implement mitigation strategies, at least for the hazards analysed here.
However, the region has experienced a concerning number of hazards since 2014. Since then, only two years have been hazard-free, with four years featuring multiple concurrent hazards. If the past seven years are indicative of the future, farmers may have to adapt to a hotter and drier climate to avoid negative impacts on coffee production.
3.3 Drivers of climate hazards and compound events
Fig 5 suggests that years with high numbers of climate hazards might be related to strong ENSO events. In particular, there were significant El Niño events in 1998, 2015, 2016 and 2019, and these years also featured high numbers of warm or dry hazards. To try and isolate any possible influence of the climate modes on the hazards from the climate change signal, we analyse their relationship using detrended data. Time series of the detrended region-level hazards are shown in Fig C in S1 Text.
Of the tropical ocean modes, ENSO is the most strongly correlated with precipitation and temperature (Fig 6). El Niño-like conditions favour warmer and drier conditions for every region except Southern Brazil, in which wetter conditions are more likely. The IOD has a relatively weak association with surface conditions. The strength of the IOD typically peaks in August to October, but the weak teleconnection is not to do with timing. This is because the growing seasons of each country, including those near the Indian Ocean where a strong teleconnection might be expected, span a wide array of months. The TNA exhibits similar (but weaker) relationships to surface variables as ENSO. The TNA and Atlantic Niño have similar teleconnection patterns to each other, and to the IOD. The exception is for Mexico, Central America and Colombia, for which the correlation sign is opposite (Figs D and E in S1 Text).
Spearman correlation between detrended growing season ENSO (measured by Niño3.4), Indian Ocean Dipole (DMI), Tropical North Atlantic index (TNA), Madden-Julian Oscillation indices for phases 1 and 4 (MJO1 or MJO4), and detrended growing season mean temperature (left column) or annual precipitation (right column). Robusta regions are to the right of the red line, plus northern Brazil. The map base layer is available from Natural Earth at https://www.naturalearthdata.com/downloads/110m-physical-vectors/110m-coastline/.
The MJO is quite strongly related to variations in temperature and precipitation, particularly in northern Brazil and Indonesia (Fig 6g–6j). The magnitude of these correlations does not vary greatly as a function of the MJO phase (also see Figs F and G in S1 Text), but the direction does. When the MJO is in phase 1 more often than normal, the majority of regions experience warmer temperatures and reduced precipitation, except in southern Brazil. Conversely, anomalously more frequent phase 4 days are linked to cooler and wetter conditions.
Our results do not feature the typical MJO teleconnection dipole pattern, whereby precipitation is enhanced over one half of the globe and suppressed in the other [35–37]. For example, we might expect MJO4, which indicates the extent to which the MJO is active (promoting precipitation) in the Maritime Continent, to be positively correlated with precipitation in Indonesia and negatively correlated with precipitation in northern Brazil. However, we see a positive correlation almost everywhere, including northern Brazil. This behaviour is potentially because the growing seasons over which we aggregate our MJO indices are much longer (four to nine months) than the time the MJO takes to traverse the planet (one to two months).
The El Niño year of 1998 is even more apparent in the detrended data than in the original series (Fig 7 and Fig C in S1 Text). Over all regions, 40 hazards occurred, 37 of which were warm or dry. As well as El Niño, that year saw positive values for the IOD, Atlantic Niño, TNA and TSA indices. In addition, there was a higher frequency of MJO activity near Africa, at the expense of activity over the Maritime Continent and western Pacific.
(a) Number of hazards per year across all regions. Circles and triangles denote years that featured above average numbers of ‘cold or wet’ and ‘warm or dry’ hazards, respectively. (b) Mode indices averaged across all regions’ growing seasons. The indices are standardised by dividing by their standard deviation in time. The mode indices are for Niño3.4 (ENSO), the Indian Ocean Dipole Mode index (IOD), the Atlantic Niño index (Atl. Niño), the Tropical North and South Atlantic indices (TNA and TSA) and the Madden-Julian Oscillation indices for each of the eight phases (MJO1 through MJO8).
Other years with a strongly positive Niño3.4 index, such as 1983, 1987, 2015 and 2016, also coincide with larger numbers of warm and dry climate hazards. Conversely, La Niña-like years such as 1989, 1999, 2000, 2008 and 2011 may be associated with some of the highest totals of cold and wet hazards.
As with so many studies analysing the ENSO teleconnection to surface climate, there is no one-to-one mapping of the ENSO phase to the climate hazards. The resultant Niño3.4 index after averaging across years with anomalously high numbers of warm/dry or cold/wet hazards (triangles and circles in Fig 7a) is around 0.6 or −0.7, respectively (Fig 8a). While the magnitude of these anomalies may seem small, they are statistically significant (according to a circular block bootstrap procedure carried out in the same manner as for the trend tests). As such, there is a clear association between ENSO phase and strength with the number and type of annual hazards (Fig 8b). An El Niño-like SST pattern tips the odds in favour of an increased number of warm or dry hazards globally, and vice versa for La Niña signature.
(a) Mean of standardised climate mode indices over years that featured above-average numbers of ‘warm or dry’ or ‘cold or wet’ hazards. Black circles indicate statistical significance. (b)-(e) Hazards per year plotted against ENSO (Niño3.4), Madden-Julian Oscillation phase 1 (MJO1), TNA and Madden-Julian Oscillation phase 4 (MJO4) indices, shaded according to ‘warm or dry’ (brown) and ‘cold or wet’ (green) hazards.
Of the remaining ocean modes, the TNA has the highest magnitude anomalies on spatially compounding years, with an average anomaly of −0.38 during cold or wet years (Fig 8a). However there is not an obvious distinction in TNA index magnitude or phase between cold/wet and warm/dry years (Fig 8d). Given that a cooler TNA has been associated with June-September droughts across tropical regions [22], we might have expected a negative relationship between the TNA index and the number of warm or dry hazards. However, we find a potentially opposite relationship, although it is weak and likely influenced by the outlier year of 1998. The differences in results we find compared to [22] are perhaps not surprising due to the different regions and seasons analysed, plus we consider temperature variables as well as precipitation.
After ENSO, the MJO exhibits the strongest changes in behaviour during spatially compounding years. When these compound years are characterised by cold and wet hazards, the MJO is more active than usual in the Maritime Continent (phases 4 and 5) at the expense of activity in the western Hemisphere and Africa (phases 1 and 8; Fig 8a, 8c and 8e). During years where hot and dry hazards dominate, this activity is reversed.
These results implicate ENSO, and to a lesser extent the MJO, as the climate modes most influencing global climate hazards important for coffee production. However, as ENSO and the MJO are correlated it is difficult to isolate the effects of each of these climate modes. Still, while the overall relationship between ENSO and the MJO is difficult to ascertain, El Niño events have been shown to modulate MJO amplitude in boreal winter [35, 61].
For the growing seasons considered here, we find Niño3.4 to be relatively strongly correlated with the MJO in phases, 1, 4 and 8 (Spearman correlation between |0.48| and |0.94|, depending on the growing season), and modestly correlated with phases 3, 5 and 7 (between |0.38| and |0.84|; Fig H in S1 Text). The MJO indices are correlated with each other, as expected given their construction. The TSA index is positively correlated with the Atlantic Niño index, which is unsurprising given the spatial overlap (Fig 4a).
We use regression analysis, excluding strongly correlated mode indices as explanatory variables, to identify the relative importance of the climate modes in explaining hazard frequencies for each region. In general, we find the Gaussian GLM performs reasonably well, as indicated by the apparent normality of the residuals (Figs I and J in S1 Text).
ENSO is the most important mode, as it has the largest absolute standardised coefficient, for explaining climate hazard occurrences in tropical South American regions (i.e. excluding southern Brazil) and for Indonesia (Fig 9). ENSO is also an important predictor for Vietnam and Mexico, although no more or less important than the IOD (for Vietnam and Mexico) and the Atlantic Niño, phase 1 MJO and phase 6 MJO indices (for Vietnam only). In some regions for which ENSO is not selected in the models (Ethiopia, Guatemala, Uganda 2 and India), MJO phases that are strongly correlated to Niño3.4 are present.
Standardised regression coefficients of explanatory variables for the best-performing Gaussian GLM models. Positive coefficient values indicate a positive relationship between the explanatory variable and the number of warm or dry hazards. Negative coefficient values indicate a positive relationship between the explanatory variable and the number of cold or wet hazards. The explanatory variables are Niño3.4 (ENSO), the Indian Ocean Dipole Mode index (IOD), the Atlantic Niño index (Atl. Niño), the Tropical North and South Atlantic indices (TNA and TSA) and the Madden-Julian Oscillation indices for each of the eight phases (MJO1 through MJO8).
Southern Brazil being relatively less influenced by ENSO may be important in mitigating systemic risk to global coffee supply. Colombia, Vietnam, northern Brazil and Indonesia are also major producers, yet are susceptible to climate hazards during El Niño-like conditions. With southern Brazil’s production power, it may be capable of making up for other regions’ shortfalls during El Niño events.
Despite the weak correlation of the DMI to temperature and precipitation (Fig 6), and its ambiguous relationship with spatially compounding hazard years (Figs 7 and 8a), the regression models imply the IOD is important in explaining climate hazard frequency variability in seven of the 15 regions (Fig 9). A positive IOD is associated with higher numbers of cold or wet hazards near the eastern Indian Ocean (Vietnam and India), plus Peru. The central and north American regions, on the other hand, may experience an increase in warm or dry hazards during a positive IOD.
4 Conclusions
We analysed the climatology, changes and drivers of climate hazards to global coffee production and cultivation suitability between 1980 and 2020. Focusing on 12 hazards identified from the literature, we showed that there is substantial variation in which hazards are experienced by the top 12 coffee-producing nations. A particular hazard may never have occurred in some regions. Conversely, in some places a ‘hazard’ may be an ever-present feature of the climate, implying cultivation practices such as irrigation or shading are used to ensure optimal productivity. Major Arabica regions in the far southeast of Brazil and southwest Ethiopia are amongst the least susceptible regions to climate hazards.
We find that spatially compounding climate hazards have increased in frequency since 1980, implying a greater systemic risk to global coffee production. The occurrence of these spatially compounding events has become particularly acute over the past decade, with 5 of the 6 most hazardous years occurring since 2010. This change is driven by the increasing frequency of temperature and VPD hazards. Furthermore, in the first half of the period 1980–2020 regions were more prone to experiencing too-cold temperatures in the flowering and growing seasons. The current climate, however, is characterised by too-hot conditions in every region.
We find that some large-scale tropical climate modes have preferential states during spatially compounding years, despite the wide range of hazards, regions and seasons analysed. ENSO is the most important climate mode in explaining yearly variations in spatially compounding events. A warm-phase ENSO is associated with an increase in warm and dry climate hazards, and a cold-phase ENSO with an increase in cold and wet climate hazards. However, this is not a one-to-one relationship, as not all El Niño- or La Niña-like years are notable for the numbers of climate hazards globally.
The MJO is also strongly associated with spatially compounding events, with the strength and direction depending on the MJO phase. An MJO active in the western Hemisphere and Africa more often than usual typically coincides with an increase in warm or dry hazards globally, and vice versa for MJO activity in the Maritime Continent.
The relationship between these climate modes and the hazards is confounded by the relatively strong correlation between ENSO and the MJO. Once this correlation is accounted for, we find ENSO to be an important predictor of hazards in tropical South America, Indonesia and Vietnam. Elsewhere, one or more phases of the MJO are the most important mode.
This includes southern Brazil, which may be a key region in offsetting El Niño-related risks to global coffee supply. As the world’s top coffee grower, the apparent lack of relationship between ENSO and climate hazards in the region could help to dampen production shocks felt elsewhere during significant ENSO events.
Over the coming decades, average temperatures and VPD are projected to increase in response to radiative forcing [62], increasing the likelihood of experiencing the related climate hazards analysed here. Moreover, climate change may alter the variability and surface teleconnections of key modes such as ENSO [63–65]. Taken together, the recent changes and future projections of the climate suggest coffee production can expect ongoing systemic shocks as a result of sub-optimal growing conditions.
Supporting information
S1 Text. Supporting text, tables and figures for the manuscript.
https://doi.org/10.1371/journal.pclm.0000134.s001
(PDF)
Acknowledgments
This work was funded by the Australian Climate Service. The work was undertaken using resources from the National Computational Infrastructure (NCI), which is supported by the Australian government. We acknowledge the Pangeo community for the open-source tools used in the analyses here. We appreciated the insightful comments of two reviewers who helped improve the manuscript.
References
- 1. DaMatta FM, Rahn E, Läderach P, Ghini R, Ramalho JC. Why could the coffee crop endure climate change and global warming to a greater extent than previously estimated? Climatic Change. 2019;152(1):167–178.
- 2.
International Coffee Organization. The Value of Coffee: Sustainability, Inclusiveness and Resilience of the Coffee Global Value Chain. International Coffee Organization; 2020. Available from: https://www.internationalcoffeecouncil.com/cdr2020.
- 3. Bunn C, Läderach P, Ovalle Rivera O, Kirschke D. A bitter cup: climate change profile of global production of Arabica and Robusta coffee. Climatic Change. 2015;129(1):89–101.
- 4. Moat J, Williams J, Baena S, Wilkinson T, Gole TW, Challa ZK, et al. Resilience potential of the Ethiopian coffee sector under climate change. Nature Plants. 2017;3(7):1–14. pmid:28628132
- 5. Pham Y, Reardon-Smith K, Mushtaq S, Cockfield G. The impact of climate change and variability on coffee production: a systematic review. Climatic Change. 2019;156(4):609–630.
- 6.
Descroix F, Snoeck J. Environmental Factors Suitable for Coffee Cultivation. In: Coffee: Growing, Processing, Sustainable Production. John Wiley & Sons, Ltd; 2004. p. 164–177. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/9783527619627.ch6.
- 7. Zullo J, Pinto HS, Assad ED, de Ávila AMH. Potential for growing Arabica coffee in the extreme south of Brazil in a warmer world. Climatic Change. 2011;109(3):535–548.
- 8. Koh I, Garrett R, Janetos A, Mueller ND. Climate risks to Brazilian coffee production. Environmental Research Letters. 2020;15(10):104015.
- 9. DaMatta FM, Ramalho JDC. Impacts of drought and temperature stress on coffee physiology and production: a review. Brazilian Journal of Plant Physiology. 2006;18(1):55–81.
- 10. Tucker CM, Eakin H, Castellanos EJ. Perceptions of risk and adaptation: Coffee producers, market shocks, and extreme weather in Central America and Mexico. Global Environmental Change. 2010;20(1):23–32.
- 11. Hoyos N, Escobar J, Restrepo JC, Arango AM, Ortiz JC. Impact of the 2010-2011 La Niña phenomenon in Colombia, South America: The human toll of an extreme weather event. Applied Geography. 2013;39:16–25.
- 12. Bacon CM, Sundstrom WA, Stewart IT, Beezer D. Vulnerability to Cumulative Hazards: Coping with the Coffee Leaf Rust Outbreak, Drought, and Food Insecurity in Nicaragua. World Development. 2017;93:136–152.
- 13. Kath J, Byrareddy VM, Craparo A, Nguyen-Huy T, Mushtaq S, Cao L, et al. Not so robust: Robusta coffee production is highly sensitive to temperature. Global Change Biology. 2020;26(6):3677–3688. pmid:32223007
- 14. Kath J, Byrareddy VM, Mushtaq S, Craparo A, Porcel M. Temperature and rainfall impacts on robusta coffee bean characteristics. Climate Risk Management. 2021;32:100281.
- 15. Craparo A, Van Asten PJ, Läderach P, Jassogne LT, Grab S. Coffea arabica yields decline in Tanzania due to climate change: Global implications. Agricultural and Forest Meteorology. 2015;207:1–10.
- 16. Iizumi T, Luo JJ, Challinor AJ, Sakurai G, Yokozawa M, Sakuma H, et al. Impacts of El Niño Southern Oscillation on the global yields of major crops. Nature Communications. 2014;5(1):3712. pmid:24827075
- 17. Anderson W, Seager R, Baethgen W, Cane M. Trans-Pacific ENSO teleconnections pose a correlated risk to agriculture. Agricultural and Forest Meteorology. 2018;262:298–309.
- 18. Sarhadi A, Ausín MC, Wiper MP, Touma D, Diffenbaugh NS. Multidimensional risk in a nonstationary climate: Joint probability of increasingly severe warm and dry conditions. Science Advances. 2018;4(11):eaau3487. pmid:30498780
- 19. Anderson WB, Seager R, Baethgen W, Cane M, You L. Synchronous crop failures and climate-forced production variability. Science Advances. 2019;5(7):eaaw1976. pmid:31281890
- 20. Vogel J, Rivoire P, Deidda C, Rahimi L, Sauter CA, Tschumi E, et al. Identifying meteorological drivers of extreme impacts: an application to simulated crop yields. Earth System Dynamics. 2021;12(1):151–172.
- 21. Kornhuber K, Coumou D, Vogel E, Lesk C, Donges JF, Lehmann J, et al. Amplified Rossby waves enhance risk of concurrent heatwaves in major breadbasket regions. Nature Climate Change. 2020;10(1):48–53.
- 22. Singh J, Ashfaq M, Skinner CB, Anderson WB, Singh D. Amplified risk of spatially compounding droughts during co-occurrences of modes of natural ocean variability. npj Climate and Atmospheric Science. 2021;4(1):1–14.
- 23. Squire DT, Richardson D, Risbey JS, Black AS, Kitsios V, Matear RJ, et al. Likelihood of unprecedented drought and fire weather during Australia’s 2019 megafires. npj Climate and Atmospheric Science. 2021;4(1):1–12.
- 24. Richardson D, Black AS, Irving D, Matear RJ, Monselesan DP, Risbey JS, et al. Global increase in wildfire potential from compound fire weather and drought. npj Climate and Atmospheric Science. 2022;5(1):23.
- 25. Mason SJ, Goddard L. Probabilistic Precipitation Anomalies Associated with EN SO. Bulletin of the American Meteorological Society. 2001;82(4):619–638.
- 26. Lyon B, Barnston AG. ENSO and the Spatial Extent of Interannual Precipitation Extremes in Tropical Land Areas. Journal of Climate. 2005;18(23):5095–5109.
- 27. Ubilava D. El Niño, La Niña, and world coffee price dynamics. Agricultural Economics. 2012;43(1):17–26.
- 28. Cobon DH, Ewai M, Inape K, Bourke RM. Food shortages are associated with droughts, floods, frosts and ENSO in Papua New Guinea. Agricultural Systems. 2016;145:150–164.
- 29. Bastianin A, Lanza A, Manera M. Economic impacts of El Niño southern oscillation: evidence from the Colombian coffee market. Agricultural Economics. 2018;49(5):623–633.
- 30. Almeida Silva K, de Souza Rolim G, Borges Valeriano TT, da Silva Cabral de Moraes JR. Influence of El Niño and La Niña on coffee yield in the main coffee-producing regions of Brazil. Theoretical and Applied Climatology. 2020;139(3):1019–1029.
- 31. Nicholson SE, Kim J. The Relationship of the El Niño Southern Oscillation to African Rainfall. International Journal of Climatology. 1997;17(2):117–135.
- 32. Ashok K, Guan Z, Saji NH, Yamagata T. Individual and Combined Influences of ENSO and the Indian Ocean Dipole on the Indian Summer Monsoon. Journal of Climate. 2004;17(16):3141–3155.
- 33. Yoon JH, Zeng N. An Atlantic influence on Amazon rainfall. Climate Dynamics. 2010;34(2):249–264.
- 34. Jimenez JC, Marengo JA, Alves LM, Sulca JC, Takahashi K, Ferrett S, et al. The role of ENSO flavours and TNA on recent droughts over Amazon forests and the Northeast Brazil region. International Journal of Climatology. 2021;41(7):3761–3780.
- 35.
Woolnough SJ. Chapter 5—The Madden-Julian Oscillation. In: Robertson AW, Vitart F, editors. Sub-Seasonal to Seasonal Prediction. Elsevier; 2019. p. 93–117. Available from: http://www.sciencedirect.com/science/article/pii/B978012811714900005X.
- 36. Donald A, Meinke H, Power B, Maia AdHN, Wheeler MC, White N, et al. Near-global impact of the Madden-Julian Oscillation on rainfall. Geophysical Research Letters. 2006;33(9).
- 37. Schreck CJ III. Global Survey of the MJO and Extreme Precipitation. Geophysical Research Letters. 2021;48(19):e2021GL094691.
- 38. Anderson W, Muñoz ÁG, Goddard L, Baethgen W, Chourio X. MJO teleconnections to crop growing seasons. Climate Dynamics. 2020.
- 39.
USDA. 2019/20 Forecast Overview. United States Department of Agriculture (USDA); 2019. Available from: https://downloads.usda.library.cornell.edu/usda-esmis/files/m900nt40f/sq87c919h/8w32rm91m/coffee.pdf.
- 40. Cannell M. Physiology of the coffee crop. Coffee. 1985; p. 108–134.
- 41.
Wintgens JN, et al. Coffee: growing, processing, sustainable production. A guidebook for growers, processors, traders, and researchers. WILEY-VCH Verlag GmbH & Co. KGaA; 2004.
- 42. DaMatta FM, Ronchi CP, Maestri M, Barros RS. Ecophysiology of coffee growth and production. Brazilian journal of plant physiology. 2007;19:485–510.
- 43. Wagner S, Jassogne L, Price E, Jones M, Preziosi R. Impact of climate change on the production of coffea arabica at Mt. Kilimanjaro, Tanzania. Agriculture. 2021;11(1):53.
- 44. Kath J, Craparo A, Fong Y, Byrareddy V, Davis AP, King R, et al. Vapour pressure deficit determines critical thresholds for global coffee production under climate change. Nature Food. 2022;3(10):871–880.
- 45. Rohde R, Muller R, Jacobsen R, Muller E, Perlmutter S, Rosenfeld A, et al. A New Estimate of the Average Earth Surface Land Temperature Spanning 1753 to 2011. Geoinformatics & Geostatistics: An Overview. 2013;01(01).
- 46. Becker A, Finger P, Meyer-Christoffer A, Rudolf B, Schamm K, Schneider U, et al. A description of the global land-surface precipitation data products of the Global Precipitation Climatology Centre with sample applications including centennial (trend) analysis from 1901-present. Earth System Science Data. 2013;5(1):71–99.
- 47. Barkhordarian A, Saatchi SS, Behrangi A, Loikith PC, Mechoso CR. A Recent Systematic Increase in Vapor Pressure Deficit over Tropical South America. Scientific Reports. 2019;9(1):15331. pmid:31653952
- 48. Hersbach H, Bell B, Berrisford P, Hirahara S, Horányi A, Muñoz-Sabater J, et al. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society. 2020;146(730):1999–2049.
- 49. Yu Q, You L, Wood-Sichra U, Ru Y, Joglekar AKB, Fritz S, et al. A cultivated planet in 2010—Part 2: The global gridded agricultural-production maps. Earth System Science Data. 2020;12(4):3545–3572.
- 50. Mann HB. Nonparametric Tests Against Trend. Econometrica. 1945;13(3):245–259.
- 51.
Kendall MG. Rank correlation methods. Rank correlation methods. Oxford, England: Griffin; 1955.
- 52.
Wilks DS. Chapter 5—Frequentist Statistical Inference. In: Wilks DS, editor. Statistical Methods in the Atmospheric Sciences (Fourth Edition). Elsevier; 2019. p. 143–207. Available from: http://www.sciencedirect.com/science/article/pii/B9780128158234000055.
- 53. Rayner NA, Parker DE, Horton EB, Folland CK, Alexander LV, Rowell DP, et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. Journal of Geophysical Research: Atmospheres. 2003;108(D14).
- 54. Adames ÁF, Wallace JM. Three-Dimensional Structure and Evolution of the MJO and Its Relation to the Mean Flow. Journal of the Atmospheric Sciences. 2014;71(6):2007–2026.
- 55. Rhemtulla M, Brosseau-Liard PÉ, Savalei V. When can categorical variables be treated as continuous? A comparison of robust continuous and categorical SEM estimation methods under suboptimal conditions. Psychological methods. 2012;17(3):354. pmid:22799625
- 56. Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography. 2013;36(1):27–46.
- 57.
Gelman A, Hill J. Data analysis using regression and multilevel/hierarchical models. Cambridge university press; 2006.
- 58. Beer J, Muschler R, Kass D, Somarriba E. Shade management in coffee and cacao plantations. Agroforestry Systems. 1997;38(1):139–164.
- 59. Amarasinghe UA, Hoanh CT, D’haeze D, Hung TQ. Toward sustainable coffee production in Vietnam: More coffee with less water. Agricultural Systems. 2015;136:96–105.
- 60. Boreux V, Vaast P, Madappa LP, Cheppudira KG, Garcia C, Ghazoul J. Agroforestry coffee production increased by native shade trees, irrigation, and liming. Agronomy for Sustainable Development. 2016;36(3):42.
- 61. Hendon HH, Zhang C, Glick JD. Interannual Variation of the Madden Julian Oscillation during Austral Summer. Journal of Climate. 1999;12(8):2538–2550.
- 62. Yuan W, Zheng Y, Piao S, Ciais P, Lombardozzi D, Wang Y, et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. Science Advances. 2019;5(8):eaax1396. pmid:31453338
- 63. Cai W, Santoso A, Collins M, Dewitte B, Karamperidou C, Kug JS, et al. Changing El Niño Southern Oscillation in a warming climate. Nature Reviews Earth & Environment. 2021;2(9):628–644.
- 64. Wengel C, Lee SS, Stuecker MF, Timmermann A, Chu JE, Schloesser F. Future high-resolution El Niño/Southern Oscillation dynamics. Nature Climate Change. 2021;11(9):758–765.
- 65. McGregor S, Cassou C, Kosaka Y, Phillips AS. Projected ENSO Teleconnection Changes in CMIP6. Geophysical Research Letters. 2022;49(11):e2021GL097511.