Impacts of land use and land cover dynamics on ecosystem services in the Yayo coffee forest biosphere reserve, southwestern Ethiopia

Land management to increase food production while conserving the environment and associated ecosystem services (ESs) is one of the major development and research challenges of the 21st Century. Any land-use practice or change to obtain a particular ecosystem service affects the other ES positively or negatively. The dynamics of these changes is more marked in biodiversity hotspot areas like UNESCO registered Yayo coffee forest biosphere reserve in southwestern Ethiopia. We used a time series InVEST modeling framework to estimate six ESs and analyze their spatial and temporal dynamics due to land-use/cover change over the last 31 years. Pearson correlation coefficients and k-mean clustering were employed to analyze tradeoffs/synergies and to cluster ESs supply spatially. The analysis also considers land-use change impact in the three management zones (core, transition and buffer) of the Yayo biosphere area. The production efficient frontier is used to identify the optimal combination of ESs and to suggest where an increase of one ES is possible without decreasing the others. Mostly, the highest change is observed in the transition zone followed by buffer zones. Positive correlation (synergies) are observed between regulating ecosystem services. Negative correlations (tradeoffs) are observed between provision ecosystem services. The clustering analysis shows that the spatial ESs can be divided in two clusters (bundle): cluster 1 with “ High regulating ESs ” that can be characterized by core zone and some forest patches in the central part of the biosphere reserve, and cluster 2 with “ High provisioning ESs areas ’’ that can be char- acterized by cultivated lands at transition and buffer zones. The result shows that the existing ES pairs are far from the Pareto efficient combination(s), confirming that landscape optimization for ES bundles are rarely possible on the ground due to many reasons and indicating the need for well thought land restoration strategies and land management practices that are forest type and context specific.


Introduction
As the world population increases, demand for food will grow at an alarming rate, and meeting this food demand would require between 320 and 850 million hectares of additional agricultural land by 2050 (Fróna et al., 2019). In recent decades, the demand for growing agricultural products has been partly met by increasing cultivated land (Boserup, 2017). In areas where population pressure outweighs resource supply, there will be expansion of cultivation and grazing areas into forests, marginal lands and other land uses to satisfy food demand (Fróna et al., 2019). However, in the future, the efficiency of agricultural production and specific yields must be increased through intensifying farming systems with the application of inputs such as fertilizers and pesticides as it will not be possible to expand cultivated lands indefinitely (Lambin and Meyfroidt, 2011). In addition, there is serious competition for other uses of land besides needs for expanding agricultural areas that may even lead to conflicts. Expansion of agricultural land can lead to significant land-use and land-cover (LULC) changes (Tilman and Clark, 2014). Currently agriculture is responsible for 80% of water use globally (FAO, 2017). In many parts of the world, the expansion of urban and rural settlements has also shown high correlation with natural resource exploitation, habitat fragmentation and biodiversity loss (Lawler et al., 2014;Lyu et al., 2018). Agricultural intensification also causes negative impacts such as increased pollution of ground water and eutrophication of rivers and lakes (Sala et al., 2000;Matson et al., 1997) as well as loss in biodiversity. Severe environmental degradation due to LULC and unsustainable agricultural intensification can result in the loss of ecosystem services (ESs) that are essential to human wellbeing. Reconciling the two needs -increasing food production and conserving the environment and associated ESs has thus become a major development and research challenge of the 21st Century (Foley et al., 2011;Brussaard et al., 2010).
Ethiopia is emblematic of a conflict between securing food selfsufficiency through enhanced agricultural production and avoiding the loss ecosystem functions and ESs due to land degradation (Hurni et al., 2010). The country is hugely affected by land degradation, as 85% of the total land is said to suffer from moderate to very serious levels of land degradation, costing about $4.3 billion per year (Gebreselassie et al., 2016). Estimates of average soil losses range between 3.4 and 84.5 tons ha-1 yr-1 with maximum rates reaching 300 tons ha-1 yr-1 (Hurni et al, 2015;Gashaw, 2015). National level nutrient depletion rates were estimated to be 122, 13 and 82 kg ha-1 yr-1 for N, P, and K, respectively (Haileslassie et al., 2005). The soil loss and nutrient depletion result in decline in agricultural productivity. Degradation of natural habitat such as conversion of pristine forest to cropland and settlements is the primary cause of biodiversity decline and loss (Fuller et al., 2007).
The Country has put various efforts in place to tackle land degradation and its associated impacts, including soil and water conservation practices since the 1960s. Recently, encouraging results have been observed in terms of improvement of ESs due to land restoration Adimassu et al., 2017). The country has also put different guidelines and policy action to tackle land degradation and support land restoration efforts. Among these is the declaration and registration of biosphere reserves by UNESCO to protect two of the biodiversity hotspots of the world that Ethiopia hosts, namely: the Eastern Afromontane and the Horn of Africa hotspots.
The Yayo Coffee Forest Biosphere Reserve (hereafter referred to as Yayo) is located in southwestern parts of Ethiopia. The main objective of establishing the biosphere reserve is to protect habitat characteristics linked to biodiversity and ESs that are vital to the maintenance of human well-being (Gole, 2003). Though some studies show the encouraging contributions of such interventions (Duguma et al., 2019), there are still gaps in that local communities have not thus far got sustained benefits (Getahun and Keno, 2019). As a result, there is conflict between the local communities and the nearby biosphere reserve (Getahun and Keno, 2019). In addition, there is expansion of agricultural land, illegal cutting of trees, and deforestation in Yayo, which are mainly prevalent in the transition zones, and fears are growing that these could expand to the buffer and interior zones (Beyene, 2014;Getahun and Keno, 2019).
Land-use and land-cover changes and management options that are intended to enhance a certain ES can change the supply of another ES (Bennett et al., 2009). The emergent interaction can be desirable (synergy) or undesirable (tradeoff). In the former case, an increase in one ecosystem service is associated with an increasing supply of another while the latter is when a gain in one ecosystem service results in a decline of the supply of another. It is thus essential to develop land-use practices and management plans that can enhance multiple ESs simultaneously (Balbi et al., 2015;Foley et al., 2011). In addition, it is necessary to consider that trade-offs between the increasing needs of the local population and maintenance of ecosystem health are minimized. Balancing the provision of ESs with human development needs is thus necessary for sustainable conservation and development pathways.
Ecosystem service analyses have been increasingly used to analyze the impacts of environmental changes on societal and economic terms to provide rational land management decisions. Spatially explicit assessment of ESs is essential for informative land management. Thus, mapping of various ESs is required to incorporate ESs information into land use related policy development. Different approaches can be used to assess ESs at different scales (e.g., MIMES (Boumans et al., 2015); Co $ting Nature (Mulligan, 2015); InVEST (Sharp et al., 2014a), SolVES (Sherrouse et al., 2014); EVT (Briceno & Kochmer, 2014), TESSA (Peh et al., 2013); ARIES (Villa et al., 2009), PA-BAT (Dudley and Stolton, 2008)). Many studies have explored the intended use, flexibility, limitations and uncertainty of these ecosystem service models (Sharp et al., 2015;Hamel and Bryant, 2017;Dennedy-Frank et al., 2016;Redhead et al., 2016;Willcock et al., 2016;Malinga et al., 2015;Maes et al., 2012). In this study, we selected InVEST model due to its generalizability and publicity by scientific literature. ESs models are less used in Sub-Saharan Africa to guide decision making and policy implementation due to lack of required data and researchers capacity to run ESs models (Willcock et al., 2016). More specifically, the applications of models for ESs estimation are few in Ethiopia (Mengist and Soromessa, 2019). The overall aim of this study is to estimate tradeoffs and synergies between agricultural production and other ESs due to the LULC changes happening in Yayo. This is useful to produce evidence and develop guidelines for optimal future investment and spatial targeting of agricultural production, and environmental conservation and management efforts. This enables to explore land use options for enhancing food security while conserving biodiversity and maintaining environmental integrity in Yayo. The specific objectives are: (1) to examine LULC change in the Yayo between 1986 to 2017, (2) to assess spatiotemporal dynamics to identify hotspots and coldspots of selected ecosystem services, and (3) to understand tradeoffs and synergies between ecosystem services and develop production frontiers curve to identify optimal ES provision. The study followed a conceptual modelling approach that focuses on obtaining relative ES values, tradeoff and synergy analysis, rather than making accurate prediction.

Study area
The Yayo Coffee Forest Biosphere Reserve is located in the southwestern Ethiopia and has been designated as UNESCO Biosphere Reserve since 2010 (Fig. 1). The mean annual temperature is about 23.76 • C, and it is located in the rainiest areas of the country with annual precipitation of about 1625 mm/year (Mulatu and Getahun, 2018). The total area of the biosphere reserve is 1670 km 2 , covering about 2 administrative zones and 5 Woredas (third level administrative units in Ethiopia, equivalent to district). Yayo includes Eastern Afromontane Biodiversity Hotspot, and montane rainforest fragments with wild Coffee arabica populations and encompasses important bird areas of international significance. Following its designation as a UNESCO biosphere reserve, three management zones were created: core, buffer and transition zones. The core zone (27,733 ha) is where there is no human intervention and is fully protected. The buffer zone (21,552 ha) is where limited human intervention is possible only if compatible with conservation objectives of the reserve. The transitional zone (117,736 ha) is found adjacent to the buffer zone and includes agricultural lands, wetlands, grasslands, settlements, and forest fragments. In this area, minimal human activities such as cultivation around homesteads are allowed. Despite its importance in terms of biodiversity, forest coffee and ESs provision, the area is threatened by deforestation, land degradation, and conflict between land uses and users (Dorresteijn et al., 2017).

Land use, land cover, and change (LULC) mapping
Ecosystem service estimation and trade-off analysis require well differentiated and spatially explicit LULC datasets. In this study, contextbased LULC mapping and change detection approaches were employed to derive LULC maps and analyze changes in the Yayo. Four years (1986, 1996, 2006 and 2017) were selected for LULC mapping. The selection of the years was dictated by the availability of cloud free satellite images from the earliest period possible (1986) and with the idea to map decade-level changes. In addition, the 31-year study period covers all the LULC dynamics due to major government regime changes (subsequently policy) and management approaches in the study area. LAND-SAT archives were useful as they captured longer-term land-use change dynamics compared to other data sources. Level II images were downloaded from USGS GLOVIS website (Http://glovis.usgs.gov). Due to the influence of the monsoon rainfall pattern in the study area, the availability of cloud and haze-free satellite data is best for acquisition dates between December and May. Moreover, LULC features are easier to distinguish from satellite images if they are captured during these dry season months. Thus, we gave preference to datasets acquired in these months of the year.
Ecosystem assessment requires both accurate classification and thematically detailed LULC maps. For this reason, a detailed second level classification scheme was chosen. First, a set of fairly broad thematic classes were identified and mapped (Level-I). Iteratively applying the classification process, each level I class was refined into subclasses (Level -II). Considering the complexity of the study area in terms of mixed land use/cover types and features, there was a need to conduct exploratory analysis before the final LULC classification efforts. Both unsupervised and supervised classification as well as object-oriented feature extraction approaches were used to produce the final LULC maps. Once the multi-temporal LULC maps were produced, change detection analysis was undertaken to understand the extent and spatial dynamics of changes between various LULC types. To estimate accuracy of the LULC maps, we have used 1700 reference points for each year, up to 100 points for one class depending on the area coverage and spatial distribution of each class over the study area, and overall accuracy value for each class in each period ranges between 86% and 92%. This accuracy assessment results fulfil accuracy value of a LULC map required for any application.

Quantifying and mapping ecosystem services
Ecosystem services are the various and diverse benefits that are provided to humans by nature, and can be modified by the humannature interactions. Generally, ESs are divided into four major categories: provisioning, regulating, supporting and cultural. Different services are grouped under these and an ideal ecosystem service accounting should consider all the elements under each category. However, time, resources and complexity do not allow assessment of all services under the four categories, thus selection of ESs can be considered based on the objectives at hand, resources available, and data availability. In this study, six ESs (i.e. total carbon (TC) sequestration, water yield or provision (WY), crop yield (CY), habitat provision quality (HQ), sediment retention (SR) and nutrient retention (NR)) were selected based on focus-group discussion with experts and local government representatives.
We used Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST version 3.5), developed by the Natural Capital Project (Nelson et al., 2009) to estimate the above ESs, except agricultural crop yield. While the InVEST model is coarse and conceptual, it provides useful information for resource management strategies and qualitative ranking of scenarios. Details on the theoretical formulations and practical applications in various parts of the world can be found in many papers Kareiva et al, 2011;Polasky et al., 2012;Tallis and Polasky, 2009). The summary of the models and data sources used are given in Table 1. More details on individual ES methods of estimation and data sources are provided in the subsequent sections. Below we present the key approaches and datasets used to derive the ESs for the different years of the Yayo.

Total carbon sequestration (TC)
Carbon storage is an important element of global climate regulation. The total carbon (TC) stored in an area depends on the storage in the aboveground biomass, belowground biomass, litter, deadwood and soil; which depends primarily on land use (e.g., forest, shrub land, grassland) and land management (Zheng et al., 2016). We used the InVEST Version 3.3.3 Carbon model to estimate the carbon storage status of a given LULC class (Jiang et al. 2017). The carbon storage is based on the longterm equilibrium conditions of a given LULC, and does not include temporal dynamics. The main input for InVEST carbon model is an empirical data for the four pools of carbon for each LULC class. If data is not available for a given pool, TC is estimated using the other pools. In this case, literature on carbon sequestration is used to estimate TC for the study area.

Water yields (WY)
The provision of freshwater is a key ES that depends on the healthy ecosystem functions of the landscape. The water yield model in InVEST is based on the Budyko hypothesis (Budyko, 1974), which is based on the simple water production function that is designed to be sensitive to the main driver of water resource changes such as land use and climate elements. While the objective of Budyko-like water yield analyses is not to estimate water quantity accurately, but to understand the relative and inter-watershed and sub-watershed comparisons. Abera et al (2019) used Budyko framework to understand how the climate and land-use change affects water yield in Ethiopia at the national scale.
The water yield on a given pixel x (WY (x)) is estimated as the difference between the annual precipitation (P(x)) and annual actual evapotranspiration (AET (x)), as follows: where the actual evapotranspiration component of the water balance is expressed by: where PET(x) is the potential evapotranspiration, ω(x)is a non-physical parameter that characterizes the climate-soil properties, ETo (x) is the reference evapotranspiration for pixel x, k c (ζ x ) is the plant evapotranspiration coefficient associated with the LULC ζ x on pixel x, AWC(x) is the volumetric plant available water content, and z is an empirical constant. AWC(x) is defined by the soil texture and effective rooting depth, and estimated as the product of the plant available water capacity (PAWC), the minimum of root restricting layer depth, and vegetation rooting depth (InVEST2.0 manual): Root restricting layer depth is the soil depth at which root penetration is inhibited because of physical or chemical characteristics. Vegetation rooting depth is often given as the depth at which 95% of a vegetation type's root biomass occurs. PAWC is the plant available water capacity, i.e. the difference between field capacity and wilting point (InVEST2.0 manual).
TerraClim precipitation and reference evapotranspiration reanalysis data (Abatzoglou et al., 2018) is used because it provides high resolution (4 km) and long term consistent data that can be used for change detection. The root restricting layer depth, plant available water content, maximum root depth is obtained from ISRIC SoilGrids dataset (Wu et al., 2008;Bao et al., 2016;Zheng et al., 2016). Evaporation coefficient for each LULC type is assigned based on Wu et al. (2008); Bao et al. (2016); Zheng et al. (2016). For both water yield and sediment retention, while it was recommended to use subwatershed scale results for easy comparison with observed data and calibration, we used pixel and administrative unit level analysis as it is more convenient to examine the spatial and temporal trend, tradeoffs and synergetic relationships between many ESs.

Habitat quality (HQ)
Habitat quality is used as an indicator for biodiversity. To assess the impacts of human activities on the biological diversity of the study area, we applied the habitat quality module of InVEST (v.2.4.4;Kareiva et al., 2011;Tallis et al., 2011). The model combines information on LULC suitability and threats to biodiversity to produce habitat quality maps. Shumi et al (2018) found that LULC is the important factor determining the vegetation diversity in southwest Ethiopia. The model is based on the hypothesis that areas with higher habitat quality support higher biodiversity and species richness.
The key inputs of habitat quality in the InVEST model are: 1) the suitability of each LU/LC type (H j ) for providing habitat for biodiversity, 2) anthropogenic threats that originates at pixel x (r x ) affecting habitat quality, and 3) the sensitivity of each LU/LC type to each threat. While the LULC maps for the four years are extracted as detailed in Section 2.2, the anthropogenic threats are also extracted from the LULC map but for specific threats. For instance, conversion to/expansion of cropland is one of the main drivers of biodiversity reduction, thus croplands are extracted as a source of threat as separate raster information. Likewise, other land uses such as urban and rural settlements are extracted for LULC maps as threats. In addition to the threats extracted from the landuse map, we have also spotted some point sources of habitat degradation such as a fertilizer company established in the study area, which is included in our analysis. The relative habitat suitability score (Hj), from 0 to 1, where 1 indicated the highest suitability to species has been assigned to LULC types. The last input of the model i.e. the sensitivity of habitat type to different threats, help to account the differentiated impacts of threats to different habitats. The impacts of the threats on the habitat is determined by: 1) the effect of the threat over space (i rxy ); 2) the relative weight of each threat's importance compared to the others (w r ); 3) the relative sensitivity of each habitat to each threat (S jr ). The stress level D xj of the grid x with land-use type j is calculated as follows: where y and Yr indicate all grid cells and the set of grid cells on r's raster map, respectively, θx is the level of accessibility in grid cell x. Finally, the habitat quality of each grid cell is determined by the cell's habitat suitability condition (H j ) and the degree of habitat degradation given in Eq. (9), as follows:  Table 3 Formulation used to derive the sediment retention services. EI 30 is rainfall erosivity of a single event; e r is the unit rainfall energy (MJ ha − 1 mm − 1 ); v r the rainfall volume (mm) during the r th time period of a rainfall event divided in k-parts; I 30 is the maximum 30-minutes rainfall intensity (mm h − 1 ); CI is connectivity index; SDRmax is the maximum theoretical SDR.
Term in eq (10) Formulation whereas References The results from the model are within the range of 0 to 1, with 1 representing the highest level of habitat quality. In other words, the impacts of the threat on habitat decreases as the distance from the degradation sources increases, threats with higher destructive values (in the scale of 0-1) has higher impacts and the more sensitive a habitat type is to a threat (higher S jr ), the more degraded the habitat type will be by the threat. In our study, we assigned the suitability of each LULC and the distance between the threats and the habitat based on literature obtained specifically from the Yayo forest area (Shumi et al., 2019;Ango et al., 2014;Mereta et al., 2013;Eshete, 2013;Woldegeorgis and Wube, 2012; Table 2). While the model is a relatively simplified representation of the issue, it is useful to understand the management approaches of core, buffer and transition zones and spatial effects of threats of habitats and which areas are most affected, and provide useful information for making a decision on initial assessment of conservation needs.

Agricultural crop yield (CY)
The InVEST Crop Production model is too coarse to capture spatial variability for agricultural yield, as the current climate data integrated in the model is based on 5 min resolution and does not capture the local biophysical variabilities such as soil fertility and agroecological zones. As a result, we used statistical data driven models, particularly random forest regression based on agronomic data collected at the national scale. The CIAT soil and agronomic dataset that is hosted by Ethiopia institute of agricultural research (EIAR) is used for predicting yield in different soil, topographic and agro-ecological conditions (see Kihara et al 2017). We developed a random forest (RF) model based on those spatially varying input data. Many covariates such as elevation, slope, and topographic position were derived from SRTM data (Rabus et al 2003), soil properties such as pH, texture, soil class, and soil organic carbon were derived from the ISRIC SoilGrids database (Hengl et al. 2017), soil moisture was derived from the TerraClim analysis dataset (Abatzoglou et al., 2018) and fertilizer application at national recommendation rate have been included as predictors. The analysis of agricultural production in this study is based on maize yield as it is the main crop of the area. The model is calibrated based on the 70% of the data points and validated for the remaining 30% of the data points. To build an optimal predictive model, we used R caret package (kuhn, 2008) to train the Random Forest model in ranger package (Wright and Ziegler, 2015). The RF model has R 2 = 0.53 and RMSE = 960 kg/ha of model performance indicating that the it can be used to provide yield estimation with good accuracy.

Sediment retention (SR)
The Sediment retention ES is estimated using the multiplication of the amount of annual soil loss on a pixel (according to RUSLE formulation) and a sediment delivery ratio, using the sediment delivery ratio (SDR) model of InVEST. Mathematically, sediment retention is estimated by (Hamel et al., 2015): where R is rainfall erosivity (MJ mm ha − 1 h − 1 y − 1 ), K is soil erodibility (ton⋅ha⋅hr (MJ⋅ha⋅mm) − 1 ), LS is slope length-gradient factor, C is cropmanagement factor, P-factor is mostly associated with 'management' and direct linkage with land use/cover can't provide accurate estimation, thus used P = 1 as the default value for all; and SDR is soil delivery ratio. The formulation and source of data for each term of Eq. (10) is given in Table 3.
The sediment retention model requires user-defined threshold for flow accumulation number that characterize the watershed drainage network; the IC0 and kb parameters are used the default values (i.e. threshold area = 1000 cells; IC0 = 0.5; kb = 2) (Sharp et al., 2014b). SR is avoided soil loss by the current land use compared to bare soil, weighted by the SDR factor. It is qualitative index as it does not account for retention from upstream sediment flowing through the given pixel. For both SR and NR estimations, the study area is based on upslope contributing area to each pixel, which basically a watershed enclosed at the lower part of Yayo biosphere reserve (Fig. A1). Once the analysis was done at watershed scale, the estimations were cropped to Yayo biosphere reserve area.

Nutrient retention (NR)
The InVEST nutrient delivery ratio (NDR) model is used to quantify the relative nutrient retention spatially under different land management scenarios (Redhead et al., 2018). The main inputs determining nutrient retention are LULC maps and a nutrient runoff proxy. We used the annual precipitation grid (TerraClim) as a proxy to capture the spatial variability in runoff potential (i.e. capacity to transport nutrient downstream) of the study area. In addition, the input data of the InVEST NDR model includes the biophysical table, threshold flow accumulation, Borselli k parameter, subsurface maximum retention efficiency of the nutrient, and subsurface critical length nutrient. Here, we focused on nitrogen as one example of a pollution indicator. The nutrient mass balance is based on: 1) nutrient sources associated with different LULC, and 2) the retention capacity that is related to LULC and geomorphology (particularly slope) along the flow path. The sources are related to nonpoint sources such as fertilizer application. The InVEST NDR model requires LULC maps, several parameter values for each distinct LULC class. These include the nutrient load applied to the land (kg ha − 1 y − 1 ), the proportional retention of that nutrient load, the length of flow path required to achieve that retention (in metres). Since there were no data on load factor, retention efficiency, and critical length of nitrogen in the Yayo area, the relevant parameters were assigned from the literature based on similar environments. Here, we assumed that the nutrient travels on the surface thus no need to assign the proportion of the nutrient load that travels via subsurface flow. Finally, the nutrient retention was estimates as a difference between total N load and N export amount for each grid.

Trade-off, synergy and clustering analysis
Pixel level tradeoffs and synergies among ES was conducted using the Pearson correlation coefficients (Vallet et al., 2018). Positive correlation means that there is a synergetic relation whereas negative correlation indicates a trade-off relationship between the two ESs. For static spatial correlation analysis, we used the correlation between pairs of ESs at the four years. For temporal correlation, we calculated the variations of each ESs between two consecutive dates (1986-1996, 1996-2005, 2005-2017) and between the start and end of the whole period studied (1986-2017).
In addition to pixel level individual ESs analysis, we aggregated the ESs estimation at the kebele level (the smallest administrative unit in Ethiopia). We prefer kebele to aggregate the ESs because the social process is an important factor that shaped the supply of and demand for ESs (Raudsepp-Hearne et al., 2010). In addition to aggregation of individual ES at kebele level, we have analyzed patterns of interaction among six ESs and identified a set of ESs that appeared together across kebele (i.e. bundles of ES). K-mean cluster analysis (MacQueen, 1967) was used to identify the distinct types of bundles of ESs and all 98  kebeles in the study areas. The k-means cluster algorithm was parameterized to build two to 20 clusters and we used the silhouette width index (Rousseeuw, 1987) to determine the optimal number of clusters. Two clusters were identified to be the optimal cluster, with the highest silhouette width value (Rousseeuw, 1987). The rose-wind diagram, based on dimensionless value was calculated based on normalized values for each ES, was used to visualize the ES bundles.

ES production possibility and efficiency frontiers
The set of combinations of ESs that can be supplied within a given landscape area is called production possibility. The boundary of the set is known as the production possibility frontiers or efficiency frontier (White et al., 2012). The production possibility set for each pair of ES is produced based on 98 Kebeles. We produced 6 by 6 (36) scatter plots based on aggregated ESs at kebele level. We then generated ESs bagplots around the scatterplot to further describe the ESs relationship in terms of dispersion, shapes, and direction of scatter plot graphically (Jopke et al., 2015). The envelope of the cloud points was computed using convex hull geometry, particularly alpha-shape (Edelsbrunner et al., 1983), and we used R package ggConvexHull for generating the envelope of cloud of points (Hijmans et al., 2017). The portion of the efficiency frontier of the envelope was then produced using the skyline function of the rPref package in R (Roocks and Roocks, 2019). The efficient frontier enables mapping pairs of ES in an efficient manner, so that an increase of service is possible without decreasing another. We developed a separate frontier for each pair of ES to inform the possibilities of synergies, tradeoffs and neutral combinations of ES supply. Synergy curve is oriented from lower left to upper right whereas the tradeoff efficiency curve is oriented from higher left to lower right (Jopke et al., 2015). If there is a tradeoff curve Fig. 3. Spatial distribution of six ES supply in Yayo in four years (1986,1996,2005,2017). Please note that the values are normalized between 0 to 1 for better visualization. Key: WYwater yield, SR-sediment retention, NR-Nitrogen retention, TC-total carbon, Yieldcrop yield, HQ -habitat quality.
between ES, then the community should know that when preferring one efficient ES combination over another, and show which combination of pairs of ES are in fact impossible.

Spatio-temporal trend of land use
Land-use/cover types for the years 1986, 1996, 2005, and 2017 are shown in Fig. 2. Spatially presenting LULC change using pixel-level multi-temporal maps makes the interpretation and visualization difficult as changes are happening at a very small scale (pixel). As a result, the multi-temporal LULC maps depicted in Fig. 2 appear to show no considerable change in the Yayo biosphere reserve, particularly in the core and buffer zones. Statistical analysis of the results however reveals substantial change in LULC in the Yayo over the last four decades (Table 4). For instance, cropland has been increasing over time apart from a short period of stability observed between 2005 and 2007. Conversely, the forest cover revealed a decreasing trend except for a small period of increase between 1996 and 2005. Similarly, grasslands also revealed an overall decreasing trend with a short period of increase between 1996 and 2005. Increment of the cultivated area in the landscape happened on the expense of forest and grassland areas in the landscape in the same periods. Generally, woodland and shrub/bush covered areas have experienced a decreasing trend. A considerable change of woodland was measured between 1996 and 2005 as this cover type decreased by about half. Both forest and cultivated landscape experienced a wavering character in terms of cover change, as sudden rise and fall is observed in between. According to the statistics presented in Table 4, the year 1996 represented the period when massive conversion of vegetation cover (both woody and non-woody) is detected and extreme increment of cultivated landscape. Given that this period represents episodes of political transitions, there was little or no law enforcement in controlling illegal human activities in forests. This is reflected in the statistics in 2005 when the government gave attention to natural resources protection and controlling of illegal land-use practices (deforestation and uncontrolled cropland expansion and settlements). Indeed, such generalizations need to be proved whether these changes are human driven or not and locate when and where they have happened. The former (driver of change) requires conducting zone level assessment defined by appropriate criteria. The latter issues (when and where) can also be verified by analyzing the type/direction of changes, which indicate the from-to changes processes taking place during a specific period of change assessment. The from-to change assessment that can indicate the time and the type as well as the location of the changes, was performed by comparing only the beginning and end of the analysis period (1986 vs. 2017).
The cumulative overall gain and loss of the major LULC classes between the 1986 and 2017 is depicted in Table A1. Out of the eight classes, woodland, forest, and grassland revealed a cumulative negative change, while the rest showed positive change. This does not mean that all of these classes have shown the same trends in the interim periods. Tables 4 and 1A show the changes between each period analyzed. According to the statistics presented in Table 4, the area lost is highest (− 105 km 2 ) in forests, followed by that in the grasslands (− 64 km 2 ). Overall (between 1986 and 2017), cultivated landscape revealed a considerable increment compared to the base year (247 km 2 ), or about 7% (132 km 2 ). This gain was largely from all kinds of vegetated landscapes. The loss of 6% forest cover, 4% of the grasslands and 3% of the woodlands over the last three decades contributed to the gain of the cropland (22%). Reid et al. (2000) found similar rapid LULC change and identified that the combined effect of drought, migration, change in settlement and land tenure policy are the main causes of change in southwestern Ethiopia. In western Ethiopia, forest cover is decreasing by 28% over the last 38 years due to agricultural land expansion (Betru et al., 2019).

Spatio-temporal patterns of individual ESs (Mapping ESs)
Generally, the ES supply map follows land-use/cover patterns (Fig. 3). For instance, regulating ES such as total carbon and habitat quality strongly follow the land-use pattern where highest values are observed in the protected forest located in the core zone and decreases to cultivated lands which are located in the transition zone (Fig. 3). Nutrient retention and sediment retention shows similar pattern but follow loose pattern of land-use maps. As expected crop yield is higher in the transition zone (Fig. 3). The highest values of water yield are observed in the northeastern part of the study area which corresponds to an open grassland and cropland mosaic. The temporal patterns of ES supply are presented in Figs. 3 and 4. The aggregate values of crop yield increased from 1986 to 2017. The increase in crop yield from 1986 to 2017 is due to an increase in cultivated lands in the transition zones (Fig. 4). Highest yield amount produced in the Yayo area was observed in 1996 due to the sudden increase in cropland (Fig. 4).
As the main inputs of the habitat quality model are the land-use and threat maps, which correspond to land uses that are related to agriculture and human settlement, it shows that the highest habitat quality is observed in the forest area (core zones) and decreases in the direction of the outer zones. The area of patches with high habitat qualities which were observed in 1986 have substantially decreased and concentrated to smaller patches of the main forest area (core zone) in 2017 (Fig. 4). Fig. 4 shows that the highest habitat quality is provided in the core forest followed by the buffer and transition zones, as would be expected, while the total pattern shows that the values of habitat quality have gradually declined between the years 1986-2017 (Fig. 4). This is consistent with the findings of Decuyper et al. (2018) where higher structural complexity (i.e. representing habitat complexity thus quality) was also found in the core forests and declined towards coffee forests, silvopasture, and plantation forests that were located in the transitional and buffer zones of the Kafa biosphere reserve (also located in SW Ethiopia). The decline in total habitat quality across the years is most likely due to the agriculture, and urban and rural settlement expansion into forest areas in the transition and buffer zones, as shown by the land-use maps (Fig. 2).
Our results show that the source of erosion is in the uplands of the study area and the highest retention service is observed around the main streams most likely due to the riparian vegetation. Nutrient leaching (nitrogen and phosphorus), the main cause of water pollution and environmental degradation, is high in the agricultural production area due to the application of fertilizers. In addition, animal manure in grazing lands can be sources of nitrate in downstream waters (Sakadevan and Nguyen, 2017). Here, the focus is on nitrogen retention service for land use.
InVest model results show that total carbon stored in forested ecosystems significantly declined from 33 Mg/pixel in 1986 to 25 Mg/pixel in 2017. The highest loss was observed in the transition zone followed by buffer zones. The carbon storage was almost constant in the core zone (Fig. 4). A similar study by Decuyper et al. (2018) in the Kafa biosphere reserve reported that higher above-ground biomass was found in the core zones of the biosphere reserve and declined towards plantation and silvopasture land uses that existed in the buffer and transition zones.
Land-use change has influenced the spatial distribution of water yield across the Yayo, as shown in Fig. 3, where water yield is higher in the northern part. Water yield increased from 1986 to 2017, and there was a particularly large increase from 630 mm in 1996 to 940 mm in 2017 for the whole study area. This increase in water yield is strongly related to conversion of forestland to cropland, as cropland has lower evapotranspiration than forest, which results in higher runoff. A recent meta-analysis study based on cases from all over the world reported that conversion of forest to non-forest land covers such as cultivation increased water yield, and vice versa (Filoso et al., 2017). This results however need to be taken carefully as the water yield in this analysis is based on all seasons, including rainy season, which is highly related to an increase in rainfall event runoff coefficients as the runoff generated in flooding events are not available for use by the community. Fig. 5a and b present spatial and temporal correlation between the pairs of ESs. Of the 36 possible pairs of ecosystem services, the highest positive correlation is observed between total carbon and habitat quality both for spatial and temporal correlations. Both total carbon and habitat quality are high in pristine and dense forest areas. Similar evidence that habitat quality and total carbon are higher in forest areas are reported in Decuyper et al. (2018). Wu et al., (2019) obtain the same highest correlation between habitat quality and total carbon in China. High positive correlation (synergies) are observed between habitat quality and nutrient retention, habitat quality and sediment retention, total carbon and sediment retention, and total carbon and nutrient retention (Fig. 5b). Generally, the positive correlation is observed between regulating services (Fig. 5).

Spatio-temporal tradeoff and synergy analysis among ESs
On the contrary, negative correlations (tradeoffs) are observed between crop yield and sediment retention, crop yield and nutrient retention, crop yield and habitat quality, as well as crop yield and total carbon. This is possibly because these pairs of ESs are associated with opposite land-use types. The two provisioning services show also negative correlation with the other ESs. However, these correlation coefficients particularly between crop yield and other ESs are very low (<0.08) (Fig. 5), suggesting that it is possible to move into synergy relationships where crop yield can be increased and other ESs can be enhanced concurrently with appropriate land management interventions.

ES bundles spatial configurations over time
A cluster analysis was performed on ES supply in the Yayo area, resulting in a partitioning of all 98 kebeles into two ES clusters (Fig. 6a). Based on the model estimations of 6 ESs, the ESs supply and their ecological features, the two ESs clusters are spatially segregated in the outer and inner part of the study area (Fig. 6a).
Cluster 1 is characterized by "high regulating services" i.e. high habitat quality, high total carbon, high sediment retention and nutrient retention, but with a lower crop yield and water yield services. Most of the areas classified into this category are the core zone and some patches in the central part of the biosphere reserve. It comprised about 39 Fig. 6. spatial clusters of ES supply determined by K-means method in Yayo: a) the spatial distribution of the clusters over time, b) ES profile of the clusters. Each slice of the pie chart represents an ES. Data was standardized to facilitate visualization. Key: WYwater yield, SR-sediment retention, NR-nitrogen retention, TC-total carbon, CYcrop yield, HQhabitat quality.
kebeles in 1986 and decreased to 36 kebeles in 2017. Cluster 2 is characterized as "high provision ES areas" i.e. high crop and water yield. This cluster has lower habitat quality, total carbon storage and sediment retention services. It is characterized by cultivated lands in the transition and buffer zones. The cluster increased from 56 kebeles in 1986 to 59 kebeles in 2017.
Generally, the cluster procedure yielded two groups and the spatial patterns were almost the same in all four years, with very few shifts of kebeles from one cluster to another (Fig. 6). For instance, one kebele in the core was grouped into cluster 1 in 1986 and shifted to cluster 2 in 2005 and 2017. There was also movement in the opposite direction with some kebeles in the southeastern part of the study changed to cluster 1 from cluster 2. Generally, the area for each cluster is consistent in the four-time steps with slight increase for ES bundle 1 over time. The clustering result in our study shows that all the six services can be nicely grouped into provisioning and regulating service zones. Yang et al. (2015) also divided 12 ESs into four groping where in one group constituted crop and livestock production separated from other regulating and cultural services. In fact, the spatial analysis of bundles of ESs made it possible to identify and locate which landscape is the dominant area for which ESs and guide land-use planning and management in different part of the world (Li et al., 2019;Dittrich et al., 2017;Turner et al., 2014).

Efficiency frontier of ES supply
The clouds of scatter points and production frontiers show diversified shapes and patterns of overall ES pairs (Fig. 7). Some pairs of ESs show dispersed clouds of points such as habitat quality and sediment retention whereas others show elongated and point of clouds such as habitat quality and total carbon; and nutrient retention and total carbon. The Pareto production frontier efficiency curve is also presented in red color in Fig. 7. The efficient frontier enables mapping pairs of ES in an efficient manner, so that an increase of service is possible without decreasing another. For most ES pairs, there are many options to achieve pareto efficiency ES productions. The length and shape of optimal combination (frontier line) determines the management decisions between the ES pairs. For most of the cases, e.g. crop production with other ESs, there are direct tradeoffs and provided long Pareto efficient combinations (Fig. 7). These pairs require management decisions where an increase in the provisioning of one service results in a proportional decrease of the other service, with no diminishing returns, and vice versa. This kind of efficiency frontier is expected as services tradeoff with each other. Some pairs of ESs show a single Pareto efficient Fig. 7. Results of ES production frontier. The dot represents the mean ES value at the kebele level (98 kebele in four years, which is about 392 observations in total). Pareto efficient production frontiers are shown in red and non-Pareto efficient portions are the cloud envelopes by the convex envelope. Key: WYwater yield, SRsediment retention, NR-Nitrogen retention, TC-total carbon, Yieldcrop yield, HQhabitat quality.
combination. For instance, between habitat quality and sediment retention; habitat quality and total carbon; habitat quality and nutrient retention, and total carbon and nutrient retention indicating that the two ESs can be managed independently as the gain in one service may not compromise the other service (Lester et al., 2013). In most of these cases, the patterns of point cloud show positive pattern suggesting landuse planning that facilitate an increase in one ecosystem service would have synergetic relationship with the others. The concave Pareto efficiency frontier line between total carbon and water yield indicates that there are scenarios that increase the short of one service substantially without a large cost to the other service. The shape for the production frontiers can also be affected by a few outliers. Generally, most observed combinations (the black dots) are far from the Pareto efficient combination(s), which confirms that landscape optimization for ecosystem service bundles rarely exist on the ground, and this is due to many reasons (Burgi et al., 2005;Nassauer, 1995;Schneeberger et al., 2007). Tallis and Polasky (2009) suggested that non-Pareto operation of land use could be related to societal choices. Land configuration that can increase ESs pairs to the efficiency frontier by identifying various scenarios of land restoration strategies and land-use management is a key for enhancing the current land use.

Caveats
Due to the lack of data, some of the model inputs and calibration parameters are derived from the official user's guide of the InVEST and literature, and it would be beneficial to improve the simulation results by using the information from the study area to calibrate the input parameters. Specifically, some of the challenges and that can be improved are: -The nitrogen retention model can be improved with the fertilizer application rate of the study area instead of using national average fertilizer application rates. -The literature reporting carbon stocks for various land uses is sparse and does not provide consistent carbon stocks for the various pools. A data collection campaign to measure precise carbon stock for all land uses/cover and carbon pools would improve the accuracy of total carbon estimations. -In addition to information regarding land use/cover related threats of habitat quality, data regarding the point-source threats distributed through the biosphere reserve such as illegal logging sites can improve the habitat quality model outputs. -Some of the datasets used for modelling the ES have uncertainty that can propagate in the modelling results. For instance, errors and uncertainties related to satellite rainfall, satellite evapotranspiration, and ISRIC soil information such as soil depth, water holding capacity can affect the water yield estimation. -The water yield model does not assess the contribution of land cover to water flows during dry season, which is a key regulating waterrelated ES provided by forests and other natural ecosystems. Also, this approach does not allow to estimate the proportion of the water yield that flows out of the pixel as runoff or as lateral flow. If runoff then it may constitute a disservice (instead of service) as it may contribute to pick flows during rainy season or to soil loss. If lateral flow, then may be interpreted as an ES as it is key to contribute with water during the dry season. A more complex approach for water yield that incorporated soil water movement and watershed-level flows is needed to understand better water yield as an ecosystem service. -This study is focused on biophysical and ES supply. Integration of these results with ESs demand from the local community and their preferences is important for developing land-use plans considering the ESs supply and preference of the local community.

Conclusions
We used the InVEST model to estimate ESs change from 1986 to 2017 in the UNESCO registered Yayo Coffee Forest Biosphere Reserve in southwest Ethiopia. Various tradeoffs, clustering, and optimal production front analysis methods were employed to understand the dynamics of ESs spatially and temporally. The production efficient frontier analysis shows that there are cases where an increase of one ES is possible without decreasing the other. The results revealed that there has been substantial change in land use/cover in the Yayo over the last four decades. Forest land has been decreasing from 1986 to 2017 whereas the cropland has shown an increasing trend and this is driving a negative overall trend in the provision of ESs. The highest conversion of forest to cultivated land is observed in 1996 in Yayo area mostly likely due to the political transition in the country and little/no law enforcement in controlling illegal human activities. Our analysis of ESs offers evidence that tradeoffs can be managed. Crop yield and on-site water yield show increasing trends, while total carbon and habitat quality shows decreasing trends, with the highest change is observed in the transition zone followed by buffer zones. We found high to moderate positive correlation (synergies) between key ESs like habitat quality, sediment retention, total carbon, and nutrient retention. We found only weakly negative correlations (tradeoffs) between crop yield and other ESs (habitat quality, total carbon, sediment and nutrient retention). This weak negative correlation (tradeoffs) between crop production and other ESs provision suggesting that it is possible, with appropriate land management interventions, to reverse the tradeoff relationships into synergy relationships where crop yield can be increased and other ESs can be enhanced Simultaneously. A strongly negative correlation between water yield and the other ESs, as intact forests with high habitat quality use more water than other land-cover types. This result has to be interpreted with cautious as this is not reflecting the role of forests in regulating water flows throughout different seasons and mitigating pick flows (and the threat of flooding) due to reduced runoff (compared to croplands). Although tradeoffs existed between increasing crop production and negative environmental effects such as carbon emission to atmosphere, erosion and water pollution, synergy can be achieved through application of crop management practices that can enhance carbon sequestration, sediment retention, nutrient retention Abegaz et al., 2020). It is then imperative to start facilitating the implementation of these practices to increase the synergies between crop production and other ES in the croplands. The study also outlined that there are many possibilities towards land use planning and crop land management options that has a potential to enhance food security while conserving biodiversity and maintaining environmental integrity in the region. This is in line with the national biodiversity conservation strategy and action plan that highlighted the value of forest and agricultural biodiversity and other ecosystem services as a core driver of economic growth and long-term food security and poverty alleviation.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table A1
Class composition based overall summary of change for the four periods considered for the analysis.