Ecosystem Services 50 (2021) 101338 Contents lists available at ScienceDirect Ecosystem Services journal homepage: www.elsevier.com/locate/ecoser Impacts of land use and land cover dynamics on ecosystem services in the Yayo coffee forest biosphere reserve, southwestern Ethiopia Wuletawu Abera a,*, Lulseged Tamene a, Tibebu Kassawmar b, Kalkidan Mulatu a, Habtemariam Kassa c, Louis Verchot d, Marcela Quintero d a International Center for Tropical Agriculture (CIAT), Addis Ababa, Ethiopia b Water and Land Resource Centre, Addis Ababa University, Addis Ababa, Ethiopia c Center for International Forestry Research (CIFOR), Ethiopia Office, Addis Ababa, Ethiopia d International Center for Tropical Agriculture (CIAT), Cali, Colombia A R T I C L E I N F O A B S T R A C T Keywords: Land management to increase food production while conserving the environment and associated ecosystem Ecosystem services services (ESs) is one of the major development and research challenges of the 21st Century. Any land-use practice Tradeoffs or change to obtain a particular ecosystem service affects the other ES positively or negatively. The dynamics of Synergies these changes is more marked in biodiversity hotspot areas like UNESCO registered Yayo coffee forest biosphere Pareto efficiency 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. 1. Introduction production and specific yields must be increased through intensifying farming systems with the application of inputs such as fertilizers and As the world population increases, demand for food will grow at an pesticides as it will not be possible to expand cultivated lands indefi- alarming rate, and meeting this food demand would require between nitely (Lambin and Meyfroidt, 2011). In addition, there is serious 320 and 850 million hectares of additional agricultural land by 2050 competition for other uses of land besides needs for expanding agri- (Fróna et al., 2019). In recent decades, the demand for growing agri- cultural areas that may even lead to conflicts. Expansion of agricultural cultural products has been partly met by increasing cultivated land land can lead to significant land-use and land-cover (LULC) changes (Boserup, 2017). In areas where population pressure outweighs resource (Tilman and Clark, 2014). Currently agriculture is responsible for 80% supply, there will be expansion of cultivation and grazing areas into of water use globally (FAO, 2017). In many parts of the world, the forests, marginal lands and other land uses to satisfy food demand expansion of urban and rural settlements has also shown high correla- (Fróna et al., 2019). However, in the future, the efficiency of agricultural tion with natural resource exploitation, habitat fragmentation and * Corresponding author at: Bole sub-city, Gurd Shola Area, P.O. Box 5689, Addis Ababa, Ethiopia. E-mail address: wuletawu.abera@cgiar.org (W. Abera). https://doi.org/10.1016/j.ecoser.2021.101338 Received 14 September 2020; Received in revised form 26 June 2021; Accepted 30 June 2021 2212-0416/© 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Fig. 1. Location of the Yayo coffee forest biosphere reserve in the southwestern Ethiopia (left) and the spatial distribution of land use and land cover (right). biodiversity loss (Lawler et al., 2014; Lyu et al., 2018). Agricultural gaps in that local communities have not thus far got sustained benefits intensification also causes negative impacts such as increased pollution (Getahun and Keno, 2019). As a result, there is conflict between the local of ground water and eutrophication of rivers and lakes (Sala et al., 2000; communities and the nearby biosphere reserve (Getahun and Keno, Matson et al., 1997) as well as loss in biodiversity. Severe environmental 2019). In addition, there is expansion of agricultural land, illegal cutting degradation due to LULC and unsustainable agricultural intensification of trees, and deforestation in Yayo, which are mainly prevalent in the can result in the loss of ecosystem services (ESs) that are essential to transition zones, and fears are growing that these could expand to the human wellbeing. Reconciling the two needs - increasing food produc- buffer and interior zones (Beyene, 2014; Getahun and Keno, 2019). tion and conserving the environment and associated ESs has thus Land-use and land-cover changes and management options that are become a major development and research challenge of the 21st Century intended to enhance a certain ES can change the supply of another ES (Foley et al., 2011; Brussaard et al., 2010). (Bennett et al., 2009). The emergent interaction can be desirable (syn- Ethiopia is emblematic of a conflict between securing food self- ergy) or undesirable (tradeoff). In the former case, an increase in one sufficiency through enhanced agricultural production and avoiding the ecosystem service is associated with an increasing supply of another loss ecosystem functions and ESs due to land degradation (Hurni et al., while the latter is when a gain in one ecosystem service results in a 2010). The country is hugely affected by land degradation, as 85% of the decline of the supply of another. It is thus essential to develop land-use total land is said to suffer from moderate to very serious levels of land practices and management plans that can enhance multiple ESs simul- degradation, costing about $4.3 billion per year (Gebreselassie et al., taneously (Balbi et al., 2015; Foley et al., 2011). In addition, it is 2016). Estimates of average soil losses range between 3.4 and 84.5 tons necessary to consider that trade-offs between the increasing needs of the ha-1 yr-1 with maximum rates reaching 300 tons ha-1 yr-1 (Hurni et al, local population and maintenance of ecosystem health are minimized. 2015; Gashaw, 2015). National level nutrient depletion rates were Balancing the provision of ESs with human development needs is thus estimated to be 122, 13 and 82 kg ha-1 yr-1 for N, P, and K, respectively necessary for sustainable conservation and development pathways. (Haileslassie et al., 2005). The soil loss and nutrient depletion result in Ecosystem service analyses have been increasingly used to analyze decline in agricultural productivity. Degradation of natural habitat such the impacts of environmental changes on societal and economic terms to as conversion of pristine forest to cropland and settlements is the pri- provide rational land management decisions. Spatially explicit assess- mary cause of biodiversity decline and loss (Fuller et al., 2007). ment of ESs is essential for informative land management. Thus, map- The Country has put various efforts in place to tackle land degra- ping of various ESs is required to incorporate ESs information into land dation and its associated impacts, including soil and water conservation use related policy development. Different approaches can be used to practices since the 1960s. Recently, encouraging results have been assess ESs at different scales (e.g., MIMES (Boumans et al., 2015); Co observed in terms of improvement of ESs due to land restoration (Abera $ting Nature (Mulligan, 2015); InVEST (Sharp et al., 2014a), SolVES et al., 2020; Adimassu et al., 2017). The country has also put different (Sherrouse et al., 2014); EVT (Briceno & Kochmer, 2014), TESSA (Peh guidelines and policy action to tackle land degradation and support land et al., 2013); ARIES (Villa et al., 2009), PA-BAT (Dudley and Stolton, restoration efforts. Among these is the declaration and registration of 2008)). Many studies have explored the intended use, flexibility, limi- biosphere reserves by UNESCO to protect two of the biodiversity hot- tations and uncertainty of these ecosystem service models (Sharp et al., spots of the world that Ethiopia hosts, namely: the Eastern Afromontane 2015; Hamel and Bryant, 2017; Dennedy-Frank et al., 2016; Redhead and the Horn of Africa hotspots. et al., 2016; Willcock et al., 2016; Malinga et al., 2015; Maes et al., The Yayo Coffee Forest Biosphere Reserve (hereafter referred to as 2012). In this study, we selected InVEST model due to its generaliz- Yayo) is located in southwestern parts of Ethiopia. The main objective of ability and publicity by scientific literature. ESs models are less used in establishing the biosphere reserve is to protect habitat characteristics Sub-Saharan Africa to guide decision making and policy implementation linked to biodiversity and ESs that are vital to the maintenance of human due to lack of required data and researchers capacity to run ESs models well-being (Gole, 2003). Though some studies show the encouraging (Willcock et al., 2016). More specifically, the applications of models for contributions of such interventions (Duguma et al., 2019), there are still ESs estimation are few in Ethiopia (Mengist and Soromessa, 2019). 2 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 The overall aim of this study is to estimate tradeoffs and synergies Table 1 between agricultural production and other ESs due to the LULC changes Carbon storage (t ha− 1) in the four carbon pools for each LULC class used in this happening in Yayo. This is useful to produce evidence and develop study and the source of the estimate. guidelines for optimal future investment and spatial targeting of agri- Land use/cover C_above C_below C_soil C_dead References cultural production, and environmental conservation and management High forest 243 45 163 0.03 Mohammed and efforts. This enables to explore land use options for enhancing food se- Bekele (2014) curity while conserving biodiversity and maintaining environmental Open degrading 151 51 111 10 Mohammed and integrity in Yayo. The specific objectives are: (1) to examine LULC forest Bekele (2014) change in the Yayo between 1986 to 2017, (2) to assess spatiotemporal Plantation 128 20 101 5 Tadesse et al. (2014), Mohammed and dynamics to identify hotspots and coldspots of selected ecosystem ser- Bekele (2014) vices, and (3) to understand tradeoffs and synergies between ecosystem Managed/ 65 40 25 6 Mohammed and services and develop production frontiers curve to identify optimal ES Coffee Forest Bekele (2014), provision. The study followed a conceptual modelling approach that Betemariyam et al. focuses on obtaining relative ES values, tradeoff and synergy analysis, (2020) Cropland with 1.82 0.0455 108 0 Mohammed and rather than making accurate prediction. trees Bekele (2014) Cropland 9.025 0.0455 50 0 Abegaz et al. (2020) 2. Materials and methods without trees Agroforestry 24 5 120 0 Abegaz et al. (2020); Betemariyam et al. 2.1. Study area (2020) Drained 15 35 74 4 Abegaz et al. (2020) The Yayo Coffee Forest Biosphere Reserve is located in the south- grassland western Ethiopia and has been designated as UNESCO Biosphere Undrained 15 35 74 4 Abegaz et al. (2020) grassland Reserve since 2010 (Fig. 1). The mean annual temperature is about Wetland-Marsh 2 2 7.5 2 Gebre (2018) 23.76 ◦C, and it is located in the rainiest areas of the country with annual Wetland- 2 2 7.5 2 Gebre (2018) precipitation of about 1625 mm/year (Mulatu and Getahun, 2018). The Swamp total area of the biosphere reserve is 1670 km2, covering about 2 Urban/Built- 5 5 15 2 InVEST manual administrative zones and 5 Woredas (third level administrative units in ups Rural 8 8 20 2 Ethiopia, equivalent to district). Yayo includes Eastern Afromontane settlement Biodiversity Hotspot, and montane rainforest fragments with wild Coffee Homestead in 11.739 10 64 6.1 Gebre (2018) arabica populations and encompasses important bird areas of interna- Rural tional significance. Following its designation as a UNESCO biosphere settlement Bare land 5 5 15 2 InVEST manual reserve, three management zones were created: core, buffer and tran- Landslide area 5 5 15 2 InVEST manual sition zones. The core zone (27,733 ha) is where there is no human Dense 30 30 54 13 Gebre (2018); InVEST intervention and is fully protected. The buffer zone (21,552 ha) is where woodland manual limited human intervention is possible only if compatible with conser- Dense shrub/ 30 30 54 13 Gebre (2018); InVEST vation objectives of the reserve. The transitional zone (117,736 ha) is bush manual found adjacent to the buffer zone and includes agricultural lands, wet- lands, grasslands, settlements, and forest fragments. In this area, mini- mal human activities such as cultivation around homesteads are between December and May. Moreover, LULC features are easier to allowed. Despite its importance in terms of biodiversity, forest coffee distinguish from satellite images if they are captured during these dry and ESs provision, the area is threatened by deforestation, land degra- season months. Thus, we gave preference to datasets acquired in these dation, and conflict between land uses and users (Dorresteijn et al., months of the year. 2017). Ecosystem assessment requires both accurate classification and thematically detailed LULC maps. For this reason, a detailed second 2.2. Land use, land cover, and change (LULC) mapping level classification scheme was chosen. First, a set of fairly broad the- matic classes were identified and mapped (Level-I). Iteratively applying Ecosystem service estimation and trade-off analysis require well the classification process, each level I class was refined into subclasses differentiated and spatially explicit LULC datasets. In this study, context- (Level –II). Considering the complexity of the study area in terms of based LULC mapping and change detection approaches were employed mixed land use/cover types and features, there was a need to conduct to derive LULC maps and analyze changes in the Yayo. Four years (1986, exploratory analysis before the final LULC classification efforts. Both 1996, 2006 and 2017) were selected for LULC mapping. The selection of unsupervised and supervised classification as well as object-oriented the years was dictated by the availability of cloud free satellite images feature extraction approaches were used to produce the final LULC from the earliest period possible (1986) and with the idea to map maps. Once the multi-temporal LULC maps were produced, change decade-level changes. In addition, the 31-year study period covers all detection analysis was undertaken to understand the extent and spatial the LULC dynamics due to major government regime changes (subse- dynamics of changes between various LULC types. To estimate accuracy quently policy) and management approaches in the study area. LAND- of the LULC maps, we have used 1700 reference points for each year, up SAT archives were useful as they captured longer-term land-use change to 100 points for one class depending on the area coverage and spatial dynamics compared to other data sources. Level II images were down- distribution of each class over the study area, and overall accuracy value loaded from USGS GLOVIS website (Http://glovis.usgs.gov). Due to the for each class in each period ranges between 86% and 92%. This accu- influence of the monsoon rainfall pattern in the study area, the avail- racy assessment results fulfil accuracy value of a LULC map required for ability of cloud and haze-free satellite data is best for acquisition dates any application. 3 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 2.3. Quantifying and mapping ecosystem services where the actual evapotranspiration component of the water balance is expressed by: Ecosystem services are the various and diverse benefits that are provided to humans by nature, and can be modified by the human- [ ( )w ]AET x PET x PET x (1/w)( ) ( ) ( ) nature interactions. Generally, ESs are divided into four major cate- = 1 + − 1 + (2) P(x) P(x) P(x) gories: provisioning, regulating, supporting and cultural. Different ser- vices are grouped under these and an ideal ecosystem service accounting PET(x) = kc(ζx).ETo(x) (3) should consider all the elements under each category. However, time, ( ) resources and complexity do not allow assessment of all services under ω AWC(x)(x) = Z + 1.25 (4) the four categories, thus selection of ESs can be considered based on the P(x) objectives at hand, resources available, and data availability. In this study, six ESs (i.e. total carbon (TC) sequestration, water yield or pro- where PET(x) is the potential evapotranspiration, ω(x)is a non-physical vision (WY), crop yield (CY), habitat provision quality (HQ), sediment parameter that characterizes the climate-soil properties, ETo (x) is the retention (SR) and nutrient retention (NR)) were selected based on reference evapotranspiration for pixel x, kc(ζx) is the plant evapotrans- focus-group discussion with experts and local government piration coefficient associated with the LULC ζxon pixel x, AWC(x) is the representatives. volumetric plant available water content, and z is an empirical constant. We used Integrated Valuation of Ecosystem Services and Tradeoffs AWC(x) is defined by the soil texture and effective rooting depth, and (InVEST version 3.5), developed by the Natural Capital Project (Nelson estimated as the product of the plant available water capacity (PAWC), et al., 2009) to estimate the above ESs, except agricultural crop yield. the minimum of root restricting layer depth, and vegetation rooting While the InVEST model is coarse and conceptual, it provides useful depth (InVEST2.0 manual): information for resource management strategies and qualitative ranking AWC(x) = min(Rest.layer.depth, root.depth)*PAWC (5) of scenarios. Details on the theoretical formulations and practical ap- plications in various parts of the world can be found in many papers Root restricting layer depth is the soil depth at which root penetra- (Daily et al., 2009; Kareiva et al, 2011; Polasky et al., 2012; Tallis and tion is inhibited because of physical or chemical characteristics. Vege- Polasky, 2009). The summary of the models and data sources used are tation rooting depth is often given as the depth at which 95% of a given in Table 1. More details on individual ES methods of estimation vegetation type’s root biomass occurs. PAWC is the plant available water and data sources are provided in the subsequent sections. Below we capacity, i.e. the difference between field capacity and wilting point present the key approaches and datasets used to derive the ESs for the (InVEST2.0 manual). different years of the Yayo. TerraClim precipitation and reference evapotranspiration reanalysis data (Abatzoglou et al., 2018) is used because it provides high resolution 2.3.1. Total carbon sequestration (TC) (4 km) and long term consistent data that can be used for change Carbon storage is an important element of global climate regulation. detection. The root restricting layer depth, plant available water con- The total carbon (TC) stored in an area depends on the storage in the tent, maximum root depth is obtained from ISRIC SoilGrids dataset (Wu aboveground biomass, belowground biomass, litter, deadwood and soil; et al., 2008; Bao et al., 2016; Zheng et al., 2016). Evaporation coefficient which depends primarily on land use (e.g., forest, shrub land, grassland) for each LULC type is assigned based on Wu et al. (2008); Bao et al. and land management (Zheng et al., 2016). We used the InVEST Version (2016); Zheng et al. (2016). For both water yield and sediment reten- 3.3.3 Carbon model to estimate the carbon storage status of a given tion, while it was recommended to use subwatershed scale results for LULC class (Jiang et al. 2017). The carbon storage is based on the long- easy comparison with observed data and calibration, we used pixel and term equilibrium conditions of a given LULC, and does not include administrative unit level analysis as it is more convenient to examine the temporal dynamics. The main input for InVEST carbon model is an spatial and temporal trend, tradeoffs and synergetic relationships be- empirical data for the four pools of carbon for each LULC class. If data is tween many ESs. 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 2.3.3. Habitat quality (HQ) the study area. Habitat quality is used as an indicator for biodiversity. To assess the impacts of human activities on the biological diversity of the study area, 2.3.2. Water yields (WY) we applied the habitat quality module of InVEST (v.2.4.4; Kareiva et al., The provision of freshwater is a key ES that depends on the healthy 2011; Tallis et al., 2011). The model combines information on LULC ecosystem functions of the landscape. The water yield model in InVEST suitability and threats to biodiversity to produce habitat quality maps. is based on the Budyko hypothesis (Budyko, 1974), which is based on Shumi et al (2018) found that LULC is the important factor determining the simple water production function that is designed to be sensitive to the vegetation diversity in southwest Ethiopia. The model is based on the main driver of water resource changes such as land use and climate the hypothesis that areas with higher habitat quality support higher elements. While the objective of Budyko-like water yield analyses is not biodiversity and species richness. to estimate water quantity accurately, but to understand the relative and The key inputs of habitat quality in the InVEST model are: 1) the inter-watershed and sub-watershed comparisons. Abera et al (2019) suitability of each LU/LC type (Hj) for providing habitat for biodiversity, used Budyko framework to understand how the climate and land-use 2) anthropogenic threats that originates at pixel x (rx) affecting habitat change affects water yield in Ethiopia at the national scale. quality, and 3) the sensitivity of each LU/LC type to each threat. While The water yield on a given pixel x (WY (x)) is estimated as the dif- the LULC maps for the four years are extracted as detailed in Section 2.2, ference between the annual precipitation (P(x)) and annual actual the anthropogenic threats are also extracted from the LULC map but for evapotranspiration (AET (x)), as follows: specific threats. For instance, conversion to/expansion of cropland is one of the main drivers of biodiversity reduction, thus croplands are WY(x) = P(x) − AET(x) (1) extracted as a source of threat as separate raster information. Likewise, other land uses such as urban and rural settlements are extracted for 4 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Table 2 Characteristics of threats to habitat quality considered to Yayo. Land use/cover Habitat Threats (r) considered as suitability score habitat [0–1] agriculture Rural settlement Urban settlement Plantation (coffee and Infrastructure (road, khat) industry, etc) Weight of Dmax Weight of Dmax Weight of Dmax Weight of Dmax Weight of Dmax threat [0–1] [km] threat [0–1] [km] threat [0–1] [km] threat [0–1] [km] threat [0–1] [km] High forest 1 0.8 0.5 1 1 0.9/1 2 0.7 0.3 0.9/1 1.5 Degraded forest 1 0.8 0.5 1 1 1 2 0.5 0.3 1 1.5 Plantation forest 0.6 0.6 0.5 0.7 1 1 2 0.4 0.3 1 1.5 Managed/coffee 0.6 0.4/0.6 0.5 1 1 1 2 0.2 0.3 1 1.5 forest Cropland with trees 0.6 0/0.3 0.5 0.5 1 0.7 2 0.2 0.3 1 1.5 Cropland without 0.5 0 0.5 0.2 1 0.6 2 0 0.3 1 1.5 trees Agroforestry 0.6 0.2/0.4 0.5 0.5 1 0.6 2 0.3 0.3 1 1.5 grasslands 0.7 0.6 0.5 0.7 1 0.9 2 0.4 0.3 1 1.5 Wetlands (marsh) 0.7 0.6/0.8 0.5 0.5 1 0.7 2 0.4 0.3 1 1.5 Wetlands (swamps) 0.5 0.6/0.8 0.5 0.5 1 0.7 2 0.4 0.3 1 1.5 Dense woodland 0.2 0.4/0.6 0.5 1 1 1 2 0.6 0.3 1 1.5 Open woodland 0.2 0.4/0.6 0.5 1 1 1 2 0.6 0.3 1 1.5 LULC maps as threats. In addition to the threats extracted from the land- ( ) R Yr use map, we have also spotted some point sources of habitat degradation ∑ ∑ wD rxj = ∑ rxirxyθxSjr (6) such as a fertilizer company established in the study area, which is Rr=1 y=1 r=1wr included in our analysis. The relative habitat suitability score (Hj), from ( ) 0 to 1, where 1 indicated the highest suitability to species has been dxy assigned to LULC types. The last input of the model i.e. the sensitivity of irxy = 1 − (7) dr max habitat type to different threats, help to account the differentiated im- pacts of threats to different habitats. The impacts of the threats on the where y and Yr indicate all grid cells and the set of grid cells on r’s raster habitat is determined by: 1) the effect of the threat over space (irxy); 2) map, respectively, θx is the level of accessibility in grid cell x. Finally, the relative weight of each threat’s importance compared to the others the habitat quality of each grid cell is determined by the cell’s habitat (wr); 3) the relative sensitivity of each habitat to each threat (Sjr). The suitability condition (Hj) and the degree of habitat degradation given in stress level Dxj of the grid x with land-use type j is calculated as follows: Eq. (9), as follows: Table 3 Formulation used to derive the sediment retention services. EI30 is rainfall erosivity of a single event; er is the unit rainfall energy (MJ ha− 1 mm− 1); vr the rainfall volume (mm) during the rth time period of a rainfall event divided in k-parts; I30 is the maximum 30-minutes rainfall intensity (mm h− 1); CI is connectivity index; SDRmax is the maximum theoretical SDR. Term Formulation whereas References in eq (10) ( ) R ∑n ∑mj (EI ∑30) EI kR j=1 k=1 k 30 = r=1erv I Panagos et al. r 30 = n [ (2017) e = 0.29 1 − 0.72e(− 0.05i ] r r ) K K = 0.1317*(0.2+ 0.3*e[ − 0.0256*SAN(1 − SIL/100)]*(SIL/(CLA + SIL 0) .3)*[1 − (0.25*TOC)/ SAN, SIL, and CLA are sand, silt and clay Yang et al. (2018) (TOC+ e((3.72 − 2.97*TOC)))]*[1 − (0.7*SN1)/(SN1+ e((22.9*SN1 − 5.51)))] weight content (%), respectively; TOC is soil organic carbon content (%); SN1 = 1 – SAN / 100 LS Developed from 90 m SRTM dem topographic data C LULC data based on review from local literature High forest = 0.01; Juliette (2008), Plantation forest = 0.05; Hurni (1985), Managed/Coffee Forest = 0.055; Kassawmar et al. Cropland with trees = 0.47; (2018) Cropland without trees = 0.45; Agroforestry/unidentified = 0.055; Drained grassland = 0.25; Undrained grassland = 0.21; Wetland = 0.001; Settlement = 0.4; Rural settlement homestead = 0.4; Bare land = 0.4; Landslide area = 0.9; Dense woodland = 0.03; Dense shrub/bush = 0.05 SDR SDRSDR max SDRmax is the maximum theoretical SDR = ( )IC − IC 1 exp 0 i+ k 5 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 ( ( ))) Dz The sediment retention model requires user-defined threshold for Qxj = Hj 1 xj − Dz (9) + kz flow accumulation number that characterize the watershed drainage xj network; the IC0 and kb parameters are used the default values (i.e. The results from the model are within the range of 0 to 1, with 1 threshold area = 1000 cells; IC0 = 0.5; kb = 2) (Sharp et al., 2014b). SR representing the highest level of habitat quality. In other words, the is avoided soil loss by the current land use compared to bare soil, impacts of the threat on habitat decreases as the distance from the weighted by the SDR factor. It is qualitative index as it does not account degradation sources increases, threats with higher destructive values (in for retention from upstream sediment flowing through the given pixel. the scale of 0–1) has higher impacts and the more sensitive a habitat type For both SR and NR estimations, the study area is based on upslope is to a threat (higher Sjr), the more degraded the habitat type will be by contributing area to each pixel, which basically a watershed enclosed at the threat. In our study, we assigned the suitability of each LULC and the the lower part of Yayo biosphere reserve (Fig. A1). Once the analysis was distance between the threats and the habitat based on literature ob- done at watershed scale, the estimations were cropped to Yayo tained specifically from the Yayo forest area (Shumi et al., 2019; Ango biosphere reserve area. et al., 2014; Mereta et al., 2013; Eshete, 2013; Woldegeorgis and Wube, 2012; Table 2). While the model is a relatively simplified representation 2.3.6. Nutrient retention (NR) of the issue, it is useful to understand the management approaches of The InVEST nutrient delivery ratio (NDR) model is used to quantify core, buffer and transition zones and spatial effects of threats of habitats the relative nutrient retention spatially under different land manage- and which areas are most affected, and provide useful information for ment scenarios (Redhead et al., 2018). The main inputs determining making a decision on initial assessment of conservation needs. nutrient retention are LULC maps and a nutrient runoff proxy. We used the annual precipitation grid (TerraClim) as a proxy to capture the 2.3.4. Agricultural crop yield (CY) spatial variability in runoff potential (i.e. capacity to transport nutrient The InVEST Crop Production model is too coarse to capture spatial downstream) of the study area. In addition, the input data of the InVEST variability for agricultural yield, as the current climate data integrated NDR model includes the biophysical table, threshold flow accumulation, in the model is based on 5 min resolution and does not capture the local Borselli k parameter, subsurface maximum retention efficiency of the biophysical variabilities such as soil fertility and agroecological zones. nutrient, and subsurface critical length nutrient. Here, we focused on As a result, we used statistical data driven models, particularly random nitrogen as one example of a pollution indicator. The nutrient mass forest regression based on agronomic data collected at the national balance is based on: 1) nutrient sources associated with different LULC, scale. The CIAT soil and agronomic dataset that is hosted by Ethiopia and 2) the retention capacity that is related to LULC and geomorphology institute of agricultural research (EIAR) is used for predicting yield in (particularly slope) along the flow path. The sources are related to non- different soil, topographic and agro-ecological conditions (see Kihara et point sources such as fertilizer application. The InVEST NDR model re- al 2017). We developed a random forest (RF) model based on those quires LULC maps, several parameter values for each distinct LULC class. spatially varying input data. Many covariates such as elevation, slope, These include the nutrient load applied to the land (kg ha− 1 y− 1), the and topographic position were derived from SRTM data (Rabus et al proportional retention of that nutrient load, the length of flow path 2003), soil properties such as pH, texture, soil class, and soil organic required to achieve that retention (in metres). Since there were no data carbon were derived from the ISRIC SoilGrids database (Hengl et al. on load factor, retention efficiency, and critical length of nitrogen in the 2017), soil moisture was derived from the TerraClim analysis dataset Yayo area, the relevant parameters were assigned from the literature (Abatzoglou et al., 2018) and fertilizer application at national recom- based on similar environments. Here, we assumed that the nutrient mendation rate have been included as predictors. The analysis of agri- travels on the surface thus no need to assign the proportion of the cultural production in this study is based on maize yield as it is the main nutrient load that travels via subsurface flow. Finally, the nutrient crop of the area. The model is calibrated based on the 70% of the data retention was estimates as a difference between total N load and N points and validated for the remaining 30% of the data points. To build export amount for each grid. 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 R2 0.53 and RMSE 960 kg/ha of model 2.4. Trade-off, synergy and clustering analysis = = performance indicating that the it can be used to provide yield estima- tion with good accuracy. Pixel level tradeoffs and synergies among ES was conducted using the Pearson correlation coefficients (Vallet et al., 2018). Positive cor- 2.3.5. Sediment retention (SR) relation means that there is a synergetic relation whereas negative The Sediment retention ES is estimated using the multiplication of correlation indicates a trade-off relationship between the two ESs. For the amount of annual soil loss on a pixel (according to RUSLE formu- static spatial correlation analysis, we used the correlation between pairs lation) and a sediment delivery ratio, using the sediment delivery ratio of ESs at the four years. For temporal correlation, we calculated the (SDR) model of InVEST. Mathematically, sediment retention is esti- variations of each ESs between two consecutive dates (1986–1996, mated by (Hamel et al., 2015): 1996–2005, 2005–2017) and between the start and end of the whole period studied (1986–2017). SR = R × K × LS(1 − CP) × SDR (10) In addition to pixel level individual ESs analysis, we aggregated the ESs estimation at the kebele level (the smallest administrative unit in where R is rainfall erosivity (MJ mm ha− 1 h− 1 y− 1), K is soil erodibility Ethiopia). We prefer kebele to aggregate the ESs because the social (ton⋅ha⋅hr (MJ⋅ha⋅mm)− 1), LS is slope length-gradient factor, C is crop- process is an important factor that shaped the supply of and demand for management factor, P-factor is mostly associated with ‘management’ ESs (Raudsepp-Hearne et al., 2010). In addition to aggregation of indi- and direct linkage with land use/cover can’t provide accurate estima- vidual ES at kebele level, we have analyzed patterns of interaction tion, thus used P = 1 as the default value for all; and SDR is soil delivery among six ESs and identified a set of ESs that appeared together across ratio. The formulation and source of data for each term of Eq. (10) is kebele (i.e. bundles of ES). K-mean cluster analysis (MacQueen, 1967) given in Table 3. was used to identify the distinct types of bundles of ESs and all 98 6 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Fig. 2. Second level classification scheme based produced LULC maps for the considered periods of analysis in 1986, 1996, 2005 and 2017. Table 4 Overall summary of change (beginning and end periods only) on major LULC types. The percentage change is calculated is in reference to the total study area. Classification classes 1986 2017 Area change between 1986 and 2017 Percent change between 1986 and 2017 Level 1 Level II km2 % km2 % km2 % Forest High forest 1022.1 61.2 545.1 32.6 − 477 − 28.6 Degraded/open forest 74.1 4.4 88.5 5.3 14.4 0.9 Plantation forest 0.4 0.0 1.1 0.1 0.7 0.1 Managed/coffee forest 8.3 0.5 364.5 21.8 356.2 21.3 Sub-total 1104.9 66.1 999.1 59.8 ¡105.8 ¡6.3 Cultivated land Cropland with trees 77.9 4.7 51.7 3.1 − 26.2 − 1.6 Cropland without trees 169.0 10.1 267.6 16.0 98.6 5.9 Agroforestry 0.2 0.0 60.4 3.6 60.2 3.6 Sub-total 247.1 14.8 379.7 22.7 132.6 7.9 Grassland Drained grasslands 83.1 5.0 89.6 5.4 6.5 0.4 undrained grasslands 74.4 4.5 0.0 0.0 − 74.4 − 4.5 Sub-total 157.6 9.4 89.6 5.4 ¡68 ¡4 Wetland Wetlands (marshes) 4.6 0.3 5.0 0.3 0.4 0 Wetlands(swamps) 4.9 0.3 6.3 0.4 1.4 0.1 Sub-total 9.5 0.6 11.2 0.7 1.7 0.1 Settlement Urban built-ups 2.7 0.2 7.5 0.4 4.8 0.2 Rural settlements 27.0 1.6 22.0 1.3 5 0.3 Homesteads 7.4 0.4 4.7 0.3 2.7 0.1 Sub-total 37.1 2.2 34.2 2.0 2.9 0.2 Barren land Bare land/surface 37.4 2.2 115.6 6.9 78.2 4.7 Sub-total 37.4 2.2 115.6 6.9 78.2 4.7 Woodland Dense woodland 16.1 1.0 26.9 1.6 10.8 0.6 Open woodland 60.8 3.6 13.9 0.8 − 46.9 − 2.8 Sub-total 76.9 4.6 40.8 2.4 − 36.1 − 2.2 Grand total 1670 100 1670 100 7 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 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: WY – water yield, SR-sediment retention, NR-Nitrogen retention, TC-total carbon, Yield – crop yield, HQ –habitat quality. kebeles in the study areas. The k-means cluster algorithm was parame- based on aggregated ESs at kebele level. We then generated ESs bagplots terized to build two to 20 clusters and we used the silhouette width around the scatterplot to further describe the ESs relationship in terms of index (Rousseeuw, 1987) to determine the optimal number of clusters. dispersion, shapes, and direction of scatter plot graphically (Jopke et al., Two clusters were identified to be the optimal cluster, with the highest 2015). The envelope of the cloud points was computed using convex hull silhouette width value (Rousseeuw, 1987). The rose-wind diagram, geometry, particularly alpha-shape (Edelsbrunner et al., 1983), and we based on dimensionless value was calculated based on normalized used R package ggConvexHull for generating the envelope of cloud of values for each ES, was used to visualize the ES bundles. points (Hijmans et al., 2017). The portion of the efficiency frontier of the envelope was then produced using the skyline function of the rPref 2.5. ES production possibility and efficiency frontiers 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 The set of combinations of ESs that can be supplied within a given is possible without decreasing another. We developed a separate frontier landscape area is called production possibility. The boundary of the set for each pair of ES to inform the possibilities of synergies, tradeoffs and is known as the production possibility frontiers or efficiency frontier neutral combinations of ES supply. Synergy curve is oriented from lower (White et al., 2012). The production possibility set for each pair of ES is left to upper right whereas the tradeoff efficiency curve is oriented from produced based on 98 Kebeles. We produced 6 by 6 (36) scatter plots higher left to lower right (Jopke et al., 2015). If there is a tradeoff curve 8 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Fig. 4. Time series dynamics of ES values in the three forest management units and average at the study area (Yayo). Key: WY – water yield, SR-sediment retention, NR-Nitrogen retention, TC-total carbon, Yield – crop yield, HQ – habitat quality. between ES, then the community should know that when preferring one enforcement in controlling illegal human activities in forests. This is efficient ES combination over another, and show which combination of reflected in the statistics in 2005 when the government gave attention to pairs of ES are in fact impossible. natural resources protection and controlling of illegal land-use practices (deforestation and uncontrolled cropland expansion and settlements). 3. Results and discussion Indeed, such generalizations need to be proved whether these changes are human driven or not and locate when and where they have 3.1. Spatio-temporal trend of land use happened. The former (driver of change) requires conducting zone level assessment defined by appropriate criteria. The latter issues (when and Land-use/cover types for the years 1986, 1996, 2005, and 2017 are where) can also be verified by analyzing the type/direction of changes, shown in Fig. 2. Spatially presenting LULC change using pixel-level which indicate the from-to changes processes taking place during a multi-temporal maps makes the interpretation and visualization diffi- specific period of change assessment. The from-to change assessment cult as changes are happening at a very small scale (pixel). As a result, that can indicate the time and the type as well as the location of the the multi-temporal LULC maps depicted in Fig. 2 appear to show no changes, was performed by comparing only the beginning and end of the considerable change in the Yayo biosphere reserve, particularly in the analysis period (1986 vs. 2017). core and buffer zones. Statistical analysis of the results however reveals The cumulative overall gain and loss of the major LULC classes be- substantial change in LULC in the Yayo over the last four decades tween the 1986 and 2017 is depicted in Table A1. Out of the eight (Table 4). For instance, cropland has been increasing over time apart classes, woodland, forest, and grassland revealed a cumulative negative from a short period of stability observed between 2005 and 2007. change, while the rest showed positive change. This does not mean that Conversely, the forest cover revealed a decreasing trend except for a all of these classes have shown the same trends in the interim periods. small period of increase between 1996 and 2005. Similarly, grasslands Tables 4 and 1A show the changes between each period analyzed. Ac- also revealed an overall decreasing trend with a short period of increase cording to the statistics presented in Table 4, the area lost is highest between 1996 and 2005. Increment of the cultivated area in the land- (− 105 km2) in forests, followed by that in the grasslands (− 64 km2). scape happened on the expense of forest and grassland areas in the Overall (between 1986 and 2017), cultivated landscape revealed a landscape in the same periods. Generally, woodland and shrub/bush considerable increment compared to the base year (247 km2), or about covered areas have experienced a decreasing trend. A considerable 7% (132 km2). This gain was largely from all kinds of vegetated land- change of woodland was measured between 1996 and 2005 as this cover scapes. The loss of 6% forest cover, 4% of the grasslands and 3% of the type decreased by about half. Both forest and cultivated landscape woodlands over the last three decades contributed to the gain of the experienced a wavering character in terms of cover change, as sudden cropland (22%). Reid et al. (2000) found similar rapid LULC change and rise and fall is observed in between. According to the statistics presented identified that the combined effect of drought, migration, change in in Table 4, the year 1996 represented the period when massive con- settlement and land tenure policy are the main causes of change in version of vegetation cover (both woody and non-woody) is detected southwestern Ethiopia. In western Ethiopia, forest cover is decreasing by and extreme increment of cultivated landscape. Given that this period 28% over the last 38 years due to agricultural land expansion (Betru represents episodes of political transitions, there was little or no law et al., 2019). 9 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Fig. 5. Pearson correlation between different ecosystem services: a) spatial correlation based on the 2017 ES maps, b) based on space and time combined corre- lations. Blue and red color shows positive and negative correlation, respectively. The correlations are calculated at pixel level. Key: WY – water yield, SR-sediment retention, NR-nitrogen retention, TC-total carbon, Yield – crop yield, HQ – habitat quality. 3.2. Spatio-temporal patterns of individual ESs (Mapping ESs) 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 Generally, the ES supply map follows land-use/cover patterns streams most likely due to the riparian vegetation. Nutrient leaching (Fig. 3). For instance, regulating ES such as total carbon and habitat (nitrogen and phosphorus), the main cause of water pollution and quality strongly follow the land-use pattern where highest values are environmental degradation, is high in the agricultural production area observed in the protected forest located in the core zone and decreases to due to the application of fertilizers. In addition, animal manure in cultivated lands which are located in the transition zone (Fig. 3). grazing lands can be sources of nitrate in downstream waters (Sakade- Nutrient retention and sediment retention shows similar pattern but van and Nguyen, 2017). Here, the focus is on nitrogen retention service follow loose pattern of land-use maps. As expected crop yield is higher in for land use. the transition zone (Fig. 3). The highest values of water yield are InVest model results show that total carbon stored in forested eco- observed in the northeastern part of the study area which corresponds to systems significantly declined from 33 Mg/pixel in 1986 to 25 Mg/pixel an open grassland and cropland mosaic. The temporal patterns of ES in 2017. The highest loss was observed in the transition zone followed by supply are presented in Figs. 3 and 4. The aggregate values of crop yield buffer zones. The carbon storage was almost constant in the core zone increased from 1986 to 2017. The increase in crop yield from 1986 to (Fig. 4). A similar study by Decuyper et al. (2018) in the Kafa biosphere 2017 is due to an increase in cultivated lands in the transition zones reserve reported that higher above-ground biomass was found in the (Fig. 4). Highest yield amount produced in the Yayo area was observed core zones of the biosphere reserve and declined towards plantation and in 1996 due to the sudden increase in cropland (Fig. 4). silvopasture land uses that existed in the buffer and transition zones. As the main inputs of the habitat quality model are the land-use and Land-use change has influenced the spatial distribution of water threat maps, which correspond to land uses that are related to agricul- yield across the Yayo, as shown in Fig. 3, where water yield is higher in ture and human settlement, it shows that the highest habitat quality is the northern part. Water yield increased from 1986 to 2017, and there observed in the forest area (core zones) and decreases in the direction of was a particularly large increase from 630 mm in 1996 to 940 mm in the outer zones. The area of patches with high habitat qualities which 2017 for the whole study area. This increase in water yield is strongly were observed in 1986 have substantially decreased and concentrated to related to conversion of forestland to cropland, as cropland has lower smaller patches of the main forest area (core zone) in 2017 (Fig. 4). evapotranspiration than forest, which results in higher runoff. A recent Fig. 4 shows that the highest habitat quality is provided in the core forest meta-analysis study based on cases from all over the world reported that followed by the buffer and transition zones, as would be expected, while conversion of forest to non-forest land covers such as cultivation the total pattern shows that the values of habitat quality have gradually increased water yield, and vice versa (Filoso et al., 2017). This results declined between the years 1986–2017 (Fig. 4). This is consistent with however need to be taken carefully as the water yield in this analysis is the findings of Decuyper et al. (2018) where higher structural based on all seasons, including rainy season, which is highly related to complexity (i.e. representing habitat complexity thus quality) was also an increase in rainfall event runoff coefficients as the runoff generated in found in the core forests and declined towards coffee forests, silvopas- flooding events are not available for use by the community. ture, and plantation forests that were located in the transitional and buffer zones of the Kafa biosphere reserve (also located in SW Ethiopia). 3.3. Spatio-temporal tradeoff and synergy analysis among ESs The decline in total habitat quality across the years is most likely due to the agriculture, and urban and rural settlement expansion into forest Fig. 5a and b present spatial and temporal correlation between the areas in the transition and buffer zones, as shown by the land-use maps pairs of ESs. Of the 36 possible pairs of ecosystem services, the highest (Fig. 2). positive correlation is observed between total carbon and habitat quality 10 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 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: WY – water yield, SR-sediment retention, NR-nitrogen retention, TC-total carbon, CY – crop yield, HQ – habitat quality. both for spatial and temporal correlations. Both total carbon and habitat (<0.08) (Fig. 5), suggesting that it is possible to move into synergy re- quality are high in pristine and dense forest areas. Similar evidence that lationships where crop yield can be increased and other ESs can be habitat quality and total carbon are higher in forest areas are reported in enhanced concurrently with appropriate land management Decuyper et al. (2018). Wu et al., (2019) obtain the same highest cor- interventions. relation between habitat quality and total carbon in China. High positive correlation (synergies) are observed between habitat quality and 3.4. ES bundles spatial configurations over time nutrient retention, habitat quality and sediment retention, total carbon and sediment retention, and total carbon and nutrient retention A cluster analysis was performed on ES supply in the Yayo area, (Fig. 5b). Generally, the positive correlation is observed between regu- resulting in a partitioning of all 98 kebeles into two ES clusters (Fig. 6a). lating services (Fig. 5). Based on the model estimations of 6 ESs, the ESs supply and their On the contrary, negative correlations (tradeoffs) are observed be- ecological features, the two ESs clusters are spatially segregated in the tween crop yield and sediment retention, crop yield and nutrient outer and inner part of the study area (Fig. 6a). retention, crop yield and habitat quality, as well as crop yield and total Cluster 1 is characterized by “high regulating services” i.e. high carbon. This is possibly because these pairs of ESs are associated with habitat quality, high total carbon, high sediment retention and nutrient opposite land-use types. The two provisioning services show also retention, but with a lower crop yield and water yield services. Most of negative correlation with the other ESs. However, these correlation the areas classified into this category are the core zone and some patches coefficients particularly between crop yield and other ESs are very low in the central part of the biosphere reserve. It comprised about 39 11 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 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: WY – water yield, SR- sediment retention, NR-Nitrogen retention, TC-total carbon, Yield – crop yield, HQ – habitat quality. kebeles in 1986 and decreased to 36 kebeles in 2017. Cluster 2 is et al., 2014). 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 3.5. Efficiency frontier of ES supply and buffer zones. The cluster increased from 56 kebeles in 1986 to 59 kebeles in 2017. The clouds of scatter points and production frontiers show diversi- Generally, the cluster procedure yielded two groups and the spatial fied shapes and patterns of overall ES pairs (Fig. 7). Some pairs of ESs patterns were almost the same in all four years, with very few shifts of show dispersed clouds of points such as habitat quality and sediment kebeles from one cluster to another (Fig. 6). For instance, one kebele in retention whereas others show elongated and point of clouds such as the core was grouped into cluster 1 in 1986 and shifted to cluster 2 in habitat quality and total carbon; and nutrient retention and total carbon. 2005 and 2017. There was also movement in the opposite direction with The Pareto production frontier efficiency curve is also presented in red some kebeles in the southeastern part of the study changed to cluster 1 color in Fig. 7. The efficient frontier enables mapping pairs of ES in an from cluster 2. Generally, the area for each cluster is consistent in the efficient manner, so that an increase of service is possible without four-time steps with slight increase for ES bundle 1 over time. The decreasing another. For most ES pairs, there are many options to achieve clustering result in our study shows that all the six services can be nicely pareto efficiency ES productions. The length and shape of optimal grouped into provisioning and regulating service zones. Yang et al. combination (frontier line) determines the management decisions be- (2015) also divided 12 ESs into four groping where in one group tween the ES pairs. For most of the cases, e.g. crop production with other constituted crop and livestock production separated from other regu- ESs, there are direct tradeoffs and provided long Pareto efficient com- lating and cultural services. In fact, the spatial analysis of bundles of ESs binations (Fig. 7). These pairs require management decisions where an made it possible to identify and locate which landscape is the dominant increase in the provisioning of one service results in a proportional area for which ESs and guide land-use planning and management in decrease of the other service, with no diminishing returns, and vice different part of the world (Li et al., 2019; Dittrich et al., 2017; Turner versa. This kind of efficiency frontier is expected as services tradeoff with each other. Some pairs of ESs show a single Pareto efficient 12 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 combination. For instance, between habitat quality and sediment 4. Conclusions retention; habitat quality and total carbon; habitat quality and nutrient retention, and total carbon and nutrient retention indicating that the We used the InVEST model to estimate ESs change from 1986 to 2017 two ESs can be managed independently as the gain in one service may in the UNESCO registered Yayo Coffee Forest Biosphere Reserve in not compromise the other service (Lester et al., 2013). In most of these southwest Ethiopia. Various tradeoffs, clustering, and optimal produc- cases, the patterns of point cloud show positive pattern suggesting land- tion front analysis methods were employed to understand the dynamics use planning that facilitate an increase in one ecosystem service would of ESs spatially and temporally. The production efficient frontier anal- have synergetic relationship with the others. The concave Pareto effi- ysis shows that there are cases where an increase of one ES is possible ciency frontier line between total carbon and water yield indicates that without decreasing the other. The results revealed that there has been there are scenarios that increase the short of one service substantially substantial change in land use/cover in the Yayo over the last four de- without a large cost to the other service. The shape for the production cades. Forest land has been decreasing from 1986 to 2017 whereas the frontiers can also be affected by a few outliers. Generally, most observed cropland has shown an increasing trend and this is driving a negative combinations (the black dots) are far from the Pareto efficient combi- overall trend in the provision of ESs. The highest conversion of forest to nation(s), which confirms that landscape optimization for ecosystem cultivated land is observed in 1996 in Yayo area mostly likely due to the service bundles rarely exist on the ground, and this is due to many political transition in the country and little/no law enforcement in reasons (Burgi et al., 2005; Nassauer, 1995; Schneeberger et al., 2007). controlling illegal human activities. Our analysis of ESs offers evidence Tallis and Polasky (2009) suggested that non-Pareto operation of land that tradeoffs can be managed. Crop yield and on-site water yield show use could be related to societal choices. Land configuration that can increasing trends, while total carbon and habitat quality shows increase ESs pairs to the efficiency frontier by identifying various sce- decreasing trends, with the highest change is observed in the transition narios of land restoration strategies and land-use management is a key zone followed by buffer zones. We found high to moderate positive for enhancing the current land use. correlation (synergies) between key ESs like habitat quality, sediment retention, total carbon, and nutrient retention. We found only weakly 3.6. Caveats negative correlations (tradeoffs) between crop yield and other ESs (habitat quality, total carbon, sediment and nutrient retention). This Due to the lack of data, some of the model inputs and calibration weak negative correlation (tradeoffs) between crop production and parameters are derived from the official user’s guide of the InVEST and other ESs provision suggesting that it is possible, with appropriate land literature, and it would be beneficial to improve the simulation results management interventions, to reverse the tradeoff relationships into by using the information from the study area to calibrate the input pa- synergy relationships where crop yield can be increased and other ESs rameters. Specifically, some of the challenges and that can be improved can be enhanced Simultaneously. A strongly negative correlation be- are: tween 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 - The nitrogen retention model can be improved with the fertilizer interpreted with cautious as this is not reflecting the role of forests in application rate of the study area instead of using national average regulating water flows throughout different seasons and mitigating pick fertilizer application rates. flows (and the threat of flooding) due to reduced runoff (compared to - The literature reporting carbon stocks for various land uses is sparse croplands). Although tradeoffs existed between increasing crop pro- and does not provide consistent carbon stocks for the various pools. A duction and negative environmental effects such as carbon emission to data collection campaign to measure precise carbon stock for all land atmosphere, erosion and water pollution, synergy can be achieved uses/cover and carbon pools would improve the accuracy of total through application of crop management practices that can enhance carbon estimations. carbon sequestration, sediment retention, nutrient retention (Abera - In addition to information regarding land use/cover related threats et al., 2020; Abegaz et al., 2020). It is then imperative to start facilitating of habitat quality, data regarding the point-source threats distributed the implementation of these practices to increase the synergies between through the biosphere reserve such as illegal logging sites can crop production and other ES in the croplands. The study also outlined improve the habitat quality model outputs. that there are many possibilities towards land use planning and crop - Some of the datasets used for modelling the ES have uncertainty that land management options that has a potential to enhance food security can propagate in the modelling results. For instance, errors and un- while conserving biodiversity and maintaining environmental integrity certainties related to satellite rainfall, satellite evapotranspiration, in the region. This is in line with the national biodiversity conservation and ISRIC soil information such as soil depth, water holding capacity strategy and action plan that highlighted the value of forest and agri- can affect the water yield estimation. cultural biodiversity and other ecosystem services as a core driver of - The water yield model does not assess the contribution of land cover economic growth and long-term food security and poverty alleviation. to water flows during dry season, which is a key regulating water- related ES provided by forests and other natural ecosystems. Also, Declaration of Competing Interest 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 The authors declare that they have no known competing financial then it may constitute a disservice (instead of service) as it may interests or personal relationships that could have appeared to influence contribute to pick flows during rainy season or to soil loss. If lateral the work reported in this paper. 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 Acknowledgements yield that incorporated soil water movement and watershed-level flows is needed to understand better water yield as an ecosystem This research was undertaken with financial support from the David service. and Lucile Packard Foundation, under the grant numbers 67,679, 2018. - This study is focused on biophysical and ES supply. Integration of The Water, Land and Ecosystems (WLE) programme of the CGIAR has these results with ESs demand from the local community and their also provided support to the research. The content is solely the re- preferences is important for developing land-use plans considering sponsibility of the authors and does not represent the official views of the ESs supply and preference of the local community. any of the donors. 13 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Appendix Table A1 Class composition based overall summary of change for the four periods considered for the analysis. Classification classes 1986 1996 2005 2017 Level 1 Level II km2 % km2 % km2 % km2 % Forest High forest 1022.1 61.2 722.6 43.2 724.2 43.3 545.1 32.6 Degraded/open forest 74.1 4.4 170.5 10.2 45.8 2.7 88.5 5.3 Plantation forest 0.4 0.0 0.2 0.0 1.1 0.1 1.1 0.1 Managed/coffee forest 8.3 0.5 26.2 1.6 247.0 14.8 364.5 21.8 Sub-total 1104.9 66.1 919.5 55.0 1018.1 60.9 999.1 59.8 Cultivated land Cropland with trees 77.9 4.7 105.2 6.3 74.6 4.5 51.7 3.1 Cropland without trees 169.0 10.1 400.2 24.0 227.1 13.6 267.6 16.0 Agroforestry 0.2 0.0 82.2 4.9 65.2 3.9 60.4 3.6 Sub-total 247.1 14.8 587.7 35.2 366.9 22.0 379.7 22.7 Grassland Drained grasslands 83.1 5.0 34.7 2.1 135.4 8.1 89.6 5.4 undrained grasslands 74.4 4.5 14.7 0.9 0.0 0.0 0.0 0.0 Sub-total 157.6 9.4 49.4 3.0 135.4 8.1 89.6 5.4 Wetland Wetlands (marshes) 4.6 0.3 2.5 0.1 6.8 0.4 5.0 0.3 Wetlands(swamps) 4.9 0.3 4.5 0.3 4.5 0.3 6.3 0.4 Sub-total 9.5 0.6 7.0 0.4 11.3 0.7 11.2 0.7 Settlement Urban built-ups 2.7 0.2 37.8 2.3 6.4 0.4 7.5 0.4 Rural settlements 27.0 1.6 0.6 0.0 27.1 1.6 22.0 1.3 Homesteads 7.4 0.4 36.9 2.2 7.9 0.5 4.7 0.3 Sub-total 37.1 2.2 75.2 4.5 41.4 2.5 34.2 2.0 Barren land Bare land/surface 37.4 2.2 13.2 0.8 67.8 4.1 115.6 6.9 Sub-total 37.4 2.2 13.2 0.8 67.8 4.1 115.6 6.9 Woodland Dense woodland 16.1 1.0 18.9 1.1 21.8 1.3 26.9 1.6 Open woodland 60.8 3.6 60.8 3.6 8.1 0.5 13.9 0.8 Sub-total 76.9 4.6 79.7 4.8 29.9 1.8 40.8 2.4 Grand total 1670 100 1731 103 1670 100 1670 100 Fig. A1. study area watershed to run the NDR and SDR models. Abera, W., Tamene, L., Abegaz, A., Solomon, D., 2019. Understanding climate and land surface changes impact on water resources using Budyko framework and remote sensing data in Ethiopia. J. Arid Environ. 167, 56–64. References Abera, W., Tamene, L., Tibebe, D., Adimassu, Z., Kassa, H., Hailu, H., Mekonnen, K., Desta, G., Sommer, R., Verchot, L., 2020. Characterizing and evaluating the impacts Abatzoglou, J.T., Dobrowski, S.Z., Parks, S.A., Hegewisch, K.C., 2018. TerraClimate, a of national land restoration initiatives on ecosystem services in Ethiopia. Land high-resolution global dataset of monthly climate and climatic water balance from Degrad. Dev. 31 (1), 37–52. 1958–2015. Sci. Data 5 (1). https://doi.org/10.1038/sdata.2017.191. Adimassu, Z., Langan, S., Johnston, R., Mekuria, W., Amede, T., 2017. Impacts of soil and Abegaz, A., Tamene, L., Abera, W., Yaekob, T., Hailu, H., Nyawira, S.S., Silva, M.D., water conservation practices on crop yield, run-off, soil loss and nutrient loss in Sommer, R., 2020. Soil organic carbon dynamics along chrono-sequence land-use Ethiopia: review and synthesis. Environ. Manage. 59 (1), 87–101. systems in the highlands of Ethiopia. Agric. Ecosyst. Environ. 300, 106997. https:// Ango, T.G., Börjeson, L., Senbeta, F., Hylander, K., 2014. Balancing ecosystem services doi.org/10.1016/j.agee.2020.106997. and disservices: smallholder farmers’ use and management of forest and trees in an agricultural landscape in southwestern Ethiopia. Ecol. Soc. 19 (1). 14 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Balbi, S., Prado, A.D., Gallejones, P., Geevan, C.P., Pardo, G., Pérez-Miñana, E., systems in Ethiopia using partial versus full nutrient balances. Agric. Ecosyst. Manrique, R., Hernandez-Santiago, C., Villa, F., 2015. Modeling trade-offs among Environ. 108 (1), 1–16. ecosystem services in agricultural production systems. Environ. Modell. Software 72, Hamel, P., Bryant, B.P., 2017. Uncertainty assessment in ecosystem services analyses: 314–326. seven challenges and practical responses. Ecosyst. Serv. 24, 1–15. Bao, Y., Li, T., Liu, H., Ma, T., Wang, H.X., Liu, K., Shen, X., Liu, X.H., 2016. Spatial and Hamel, P., Chaplin-Kramer, R., Sim, S., Mueller, C., 2015. A new approach to modeling temporal changes of water conservation of Loess Plateau in northern Shaanxi the sediment retention service (InVEST 3.0): Case study of the Cape Fear catchment, province by InVEST model. Geogr. Res. 35, 664–676. North Carolina, USA. Sci. Total Environ. 524-525, 166–177. Bennett, E.M., Peterson, G.D., Gordon, L.J., 2009. Understanding relationships among Hengl, T., Mendes de Jesus, J., Heuvelink, G.B.M., Ruiperez Gonzalez, M., Kilibarda, M., multiple ecosystem services. Ecol. Lett. 12 (12), 1394–1404. Blagotić, A., Shangguan, W., Wright, M.N., Geng, X., Bauer-Marschallinger, B., Betemariyam, M., Negash, M., Worku, A., 2020. Comparative analysis of carbon stocks in Guevara, M.A., Vargas, R., MacMillan, R.A., Batjes, N.H., Leenaars, J.G.B., home garden and adjacent coffee based agroforestry systems in Ethiopia. Small-Scale Ribeiro, E., Wheeler, I., Mantel, S., Kempen, B., Bond-Lamberty, B., 2017. For. SoilGrids250m: Global gridded soil information based on machine learning. PLoS Betru, T., Tolera, M., Sahle, K., Kassa, H., 2019. Trends and drivers of land use/land ONE 12 (2) e0169748. cover change in Western Ethiopia. Appl. Geogr. 104, 83–93. Hijmans, R.J., Phillips, S., Leathwick, J., Elith, J., Hijmans, M.R.J., 2017. Package Beyene, D.L., 2014. Assessing the impact of UNESCO biosphere reserves on forest cover ‘dismo’. Circles 9 (1), 1–68. change: The case of Yayu coffee forest biosphere reserve in Ethiopia. Unpublished Hurni, H., 1985. Erosion – Productivity – Conservation Systems in Ethiopia. Proceedings thesis. Wageningen, The Netherlands. IV International Conference on Soil Conservation pp. 654–674. Maracay, Venezuela. Boserup, E., 2017. The Conditions of Agricultural Growth: The Economics of Agrarian Hurni, H., Abate, S., Bantider, A., Debele, B., Ludi, E., Portner, B., Yitaferu, B. and Zeleke, Change under Population Pressure. Routledge, London, UK. G., 2010. Land degradation and sustainable land management in the highlands of Boumans, R., Roman, J., Altman, I., Kaufman, L., 2015. The Multiscale Integrated Model Ethiopia. of Ecosystem Services (MIMES): simulating the interactions of coupled human and Hurni, K., Zeleke, G., Kassie, M., Tegegne, B., Kassawmar, T., Teferi, E., Moges, A., natural systems. Ecosyst. Serv. 12, 30–41. Tadesse, D., Ahmed, M., Degu, Y., Kebebew, Z., 2015. Economics of land degradation Briceno, T., Kochmer, J., 2014. Ecosystem valuation toolkit: custom data pull for the (ELD) ethiopia case study. In: Soil Degradation and Sustainable Land Management in Chesapeake Bay and the Southern Appalachians regions. Earth Economics. the Rainfed Agricultural Areas of Ethiopia: An Assessment of the Economic Brussaard, L., Caron, P., Campbell, B., Lipper, L., Mainka, S., Rabbinge, R., Babin, D., Implications. Report for the Economics of Land Degradation Initiative, p. 94. Pulleman, M., 2010. Reconciling biodiversity conservation and food security: Jiang, W., Deng, Y., Tang, Z., Lei, X., Chen, Z., 2017. Modelling the potential impacts of scientific challenges for a new agriculture. Curr. Opin. Environ. Sustainability 2 (1- urban ecosystem changes on carbon storage under different scenarios by linking the 2), 34–42. CLUE-S and the InVEST models. Ecol. Model. 345, 30–40. Budyko, M.I., 1974. Climate and life. Academic press. Jopke, C., Kreyling, J., Maes, J., Koellner, T., 2015. Interactions among ecosystem Burgi, M., Hersperger, A.M., Schneeberger, N., 2005. Driving forces of landscape change- services across Europe: Bagplots and cumulative correlation coefficients reveal current and new directions. Landscape Ecol. 19 (8), 857–868. synergies, trade-offs, and regional patterns. Ecol. Ind. 49, 46–52. Daily, G.C., Polasky, S., Goldstein, J., Kareiva, P.M., Mooney, H.A., Pejchar, L., Juliette, K., 2008. Adaptation and validation of the universal soil loss equation (USLE) Ricketts, T.H., Salzman, J., Shallenberger, R., 2009. Ecosystem services in decision for the Ethiopian-Eritrean Highlands. Doctoral dissertation Verlag nicht ermittelbar. making: time to deliver. Front. Ecol. Environ. 7 (1), 21–28. Kareiva, P., Tallis, H., Ricketts, T.H., Daily, G.C., Polasky, S., 2011. Natural Capital: Decuyper, M., Mulatu, K.A., Brede, B., Calders, K., Armston, J., Rozendaal, D.M.A., Theory and Practice of Mapping Ecosystem Services: Theory and Practice of Mora, B., Clevers, J.G.P.W., Kooistra, L., Herold, M., Bongers, F., 2018. Assessing the Mapping Ecosystem Services. Oxford Press. structural differences between tropical forest types using terrestrial laser scanning. Kassawmar, T., Gessesse, G.D., Zeleke, G., Subhatu, A., 2018. Assessing the soil erosion For. Ecol. Manage. 429, 327–335. control efficiency of land management practices implemented through free Dennedy-Frank, P.J., Muenich, R.L., Chaubey, I., Ziv, G., 2016. Comparing two tools for community labor mobilization in Ethiopia. Int. Soil Water Conserv. Res. 6 (2), ecosystem service assessments regarding water resources decisions. J. Environ. 87–98. Manage. 177, 331–340. Kihara, J., Tibebe, D., Gurmensa, B., Desta, L., 2017. Towards understanding fertilizer Dittrich, A., Seppelt, R., Václavík, T., Cord, A.F., 2017. Integrating ecosystem service responses in Ethiopia, https://doi.org/10.7910/DVN/RKUMXB, Harvard Dataverse, bundles and socio-environmental conditions–a national scale analysis from V2. International Center for Tropical Agriculture (CIAT). Germany. Ecosyst. Serv. 28, 273–282. Kuhn, M., 2008. Building predictive models in R using the caret package. J. Stat. Softw. Dorresteijn, I., Schultner, J., Collier, N.F., Hylander, K., Senbeta, F., Fischer, J., 2017. 28 (5), 1–26. Disaggregating ecosystem services and disservices in the cultural landscapes of Lambin, E.F., Meyfroidt, P., 2011. Global land use change, economic globalization, and southwestern Ethiopia: a study of rural perceptions. Landscape Ecol. 32 (11), the looming land scarcity. Proc. Natl. Acad. Sci. 108 (9), 3465–3472. 2151–2165. Lawler, J.J., Lewis, D.J., Nelson, E., Plantinga, A.J., Polasky, S., Withey, J.C., Helmers, D. Dudley, N., Stolton, S., 2008. The Protected Areas Benefits Assessment Tool: A P., Martinuzzi, S., Pennington, D., Radeloff, V.C., 2014. Projected land-use change Methodology. WWF International. impacts on ecosystem services in the United States. Proc. Natl. Acad. Sci. 111 (20), Duguma, M.S., Feyssa, D.H., Biber-Freudenberger, L., 2019. Agricultural biodiversity and 7492–7497. ecosystem services of major farming systems: a case study in Yayo Coffee Forest Lester, S.E., Costello, C., Halpern, B.S., Gaines, S.D., White, C., Barth, J.A., 2013. Biosphere Reserve, Southwestern Ethiopia. Agriculture 9 (3), 48. Evaluating tradeoffs among ecosystem services to inform marine spatial planning. Edelsbrunner, H., Kirkpatrick, D., Seidel, R., 1983. On the shape of a set of points in the Mar. Policy 38, 80–89. plane. IEEE Trans. Inf. Theory 29 (4), 551–559. Li, T., Lü, Y., Fu, B., Hu, W., Comber, A.J., 2019. Bundling ecosystem services for Eshete, G.T., 2013. Biodiversity and livelihoods in southwestern Ethiopia: forest loss and detecting their interactions driven by large-scale vegetation restoration: enhanced prospects for conservation in shade coffee agroecosystems. Doctoral dissertation. UC services while depressed synergies. Ecol. Ind. 99, 332–342. Santa Cruz. Lyu, R., Zhang, J., Xu, M., Li, J., 2018. Impacts of urbanization on ecosystem services and FAO. 2017. Water for sustainable food and agriculture a report produced for the G20 their temporal relations: a case study in Northern Ningxia, China. Land Use Policy presidency of Germany. 77, 163–173. Filoso, S., Bezerra, M.O., Weiss, K.C.B., Palmer, M.A., Silva, L.C.R., 2017. Impacts of MacQueen, J., 1967. Some methods for classification and analysis of multivariate forest restoration on water yield: a systematic review. PLoS ONE 12 (8) e0183210. observations. No. 14. In: Proceedings of the fifth Berkeley symposium on Foley, J.A., Ramankutty, N., Brauman, K.A., Cassidy, E.S., Gerber, J.S., Johnston, M., mathematical statistics and probability, pp. 281–297. Mueller, N.D., O’Connell, C., Ray, D.K., West, P.C., Balzer, C., 2011. Solutions for a Maes, J., Egoh, B., Willemen, L., Liquete, C., Vihervaara, P., Schägner, J.P., Grizzetti, B., cultivated planet. Nature 478 (7369), 337–342. Drakou, E.G., La Notte, A., Zulian, G., Bouraoui, F., 2012. Mapping ecosystem Fróna, D., Szenderák, J., Harangi-Rákos, M., 2019. The Challenge of Feeding the World. services for policy support and decision making in the European Union. Ecosyst. Sustainability 11 (20), 5816. Serv. 1 (1), 31–39. Fuller, R.A., Irvine, K.N., Devine-Wright, P., Warren, P.H., Gaston, K.J., 2007. Malinga, R., Gordon, L.J., Jewitt, G., Lindborg, R., 2015. Mapping ecosystem services Psychological benefits of greenspace increase with biodiversity. Biol. Lett. 3 (4), across scales and continents–A review. Ecosyst. Serv. 13, 57–63. 390–394. Matson, P.A., Parton, W.J., Power, A.G., Swift, M.J., 1997. Agricultural intensification Gashaw, T., 2015. The implications of watershed management for reversing land and ecosystem properties. Science 277 (5325), 504–509. degradation in Ethiopia. Res. J. Agric. Environ. Manag. 4 (1), 5–12. Mengist, W., Soromessa, T., 2019. Assessment of forest ecosystem service research trends Gebre, A.B., 2018. Comparison of bulk density methods in determining soil organic and methodological approaches at global level: a meta-analysis. Environ. Syst. Res. 8 carbon storage under different land use types. J. Soil Sci. Environ. Manag. 9 (1), (1), 1–18. 13–20. Mereta, S.T., Boets, P., De Meester, L., Goethals, P.L.M., 2013. Development of a Gebreselassie, S., Kirui, O.K., Mirzabaev, A., 2016. Economics of land degradation and multimetric index based on benthic macroinvertebrates for the assessment of natural improvement in Ethiopia. In: Economics of Land Degradation and Improvement–A wetlands in Southwest Ethiopia. Ecol. Ind. 29, 510–521. Global Assessment for Sustainable Development. Springer, Cham, pp. 401–430. Mohammed, A., Bekele, L., 2014. Changes in carbon stocks and sequestration potential Getahun, D., Keno, E.T., 2019. Attitudes and perceptions of the local community towards under native forest and adjacent land use systems at Gera, South-Western Ethiopia. Yayo Coffee Forest Biosphere Reserve, Ilu Abba Bora Zone of Oromia National Global J. Sci. Front. Res. Agric. Vet. 14 (10). Regional State. Ethiopian J. Sci. Sustainable Dev. 6 (1), 79–90. Mulligan, M., 2015. Trading off agriculture with nature’s other benefits, spatially. In: Gole, T.W., 2003. Vegetation of the Yayu Forest in SW Ethiopia: Impacts of Human Use Impact of Climate Change on Water Resources in Agriculture. CRC Press, and Implications for In situ Conservation of Wild Coffea arabica L populations, . pp. 184–204. Cuvillier Verlag. Mulatu, T., Getahun, A., 2018. Diversity of anurans in forest fragments of southwestern Haileslassie, A., Priess, J., Veldkamp, E., Teketay, D., Lesschen, J.P., 2005. Assessment of Ethiopia: The case of the Yayu Coffee Forest Biosphere Reserve (YCFBR). Amphibian soil nutrient depletion and its spatial variability on smallholders’ mixed farming Reptile Conserv. 12 (2), 30–40. 15 W. Abera et al. E c o s y s t e m S e r v i c es 50 (2021) 101338 Nassauer, J.I., 1995. Culture and changing landscape structure. Landscape Ecol. 10 (4), Sherrouse, B.C., Semmens, D.J., Clement, J.M., 2014. An application of Social Values for 229–237. Ecosystem Services (SolVES) to three national forests in Colorado and Wyoming. Nelson, E., Mendoza, G., Regetz, J., Polasky, S., Tallis, H., Cameron, D., Chan, K.M., Ecol. Ind. 36, 68–79. Daily, G.C., Goldstein, J., Kareiva, P.M., Lonsdorf, E., 2009. Modeling multiple Shumi, G., Rodrigues, P., Schultner, J., Dorresteijn, I., Hanspach, J., Hylander, K., ecosystem services, biodiversity conservation, commodity production, and tradeoffs Senbeta, F., Fischer, J., 2019. Conservation value of moist evergreen Afromontane at landscape scales. Front. Ecol. Environ. 7 (1), 4–11. forest sites with different management and history in southwestern Ethiopia. Biol. Panagos, P., Borrelli, P., Meusburger, K., Yu, B., Klik, A., Jae Lim, K., Yang, J.E., Ni, J., Conserv. 232, 117–126. Miao, C., Chattopadhyay, N., Sadeghi, S.H., Hazbavi, Z., Zabihi, M., Larionov, G.A., Shumi, G., Schultner, J., Dorresteijn, I., Rodrigues, P., Hanspach, J., Hylander, K., Krasnov, S.F., Gorobets, A.V., Levi, Y., Erpul, G., Birkel, C., Hoyos, N., Naipal, V., Senbeta, F., Fischer, J., Essl, F., 2018. Land use legacy effects on woody vegetation in Oliveira, P.T.S., Bonilla, C.A., Meddi, M., Nel, W., Al Dashti, H., Boni, M., agricultural landscapes of south-western Ethiopia. Divers. Distrib. 24 (8), Diodato, N., Van Oost, K., Nearing, M., Ballabio, C., 2017. Global rainfall erosivity 1136–1148. assessment based on high-temporal resolution rainfall records. Sci. Rep. 7 (1) Tadesse, G., Zavaleta, E., Shennan, C., 2014. Effects of land-use changes on woody https://doi.org/10.1038/s41598-017-04282-8. species distribution and above-ground carbon storage of forest-coffee systems. Agric. Peh, K.-H., Balmford, A., Bradbury, R.B., Brown, C., Butchart, S.H.M., Hughes, F.M.R., Ecosyst. Environ. 197, 21–30. Stattersfield, A., Thomas, D.H.L., Walpole, M., Bayliss, J., Gowing, D., Jones, J.P.G., Tallis, H., Polasky, S., 2009. Mapping and valuing ecosystem services as an approach for Lewis, S.L., Mulligan, M., Pandeya, B., Stratford, C., Thompson, J.R., Turner, K., conservation and natural-resource management. Ann. N. Y. Acad. Sci. 1162 (1), Vira, B., Willcock, S., Birch, J.C., 2013. TESSA: A toolkit for rapid assessment of 265–283. ecosystem services at sites of biodiversity conservation importance. Ecosyst. Serv. 5, Tallis, H.T., Ricketts, T., Guerry, A.D., Nelson, E., Ennaanay, D., Wolny, S., Olwero, N., 51–57. Vigerstol, K., Pennington, D., Mendoza, G. and Aukema, J., 2011. InVEST 2.1 beta Polasky, S., Johnson, K., Keeler, B., Kovacs, K., Nelson, E., Pennington, D., Plantinga, A. user’s guide. the natural capital project. J., Withey, J., 2012. Are investments to promote biodiversity conservation and Tilman, D., Clark, M., 2014. Global diets link environmental sustainability and human ecosystem services aligned? Oxford Rev. Econ. Policy 28 (1), 139–163. health. Nature 515 (7528), 518–522. Rabus, B., Eineder, M., Roth, A., Bamler, R., 2003. The shuttle radar topography Turner, K.G., Odgaard, M.V., Bøcher, P.K., Dalgaard, T., Svenning, J.C., 2014. Bundling mission—a new class of digital elevation models acquired by spaceborne radar. ecosystem services in Denmark: trade-offs and synergies in a cultural landscape. ISPRS J. Photogramm. Remote Sens. 57 (4), 241–262. Landscape Urban Plann. 125, 89–104. Raudsepp-Hearne, C., Peterson, G.D., Bennett, E.M., 2010. Ecosystem service bundles for Vallet, A., Locatelli, B., Levrel, H., Wunder, S., Seppelt, R., Scholes, R.J., Oszwald, J., analyzing tradeoffs in diverse landscapes. Proc. Natl. Acad. Sci. 107 (11), 2018. Relationships between ecosystem services: comparing methods for assessing 5242–5247. tradeoffs and synergies. Ecol. Econ. 150, 96–106. Redhead, J.W., May, L., Oliver, T.H., Hamel, P., Sharp, R., Bullock, J.M., 2018. National Villa, F., Ceroni, M., Bagstad, K., Johnson, G., Krivov, S., 2009. September. ARIES scale evaluation of the InVEST nutrient retention model in the United Kingdom. Sci. (Artificial Intelligence for Ecosystem Services): a new tool for ecosystem services Total Environ. 610, 666–677. assessment, planning, and valuation. In: Proceedings of the 11th Annual BIOECON Redhead, J.W., Stratford, C., Sharps, K., Jones, L., Ziv, G., Clarke, D., Oliver, T.H., Conference on Economic Instruments to Enhance the Conservation and Sustainable Bullock, J.M., 2016. Empirical validation of the InVEST water yield ecosystem Use of Biodiversity, pp. 21–22. service model at a national scale. Sci. Total Environ. 569, 1418–1426. White, C., Halpern, B.S., Kappel, C.V., 2012. Ecosystem service tradeoff analysis reveals Reid, R.S., Kruska, R.L., Muthui, N., Taye, A., Wotton, S., Wilson, C.J., Mulatu, W., 2000. the value of marine spatial planning for multiple ocean uses. Proc. Natl. Acad. Sci. Land-use and land-cover dynamics in response to changes in climatic, biological and 109 (12), 4696–4701. socio-political forces: the case of southwestern Ethiopia. Landscape Ecol. 15 (4), Willcock, S., Hooftman, D., Sitas, N., O’Farrell, P., Hudson, M.D., Reyers, B., 339–355. Eigenbrod, F., Bullock, J.M., 2016. Do ecosystem service maps and models meet Roocks, P., & Roocks, M. P. 2019. Package ‘rPref’. stakeholders’ needs? A preliminary survey across sub-Saharan Africa. Ecosyst. Serv. Rousseeuw, P.J., 1987. Silhouettes: a graphical aid to the interpretation and validation of 18, 110–117. cluster analysis. J. Comput. Appl. Math. 20, 53–65. Woldegeorgis, G., Wube, T., 2012. A survey on mammals of the Yayu forest in Southwest Sakadevan, K., Nguyen, M.L., 2017. Livestock production and its impact on nutrient Ethiopia. SINET: Ethiopian J. Sci. 35 (2), 135–138. pollution and greenhouse gas emissions. In: Advances in agronomy. Academic Press, Wright, M.N., Ziegler, A., 2015. ranger: A fast implementation of random forests for high pp. 147–184. dimensional data in C++ and R. arXiv preprint arXiv:1508.04409. Sala, O.E., Chapin, F.S., Armesto, J.J., Berlow, E., Bloomfield, J., Dirzo, R., Huber- Wu, B.F., Xiong, J., Yan, N.N., Yang, L.D., Du, X., 2008. ETWatch for monitoring regional Sanwald, E., Huenneke, L.F., Jackson, R.B., Kinzig, A., Leemans, R., 2000. Global evapotranspiration with remote sensing. Adv. Water Sci. 19 (5), 671–678. biodiversity scenarios for the year 2100. Science 287 (5459), 1770–1774. Wu, Y., Tao, Y., Yang, G., Ou, W., Pueppke, S., Sun, X., Chen, G., Tao, Q., 2019. Impact of Schneeberger, N., Bürgi, M., Hersperger, A.M., Ewald, K.C., 2007. Driving forces and land use change on multiple ecosystem services in the rapidly urbanizing Kunshan rates of landscape change as a promising combination for landscape change City of China: past trajectories and future projections. Land Use Policy 85, 419–427. research—an application on the northern fringe of the Swiss Alps. Land Use Policy Yang, W., Dietz, T., Kramer, D.B., Ouyang, Z., Liu, J., 2015. An integrated approach to 24 (2), 349–361. understanding the linkages between ecosystem services and human well-being. Sharp, R., Ricketts, T., Guerry, A.D., 2015. InVEST 3.2. 0 Beta User’s Guide. The Natural Ecosyst. Health Sustainability 1 (5), 1–12. Capital Project. Stanford University, Stanford, CA, USA. Yang, Y., Zhao, R., Shi, Z., Rossel, R.A.V., Wan, D., Liang, Z., 2018. Integrating multi- Sharp, R., Tallis, H.T., Ricketts, T., Guerry, A.D., Wood, S.A., Chaplin-Kramer, R., Nelson, source data to improve water erosion mapping in Tibet, China. Catena 169, 31–45. E., Ennaanay, D., Wolny, S., Olwero, N., Vigerstol, K., 2014a. InVEST user’s guide. Zheng, H., Li, Y., Robinson, B.E., Liu, G., Ma, D., Wang, F., Lu, F., Ouyang, Z., Daily, G.C., The natural capital project. The Natural Capital. 2016. Using ecosystem service trade-offs to inform water conservation policies and Sharp, R., Tallis, H.T., Ricketts, T., Guerry, A.D., Wood, S.A., Chaplin-Kramer, R., Nelson, management practices. Front. Ecol. Environ. 14 (10), 527–532. E., Ennaanay, D., Wolny, S., Olwero, N., Vigerstol, K., 2014b. InVEST user’s guide. The Natural Capital Project: Stanford, CA, USA. 16