Impacts of large-scale landscape restoration on spatio-temporal dynamics of ecosystem services in the Chinese Loess Plateau

From 1999 onwards, China has initiated a large-scale landscape restoration project in the Chinese Loess Plateau, which has had profound but variable impacts on the local ecosystem services supply. The dynamics of ecosystem services throughout the restoration process remain poorly understood. To analyze the spatial and temporal dynamics in ecosystem services before and after the implementation of the land restoration project, and to understand trade-offs and synergies between multiple ecosystem services. We used the InVEST model and statistical yearbook data to quantify the ecosystem services over the period 1990–2018 for the Yan’an area and applied the concept of ecosystem service bundles to understand the dynamics of 11 ecosystem services over its 13 constituent counties. A significant increase of fruit production, sediment retention, habitat quality, aesthetic landscape value, and learning and inspiration value was found over time in the Yan’an area, while a decrease of timber production and water yield was also observed. The majority of the county-level ecosystem service bundles were transformed from having a focus on timber production to aesthetic landscape value. The dynamics of ecosystem services change induced by land restoration was discovered to start with increasing regulating services at the expense of provisioning services, while cultural services exceeded regulating services and occupied the main proportion subsequently. Both trade-offs and synergies were found between provisioning, regulating and cultural services. Implementation of the large-scale restoration project is recognized as a key driving force inducing change of ecosystem services, starting with an improvement of regulating services followed by a gradually evolving prominence of cultural services.


Introduction
Over 40% of the world's land surface are arid and semi-arid areas, which are ecologically vulnerable, sensitive to erosion and facing deterioration risks (Allan et al. 2013). In order to manage and address land degradation and ecological deterioration issues, a number of large-scale land restoration programs have been implemented worldwide, which have significantly altered land use and ecosystem services (Benayas et al. 2009). China is no exception: especially the Chinese Loess Plateau, one of the most severely eroded regions in the world, has been given a lot of attention in land restoration policies (Sun et al. 2014). Starting from the 1970s, the Chinese government has implemented several small-scale land restoration programs in the Chinese Loess Plateau to rehabilitate vegetation cover, combat desertification and reduce soil and water loss (Chen et al. 2015). In 1998, Wuqi county in the Yan'an area in the Chinese Loess Plateau started a pioneer land restoration program to reverse the ecological degradation by stopping cultivation of steep slopes and converting cropland and bare land to forest and grassland. One year later, based on the experiences in 1999, one of the world's largest-scale land restoration projects, the Grain for Green project (GGP), was initiated nationally covering more than 20 million hectares of cropland and bare land (Persson et al. 2013). One of the main purposes of the GGP is to maintain soil fertility and combat soil and water losses (Deng et al. 2019). The GGP brought a dramatic alteration of land use and a transformation in ecosystem services delivery (Chen et al. 2015).
Ecosystem services are defined as flows of materials, energy and information which are directly or indirectly provided by ecosystems to human society, including provisioning, regulating, cultural and supporting services (Costanza et al. 1997). Provisioning services include goods and products that we physically obtain from ecosystems, for example, food, water, raw materials etc.; Regulating services are necessary services to maintain the ecosystem functions, for instance, erosion control, sediment retention, habitat quality etc.; Cultural services like aesthetic landscape value provide spiritual pleasures to human beings (MEA 2005). Multiple ecosystem services can be provided by an ecosystem at the same time, but some ecosystem services cannot be supplied to society simultaneously . Any of these ecosystem services is associated with other services as either ''trade-offs'' or ''synergies'' (Bennett and Balvanera 2007). Trade-offs between ecosystem services can be comprehended as an increase of a (set of) specific ecosystem service(s) at the expense of other ecosystem services (Raudsepp-Hearne et al. 2010). Synergies, which are the opposite to trade-offs, are characterized as ecosystem services that either increase or decrease jointly (Bennett et al. 2009). Meanwhile trade-offs and synergies may appear diversely in one ecosystem at different temporal and spatial scales (Power 2010). Understanding the dynamics of ecosystem services is thus essential to comprehend the possible formation of trade-offs and synergies (Dade et al. 2019). Ignoring dynamics may increase the risk of unexpected changes in ecosystem services (Gordon et al. 2008). Human activity is a major factor affecting ecosystem service trade-offs and synergies through changing land use, by scale, type and intensity (Tolessa et al. 2017;Chen et al. 2020). For example, urbanization, ecological engineering and landscape restoration are often accompanied by a shift in land use for the purpose of (re-) generating a single or multiple ecosystem services. The implementation of these land use changing activities could cause ecological degradation if the trade-offs with other ecosystem services are ignored (Groot et al. 2011). Many previous studies were focused on simple trade-off and synergy relations between ecosystem services and ignore exploration of drivers and mechanisms. Ecosystem service bundles are defined as a mix of correlated ES provided at the same location and at the same time, though they may not have any direct causal relationships (Renard et al. 2015). The application of the ecosystem service bundles concept is helpful to understand the provisioning mechanisms of ES and the dynamics among multiple ecosystem services (Saidi and Spray 2018). Understanding of ecosystem service bundles is essential for achieving better management of multi-functional landscapes and minimizing costly trade-offs (Spake et al. 2017). Therefore, in this study, the ecosystem service bundles approach is used to understand the GGP impacts.
Impacts of the GGP on trade-offs and synergies between multiple ecosystem services have been investigated in multiple scientific studies. Many previous studies have analyzed ecosystem services supply in the Chinese Loess Plateau, and several address the dynamics and relations between different ecosystem services. For example, Lü et al. (2012) discovered that the entire Chinese Loess Plateau was transformed from a carbon source to a carbon sink by mapping carbon sequestration dynamics from 2000 to 2008. Feng et al. (2017) found out that vegetation coverage and types are the main factors that affect soil erosion control, soil moisture conservation and carbon sequestration based on field experiments in 2014. However, the majority of previous studies mainly put emphasis on single trade-off and synergy relations between regulating services, such as soil retention, water retention and water purification, ignoring the changes of other ecosystem services and the driving forces behind such changes. Meanwhile, some researchers focus on comparing ecosystem services after a certain time period against a baseline, but neglect the dynamics during that time span . For example, Li et al. (2019) mapped changes of ecosystem services in the entire Loess Plateau from year 2000 to 2015, without describing the fluctuation of ecosystem services within these 15 years. Besides, cultural services, which are defined as the nonmaterial benefits people obtain from the ecosystem, were not taken into account in the studies on the Loess Plateau. Furthermore, due to vegetation growth and continuance of the GGP, there is a lack of research considering the most recent impacts of the GGP on the ecosystem services on the Loess Plateau. We identified a knowledge gap in the long-term evolution of ecosystem services due to large-scale land restoration. Thus, in order to monitor the impacts of the GGP across various categories of ecosystem services, we considered the time period 1990-2018 (i.e., including ex-ante and ex-post phases of the GGP project) and selected 11 ecosystem services covering four provisioning, four regulating and three cultural services.
The implementation of the GGP is expected to have affected a range of ecosystem services on the Loess Plateau. The GGP proposed a reduction in cultivated area in return for an increase in forest and grassland area. Provisioning services, such as grain, livestock, fruit and timber were assessed in order to quantify the impacts from GGP land restoration measures. The main goal of the GGP is to prevent soil and water loss and maintain soil quality; thus, we included sediment retention and carbon sequestration as ecosystem indicators in our analysis. Additionally, it has been found that land restoration plays an important role in the reduction of surface streamflow in the Chinese Loess Plateau (Chen et al. 2020). Therefore, seasonal water yield, as an indicator for water supply, was also considered in this study. Furthermore, a primary goal of restoration is the protection of biodiversity, including genes, species, populations, habitats and ecosystems (Hector and Bagchi 2007); therefore, habitat quality is also quantified. Here, we define habitat as ''the resources and conditions present in an area that produce occupancy-including survival and reproduction-by a given organism'' (Hall et al. 1997, p. 175). Cultural services, like all other ecosystem services, must demonstrate unique relations between ecosystem structures and meeting the satisfaction of human needs (Daniel et al. 2012). Cultural services, including outdoor recreation, aesthetic value of the landscape and learning and inspiration values were considered.
The main objectives of this study are (a) to analyze the spatial and temporal dynamics in ecosystem services before and after the implementation of the GGP using ecosystem service bundles; and (b) to understand trade-offs and synergies between multiple ecosystem services.

Study area
The study area of Yan'an is located in the northern Shaanxi province on the south-central part of the Chinese Loess Plateau at latitude 35°21 0 -37°31 0 N and longitude 107°41 0 -110°31 0 E. Yan'an is a prefecturallevel municipality covering an area of 37,030 km 2 . It is a typical hilly area on the Loess Plateau that consist of multiple deeply incised valleys. The main soil type is Calcareous Cinnamon Soil (Xu et al. 2020). The terrain of Yan'an is higher in the northwest (highest point: 1795 m) and lower in the southeast (lowest point: 353 m), having an average elevation of around 1200 m (Fig. 1). Yan'an belongs to a semi-humid, warm temperate climate zone with continental monsoon circulation. The average annual temperature is 9.9°C and annual precipitation is 511 mm. In 1998, Yan'an area was selected as the first experimental site to start the GGP land restoration project in its northwestern Wuqi county. Up to now, Yan'an has implemented vegetation restoration for nearly 20 years and restored around 7200 km 2 of degraded land (Guo and Gong 2016).

Data sources
Land use and land cover (LULC) data of Yan'an area at a 30 m resolution for the years 1990, 1995, 2000, 2005, 2010, 2015 and 2018 was provided by the Data Center for Resources and Environmental Sciences of the Chinese Academy of Science (http://www.resdc. cn). This data was extracted from remote sensing data of Landsat-TM/ETM and Landsat 8. LULC data was classified into six classes: cropland, forest, grassland, water body, urban land and bare land. Meteorological data from 1990 to 2018, including precipitation, solar radiation, temperature, humidity and evapotranspiration, was obtained from the National Meteorological Administration of China (http://data.cma.cn) for 21 meteorological stations (see Fig. 1). A 30 m resolution DEM of Yan'an was obtained from the ASTER Global Digital Elevation Model (ASTER GDEM) from the Geospatial Data Cloud site of the Computer Network Center of the Chinese Academy of Science (http://www.gscloud.com). A soil erodibility map of Shaanxi province was obtained from the National Earth System Science Data Center (http://geodata.cn) and a world rainfall erosivity index map was acquired from the European Soil Data Center (ESDAC); http:// esdac.jr.ec.europa.eu). Additionally, world soil group data was obtained from EARTHDATA from NASA (http://earthdata.nasa.gov). Statistical data of the 13 counties in Yan'an was derived from the Statistical Yearbook of Yan'an from the Yan'an Statistical Bureau (http://tjj.yanan.gov.cn/). More details on parameters used in the InVEST model can be found in Appendix 1.

Quantification of ecosystem services
Eleven ecosystem services were selected to monitor the impacts of the GGP land restoration project in the 13 counties of Yan'an area (Table 1). Each ES was quantified in a biophysical way for the 13 individual Indicators for the four provisioning services were derived from the statistical yearbooks. As an indicator for grain production (GAP), the average yield of wheat and corn of each county (in t/km 2 ) was used. Apple yield (t/km 2 ) was used as an indicator for fruit production (FUP). Livestock production (LVP) was indicated by pork, beef and mutton meat productivity (t/km 2 ). Timber production (TBP) was indicated by the weight of timber produced per hectare (t/km 2 ). Four regulating services, including carbon sequestration (Mg/ha), sediment retention (t/ha), seasonal water yield (mm of base flow) and habitat quality (index from 0-1), were assessed by the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model (Nelson et al. 2018), which is explained in detail in ''Carbon sequestration'', ''Seasonal water yield'', ''Sediment retention'', and ''Habitat quality''. Indicators for the three cultural services were obtained from the statistical yearbooks and the LULC maps, respectively. Terrace area (%) was used as an indicator for the aesthetic value of the landscape (AVL), forest area (%) offered an underpinning for outdoor recreation (OR) and the number of local cultural institutes (n/1000 km 2 ) for entertainment and cultural education as an indicator for learning and inspiration (LAI). Additionally, the gross value of agriculture, industry and forestry (in USD/km 2 ) as well as population density (in person/km 2 ) were calculated from the statistical yearbooks as covariables.

Carbon sequestration
Carbon sequestration (CAS) was calculated based on the carbon storage and sequestration model from InVEST (version 3.7.0). This model is composed of three parts to calculate the carbon storage (Eq. 1): (1) carbon from plants including aboveground biomass and belowground biomass; (2) carbon from soil; (3) carbon from dead litter. Based on this calculation, land use and land cover change contribute mostly to changes in carbon storage due to changes in vegetation types.
To run this model, land use maps and carbon pools which indicate carbon storage values of different land use types are required. In this study, carbon sink data is based on experimental field data collected in Yan'an area: aboveground biomass data was obtained from  Zhang et al. (2001).

Seasonal water yield
Because Yan'an has a typical seasonal climate where precipitation is usually concentrated between July and September (Yang et al. 2018), the seasonal water yield model from InVEST was used to estimate water yield of the 13 counties in Yan'an. This model represents seasonal water yield (SWY) using two indices: quick flow and base flow. Quick flow indicates the generation of streamflow of hours to days, whereas base flow is defined as the generation of streamflow of months to years (Nelson et al. 2018). In order to monitor yearly water yield and reduce the climate variability impacts from fluctuating precipitations, base flow (in mm) was used while quick flow was excluded in this study.
The SWY model requires a series of monthly evapotranspiration (ET 1 -ET 12 ) maps, monthly precipitation (P 1 -P 12 ) maps, DEM, LULC maps, soil groups and integer Curve Numbers (CN). Monthly evapotranspiration was calculated with the R-package evapotranspiration (version 3.6.2) using meteorological data. Raster maps for monthly evapotranspiration and precipitation were created using the kriging tool in ArcGIS (version 10.5), based on the locations of the 21 meteorological stations within and surrounding the Yan'an area. CN data was obtained from the Hydrology Nation Engineering Handbook of United States Department of Agricultural (https://directives.sc.egov. usda.gov/17758.wba). To calibrate the model, we used streamflow data from Yan River (the river passing through the study area) as a reference to adjust the Threshold Flow Accumulation parameter, and this parameter was set to 800 by calibration. Default settings were used for other model parameters (alpha_m, beta_i and gamma).

Sediment retention
Sediment retention (SDR) in Yan'an area was calculated using the sediment delivery ratio model from InVEST. This model is a spatially explicit model based on the spatial resolution of the input DEM raster map. The calculation of the sediment delivery ratio consists of two parts (Nelson et al. 2018). The first part computes the annual soil loss from each pixel in the raster map based on the Revised Universal Soil Loss Equation model (RUSLE; Renard 1997), the RUSLE model is explained as below: where usle i is the amount of annual soil loss in one pixel, R i is the rainfall erosivity which is derived from the world erosivity map, K i is derived from the soil erodibility map, LS i is the length-gradient factor (calculated from the DEM), and C i and P i are the crop management and support practice factors, respectively, which were obtained from Fu et al. (2005). The second part generates the portion of soil loss that eventually reaches the stream and accounts for the final water yield results (Bhattarai and Dutta 2007).

Habitat quality
Habitat quality (HBQ) was quantified using the InVEST habitat quality model. This model combines information from the LULC map and disturbances to biodiversity to generate a habitat quality map (on an indexed scale between 0-1, where 1 indicates a perfect habitat to live). Both the impacts from biodiversity disturbances and the distance between the habitat and the threat sources are considered in the model. Biodiversity disturbances of both negative and positive induced sources were accounted. Negative sources include mining areas, roads, railways, urban areas and other populated areas, whereas positive sources contain natural reserves and national parks. The dynamics in biodiversity threats over the 1990-2018 time period were presented by threat maps that varied over time. These threat maps were obtained from the Worldmap dataset of Harvard University (https://worldmap.maps.arcgis.com/).

Statistical methods
Based on the above models and statistical data, we quantified 11 ES, including four provisioning services, three regulating services and four cultural services. In order to understand the spatial and temporal dynamics of each ES, we used the Space-Time Interaction (STI) method from Legendre et al. (2010). This method tests for space-time interactions in repeated ecological data, where there are no replications at the level of individual sampling units. In STI, variables of time and space are coded by principle coordinates of neighbor matrices into a two-way analysis of variance (ANOVA) model (Renard et al. 2015). A significant result of STI (p \ 0.001) indicates that the spatial distribution of an ES has changed over time. In our study, STI was processed using the package adespatial in R (version 3.6.3), setting each STI test at 999 permutations.
Each ES was calculated based on its mean ± SD across all 13 counties and taking the average value based on county area. Synergies and trade-offs between various ecosystem services were analyzed using Pearson correlation analysis in R (version 3.6.1). For every research year, the average value of each ecosystem service was defined at county level and each ecosystem service was standardized to a comparable unit scale from -1 to 1. Correlations between different ES were determined for the study period 1990 to 2018. Ecosystem service bundles were subsequently defined to assess the dynamics of multiple ecosystem services jointly. The 11 ecosystem services from one specific year and county were considered as an entity to be categorized across seven time intervals and 13 counties (i.e. 91 entities in total). Ecosystem service bundles were categorized using k-means cluster analysis from the package cluster in R (version 3.6.3). Clustering was attempted with four, five and six clusters, and the best result was kept. Maps with ecosystem service bundles were visualized using ArcGIS (version 10.5). Additionally, in order to understand the dominant patterns of ecosystem services values among different temporal and spatial scales, principle component analysis (PCA) was applied through R package ggplot (Wold et al. 1987).

Spatial and temporal dynamics in ecosystem services
In Fig. 2, the spatial and temporal dynamics of the 11 ES are presented that resulted from the STI analysis. For provisioning services, an obvious increase in fruit production was observed from 2005 to 2015. Results from the STI analysis (p \ 0.001) indicate that this increase only happened in a few specific counties. Livestock production almost doubled from 1990 to 2018 and this increase occurred across all counties.
Grain production fluctuated in all counties during the research period. A drop in grain production was observed in 2000, followed by an increase. In contrast, timber production showed a clear drop starting from 1995 up until 2005. During this time period, timber production decreased with almost 80%. After 2005, the production level tended to stabilize.
For regulating services, a gradual increase was generally observed in sediment retention, carbon sequestration and habitat quality. This gradual increase was not covering all counties, but took place in several specific counties (see Appendix 1-3); a significant p \ 0.001* from STI test results was found for three regulating services (CAS, SDR and HBQ). Meanwhile the highest increase was determined in 2005 in both regulating services CAS and HBQ. Trends for water yield were fluctuating. Water yield dropped in 1995 and increased again in 2010 and 2018. These fluctuations in seasonal water yield occurred in all counties from 1990 to 2018 as was illustrated by the STI test results (p value = 0.052; see Appendix 4).
Cultural services, such as habitat quality and outdoor recreation, showed similar increasing trends as the three regulating services. The values for outdoor recreation, aesthetic landscape value, and learning and inspiration all increased from 1990 to 2018. Results from the STI test (p value \ 0.001*) indicate that changes in outdoor recreation and aesthetic value of the landscape only occurred in specific counties of Yan'an area (see Appendix 5 for outdoor recreation), while learning and inspiration improved in all 13 counties. Overall, the spatial and temporal dynamics of the 11 ecosystem services indicated that the majority of the selected services showed an increasing trend. Only trends for timber production decreased clearly, while water yield decreased from 1990 to 1995 and increased from 2005 onwards.

Trade-offs and synergies between ecosystem services
In the trade-offs and synergies analysis of the ecosystem services, we found that the majority of ecosystem services showed synergistic relations. We used the average value of 11 ecosystem services in each research year at county scale. In Fig. 3, linear correlations between all ecosystem services are displayed, ordered by size of the Pearson correlation coefficient (r). Positive correlations indicate a synergy between services (0 \ r \ 1; displayed in blue color in Fig. 3), while negative correlations indicate a tradeoff between services (-1 \ r \ 0; displayed in red color in Fig. 3). In general, the figure shows that the majority of the correlations are positive, indicating synergies between those ecosystem services. For instance, there are strong synergies between aesthetic landscape value, learning and inspiration, livestock production, carbon sequestration, outdoor recreation, fruit production, sediment retention and habitat quality.
For provisioning services, fruit production showed a strong synergy with the majority of other ecosystem services, except for a trade-off with timber production. Also, livestock production had a trade-off with timber production. Timber production had trade-offs with the majority of the other services, except for a synergy with grain production. Grain production showed no significant correlation with the majority of other services, besides a slight synergy with timber production. Regulating services, including carbon sequestration, habitat quality and sediment retention had synergies between each other. Water yield showed Fig. 2 Space and time interaction (STI) test results for 11 ecosystem services across the 13 counties in Yan'an from 1990 to 2018. Note: error bars indicate standard deviation, calculated based on the average value of ecosystem services from 13 counties in one specific year. GAP grain production; FUP fruit production; LVP livestock production; TBP timber production; SDR sediment retention; CAS carbon sequestration; HBQ habitat quality; SWY seasonal water yield; OR outdoor recreation; AVL aesthetic landscape value; LAI learning and inspiration trade-off correlations with the aesthetic landscape value and livestock production. As for cultural services, outdoor recreation showed synergies with the majority of the regulating services, while tradeoffs with provisioning services were observed. Learning and inspiration and aesthetic landscape value had similar correlations with other ecosystem services. Additionally, the aesthetic value of the landscape had trade-offs with water yield and timber production, while learning and inspiration only showed a significant trade-off with timber production. The highest synergies were found between carbon sequestration and outdoor recreation (r = 0.99), and between sediment retention and habitat quality (r = 0.98). Additionally, carbon sequestration showed very strong synergies with sediment retention and habitat quality (r = 0.98 and r = 0.97, respectively). Highest tradeoffs were found between timber production and outdoor recreation (r = -0.89), followed by timber production and carbon sequestration (r = -0.85).

Principle component analysis between ecosystem services
By a combined analysis of a PCA and the Pearson correlations, the internal structure and explained variance of the trade-offs and synergies between different ecosystem services was investigated. The result of the PCA can be found in Fig. 4. Component 1 (PC1) explained 70.1% of the total variance while component 2 (PC2) explained 11.4%, apparently the summed variance of PC1 and PC2 had met the 60% threshold. The PCA result indicates how data is distributed along different dimensions (PC1 and PC2) and illustrates the direction of variances in the data: a smaller angle between variables signifies a more positively correlated data relation (Yang et al. 2004). PC1 occupied a major portion of the PCA test. Within PC1, besides timber production, gross forestry value shows a negative correlation to other ecosystem services as well. Additionally, in PC2, we found a Abbreviations can be found in Fig. 2 negative correlation between water yield and grain production, which was not observed in the Pearson correlation test.

Ecosystem service bundles
Seven time intervals in 13 counties of Yan'an area were considered in the categorization of ecosystem service bundles. Based on the cluster plots generated from R, we selected five clusters where the ecosystem service bundles were identified most clearly compared to four or six clusters. In Fig. 5, the specific components of the five ecosystem service bundles are displayed. For each bundle, the dominant ecosystem services were used to name the bundles. The five bundles are labeled as food production, sediment retention, forest habitat, landscape value and timber production. These five ecosystem service bundles indicate the value distribution of 11 ecosystem services at county-level and in a certain research year. According to the value of ecosystem services in each bundle, bundle 1 Food production was dominated by provisioning services, led by fruit production, followed by grain and livestock production. Bundle 2 Sediment retention had the highest sediment retention value while the remaining 10 ecosystem services fluctuated. Carbon sequestration, habitat quality and outdoor recreation were the focal ecosystem services in bundle 3. The cultural services aesthetic landscape value and learning and inspiration were well represented in bundle 4 labeled as landscape value. Bundle 5 Timber production was led by timber production, followed by water yield and grain production.
In Fig. 6 the spatial and temporal distribution of the five ecosystem service bundles in 13 counties across 7 time intervals is displayed. In general, we can observe a change of overall color from 1990 to 2018: in 1990 the dominant ecosystem service bundles was Timber production, while in 2010 the bundles evolved to Sediment retention and changed to Landscape value in 2018 eventually. Starting from 1990, in the Northern part of Yan'an area, Timber production was the major ecosystem service bundle. In 1995, the distribution of the ecosystem service bundles remained almost the same. However, from 2000 we observe a transformation from Timber production to Sediment retention in Wuqi, Baota and Yichuan, and to Landscape value in Zhidan county. The distribution of ecosystem service  Table 2.

Discussion
Eleven ecosystem services in Yan'an area were quantified and their spatial and temporal changes were estimated over the period 1990-2018, enabling to assess the impact of the large-scale GGP land restoration project that was implemented from 2000 onwards. The trade-offs and synergies between these ESs were analyzed, and ecosystem service bundles were assessed. Based on the results, we observed increases in the majority of the ecosystem services from 1990 to 2018, and particularly dramatic increases of fruit production, habitat quality, carbon sequestration, learning and inspiration and outdoor recreation that occurred since 2000, coinciding with the start of the GGP project. Correlation analysis revealed relations between specific ecosystem services, and both trade-offs and synergies were observed. Synergies were found between sediment retention, carbon sequestration, outdoor recreation, fruit production, habitat quality, learning and inspiration, livestock production and aesthetic value of landscape, while a trade-off was found between timber production and water yield. The ecosystem service bundles results showed an obvious change since 2000, as the majority of the ecosystem bundles changed from timber production to landscape value.
The results of ES quantification were similar to previous studies on the Loess Plateau, and confirm that there were increasing trends of sediment retention and  (Yang et al. 2018). In the results of sediment retention and carbon sequestration in Fig. 2, we observe a drop from 1990 to 1995. This indicates an ordinary trend in the Yan'an area before restoration implementation, representing a general degradation trend on the Loess Plateau. Shortly after, the implementation of the GGP started and since 2000 both of these regulating services slightly increased. From the collected 4 carbon input indices (carbon in above- ground biomass, below-ground biomass, litter biomass and carbon in soil), we observe huge differences in aboveground and belowground biomass between cropland and forest: forest contains 10 times higher biomass values than cropland on the Loess Plateau (Xiao et al. 2016;Feng et al. 2017). Due to the traditional practice of removing crop residues after harvest, in the Carbon model the carbon content of the litter layer of cropland was set at 0. Hence the introduction of the GGP, through an increase of forest and reduction of cropland, has increased the local carbon storage of the Yan'an area.
According to the observational evidence from many regions in the world, land use and climate change are recognized to be two majors drivers affecting baseflow (Price 2011). In the Chinese Loess Plateau, we observed a dramatic drop of water yield in 1995, and the value kept being consistently low compared to 1990 until the end of the assessment period in 2018. The dropped water yield can be explained as the newly planted forest and grassland have caused an increase of both evapotranspiration and net primary productivity (Feng et al. 2016). Additionally, in recent decades a significant increase of extreme warm surface temperature and a decrease of average daily precipitation were observed in the Chinese Loess Plateau (Sun et al. 2016). Wang et al. (2015) initiated a research of human activity impacts on runoff and sediment transport in Yan River, which is the main river in the Yan'an area, and concluded that human activity is a main reason for runoff decline by changing the land cover. Meanwhile, according to the algorithm of the Seasonal water yield model, decline of precipitation and increase of temperature and evapotranspiration could be key factors to cause a decrease in water yield. HBQ increased around 7% from 2000 to 2010 and showed a slight decreasing trend between 2010 and 2018. In the InVEST model, habitat quality is calculated based on distance and the area of disturbances from the habitat, as well as sensitivity of land cover types to threats. In comparison, forest is less sensitive to threats than cropland and grassland (Nelson et al. 2018). From the land use change table in Appendix 6 we found that the urban area expanded more than two times by 2018 compared to 1990, while forest land continually increased from 1995 to 2010 and maintained almost the same value after 2010. This trend could be explained as land restoration leading to an expansion of forest area and increase of HBQ from 2000 to 2010; after 2010, reforestation stagnated while urban area expansion caused a slight decrease of HBQ.
Cultural services often relate to spiritual significance and landscape aesthetics (Daniel et al. 2012). It is hard to quantify and monitor the cultural services, especially when crossing a huge time span, due to the difficulties in understanding human emotions from the past. However, in this study, despite of data deficiency about cultural services, we quantified the amount and monitored the dynamics of cultural services in terms of outdoor recreation, aesthetic landscape value and learning and inspiration on the Loess Plateau. In previous publications researching ecosystem services' dynamics on the Loess Plateau region, cultural services were frequently neglected (Feng et al. 2017;Yang et al. 2018). Only a few studies have investigated dynamics of cultural services during the implementation of the GGP. Hou et al. (2017) only recorded a slight increase of recreation capacity from 2000 till 2010 in Baota distinct in Yan'an area. Similar results of outdoor recreation have been found in our study while, additionally, a decrease in 1995 had been observed. Tourism is one essential indicator of cultural services indicating the attractiveness of a landscape, which has been studied by many ecosystem researches (Raudsepp-Hearne et al. 2010;Remme et al. 2015). In this study, there is insufficient tourism data when tracing back to 1990; however, it is believed that in Yan'an area tourism coincides with outdoor recreation. According to the tourist numbers from the recent five years, Huanglong and Huangling counties received the most tourists among other counties and had the highest forest cover, suggesting that outdoor recreation is correlated with forest cover. An increase of aesthetic landscape value was observed from 1990 till 2018 indicating an expansion of terrace area.
Terraces not only bring unique scenery to the local landscape, but also stimulate crop yield. A field experiment on the Loess Plateau found that the yield of a 3-year-old terraced land was 27% higher than sloping farmlands (Liu et al. 2011). Based on the dramatic increase of learning and inspiration value from 1990 to 2015, we can speculate local people had paid more attention to their indigenous cultural learning and entertainment. Results of ecosystem service bundles displayed the temporal and spatial dynamics of ecosystem services before and after GGP implementation. From Fig. 6, it can be observed that starting from 2000, there was a transformation of ecosystem service bundles from Timber production to Sediment retention and Landscape value in northern Yan'an (particularly in Wuqi, Zhidan, Ansai, Baota, Zichang, Guanquan and Yanchang counties). After 2015, since the GGP policy in Yan'an area had altered from mainly reforesting land to maintaining the reforested land, change of land use types was minimized. From the ecosystem service bundles maps, we observe the general process of ecosystem services components change by land restoration in the Loess Plateau. In the first 10 years of the GGP from 2000 to 2010, there was an increase of regulating services at the expense of provisioning services. After 2010, cultural services surpassed regulating services to occupy the majority of the ecosystem service bundles. While during 2010-2018, based on Fig. 2, regulating services were not decreasing, it was the proportion of cultural services in ecosystem service bundles that increased. Meanwhile we observe that the ecosystem service bundles became stable after 2015 since there was no change in ecosystem service bundles between 2015 and 2018.
According to the land use change map from Appendix 7, it can be observed that the majority of the land use change occurred in the northern part of Yan'an while Huanglong and Huangling counties feature much less land use changes. During the GGP implementation, there was a decrease in cultivated land and an increase in grassland and forest area in return (Appendix 6 and 7). Therefore, it could be expected that grain production would be reduced due to the shrinkage of cropland. However, according to the results in Fig. 3 there was an increase of average grain production from 2000 to 2010. One explanation is that there has been an increase of grain productivity due to the improvement of agricultural technology as well as the terrace expansion; for instance, there has been an increased utilization of fertilizer on the Loess Plateau (Fan et al. 2005), and from 2000 till 2008, grain yields increased from 3.0 t/ha to 3.9 t/ha in the Chinese Loess Plateau (Lü et al. 2012).
Changes of economic factors in terms of gross agricultural value (GAV), gross industrial value (GIV) and gross forestry value (GFV) as well as population density (POD) were also included in the STI test (Appendix 8). Timber production plummeted after 1995 and was almost 5 times lower in 2005. This change may be due to the introduction of GGP policy that banned all tree felling activities, thereupon triggering a decrease of GFV from 1990 to 2005. Meanwhile, the fruit industry has blossomed from 2000 in parts of Yan'an area, especially in Luochuan county, which is famous for its high quality and quantity of apple production (Ma et al. 2015). Additionally, according to the GGP strategy there were two types of forests restored from cropland and bare land: economic forest and ecological forest. Economic forest contains various species of fruit trees, nut trees and pulpwood which support local farmers' income, for instance, apple, pear, red dates and walnuts, whereas for ecological forest restoration usually drought-enduring trees and shrubs are selected, such as Robinia pseudoacacia, Hippophae rhamnoides and Platycladus orientalis (Deng et al. 2014). Therefore, an increasing area of restored economic forest expanded the fruit tree area simultaneously and improved fruit production as a result.
To sum up, implementation of the Chinese land restoration project GGP not only improved the majority of ecosystem services in the Chinese Loess Plateau, but also lead to local economic growth through subsidies and agricultural products. Results of this study are coherent with the ''4 Returns approach'', since landscape restoration is expected to enhance and restore ecosystem functions which leads to improved delivery of ecosystem services and the returns of natural capital, social capital, financial capital and the return of inspiration (Moolenaar 2016). In this study, we applied ecosystem service bundles to understand the dynamics of ecosystem services in both temporal and spatial scales, offering scientists and environment managers a method to systematically monitor changes of ecosystem services. Besides, we integrated three cultural services in the ecosystem services quantification to provide a more comprehensive view of ecosystem services in the study area, including both trade-offs and synergies of ecosystem services that nature provides to people. The latter echoes the concept of Nature's Contribution to People (NCP), which posits that both positive and negative contributions of nature to people should be considered with the aim to improve incorporation of NCP into policy and practice (Díaz et al. 2018). According to the guidelines of the ESP (Ecosystem Services Partnership), impact assessment and integrated cost-benefit analysis of land restoration are essential procedures to achieve sustainable landscape management and support land use planning (Groot et al. 2018). For instance, Groot et al. (2020) undertook an integrated cost-benefit analysis of large-scale landscape restoration in Spain. It is therefore suggested to conduct an integrated cost-benefit analysis of the GGP to unravel the social-economic and environmental impacts of land restoration in the Chinese Loess Plateau in a structured and coherent way. Such an analysis would allow to value the observed trade-offs in ecosystem services, as well as assess the overall net benefits in relation to restoration costs.
This research covered an area of 37,000 km 2 and considered an assessment period of almost 30 years. Therefore, there are many factors that could change ecosystem services, for instance, population increase, urban expansion, climate change, etc. However, according to the discussion above, the GGP is understood as a major driver that changed the land cover and ecosystem services simultaneously. Overall, the GGP implementation has had positive impacts on enhancing a majority of provisioning, regulating and cultural services, while the GGP shows negative impacts on timber production and water yield. Decrease of timber production was mainly due to land management policy but may not be a severe issue, as it could be managed by timber import from other provinces. Another concern is the decrease of water yield, although due to the shrinkage of cropland area by the GGP project, the demand for agricultural water use has decreased at the same time. Liang et al. (2018) reported a decline of soil moisture after GGP implementation, with forest featuring lower moisture content than cropland and grassland, while revegetation on the Loess Plateau is considered as a main driver for the moisture decrease. Clearly, forest expansion has brought more pressure on water supply than grassland on the Loess Plateau with an average annual precipitation from 250 to 600 mm. In order to maintain local water supply, it is recommended that further landscape restoration plans balance the revegetation area of forest and grassland.

Conclusion
In this study, the dynamics of 11 ecosystem services in 13 counties from Yan'an area were quantified within a time range from 1990 to 2018 in order to assess the impact of the large-scale GGP land restoration project implemented from 2000. An increasing trend was found in the majority of the provisioning, regulating and cultural services including fruit production, livestock production, sediment retention, carbon sequestration, habitat quality, aesthetic landscape value, learning and inspiration and outdoor recreation while seasonal water yield and timber production showed decreasing trends. We observed synergies between regulating and cultural services, including SDR, CAS, HBQ and OR, while both trade-offs and synergies were found in provisioning services. TBP was negatively correlated with CAS, SDR, OR, HBQ and LAI whereas GAP showed synergies. Ecosystem service bundles revealed temporal differences from 2000 until 2015 as well as spatial differences between northern and southern Yan'an. The process of ecosystem services components change by the GGP was discovered to start with increasing regulating services at the expense of provisioning services, followed by cultural services exceeding regulating services and occupying the main proportion of ecosystem services. Implementation of the GGP is recognized as a key factor changing the land use and affecting ecosystem service bundles. To conclude, GGP implementation has improved the majority of regulating and cultural services whereas it constrained timber production and local water yields. This study reveals the dynamics of ecosystem services while land restoration occurred; this knowledge supports future land use planning and helps to maintain a balance between different ecosystem services. From this study, it is suggested for the Yan'an government to pay attention to local timber products balance, as well as balancing forest and grassland area to maintain sustainable water supply.
Acknowledgements This research was partially funded by Commonland Foundation; and the Chinese Scholarship Council [Grant Number 201807720055].

Declarations
Conflict of 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.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.