Ecosystem service bundle index construction, spatiotemporal dynamic display, and driving force analysis

ABSTRACT Introduction: Existing studies on ecosystem service relationships are mainly qualitative or semi-quantitative assessments, but lack of quantitative exploration of aggregated ecosystem services and their influencing factors. We mapped the distributions of 12 ecosystem services of Zhejiang Province in 2000 and 2015 at the district and county level, analyzed their relationships using Spearman’s correlation analysis, constructed ecosystem service bundle index (ESBI) for each district and county by structural equation model, and then through multiple linear regression, we explored factors associated with ESBI variations. Outcomes: Our results showed that (1) most ecosystem services were spatially clustered. There were synergies between individual ecosystem services in categories of provisioning and regulating services, respectively; (2) our proposed ESBI index system consists of overall index and sub-indices of provisioning, regulating, and cultural services. The higher the ESBI value, the more important the corresponding place for multiple aggregated ecosystem service provision. Compared to 2000, ESBI in 2015 distributed more unevenly, and the average dropped by 3.10%; and (3) the increase of ESBI was associated with its initial value, and four socioeconomic and natural factors; the decrease of ESBI was influenced by the initial value and six key socioeconomic factors. Discussion and Conclusion: Our proposed ESBI system has several advantages (e.g., scale free, flexible weighting, quantitative and continuous indices for further analyses, and alternative non-monetary solution) in understanding and managing relationships among multiple ecosystem services.


Introduction
In the past four decades, China has made marvelous progress in economy. From 1978 to 2018, gross domestic product (GDP) increased from 367.87 to 90,030.95 billion yuan, and urban population rate increased by 41.66% (from 17.92% to 59.58%). However, the roaring economy and accelerated urbanization based on intensive human activities such as land use change have induced a series of ecological and environmental problems (Seto, Guneralp, and Hutyra 2012). For example, excess reclamation and inadequate flood control services resulted in the 1998 catastrophic flood in Yangtze River basin (Wu et al. 2013). Overgrazing and decline of sand fixation service led to severe sandstorms in northern China (Yang and Tuo 2002). Flood control and sand fixation (Haines-Young and Potschin-Young 2018) are ecosystem services that human beings benefit from ecosystems (MA 2005). The deficiency of them can induce various disasters, resulting in dramatic losses of human well-being (Yang et al. 2015b).
So far, scholars have done plenty of work on ecosystem service assessment and mapping. Raudsepp-Hearne, Peterson, and Bennett (2010), Turner et al. (2014), and Roces-Diaz et al. (2018) evaluated ecosystem services in Canada, Denmark, and Spain, respectively, and declared that most services were spatially clustered, especially for provisioning and regulating services. Qiu and Turner (2013) also stated in their study that individual ecosystem service distributions were spatially aggregated. A number of Chinese studies (Jiang and Zhang 2016;Lyu et al. 2019;Wang et al. 2017;Wu et al. 2013;Yang et al. 2015a) mapped individual ecosystem service patterns as well. Some scholars proposed that ecosystem service patterns were spatially clustered (Lyu et al. 2019;Yang et al. 2015a) and found high levels of spatial similarity of some ecosystem services (Liu et al. 2019;Qin, Li, and Yang 2015), such as water yield, water interception, and soil retention in Guanzhong-Tianshui economic zone (Qin, Li, and Yang 2015).
Relationships between multiple ecosystem services can be negative (potential trade-offs) or positive (potential synergies). The trade-off occurs when one service increases at the cost of reducing the provision of another, and the synergy arises when multiple services strengthen or weaken simultaneously (Raudsepp-Hearne, Peterson, and Bennett 2010). In general, regulating services have synergies with each other (Dittrich et al. 2017;Li, Li, and Zhang 2016b;Qin, Li, and Yang 2015;Qiu and Turner 2013;Raudsepp-Hearne, Peterson, and Bennett 2010;Yang et al. 2015a) and have trade-offs with provisioning services Li, Li, and Zhang 2016b;Qin, Li, and Yang 2015;Raudsepp-Hearne, Peterson, and Bennett 2010;Turner et al. 2014;Wu et al. 2013;Yang et al. 2015a). However, Lyu et al. (2019) suggested different views that sand fixation had trade-offs with carbon sequestration and nutrient retention, and the latter two had synergies with crop production. There is no agreement on relationships between cultural services and other services. Some scholars found trade-offs between individual cultural services and food production (Raudsepp-Hearne, Peterson, and Bennett 2010;Renard, Rhemtulla, and Bennett 2015;Turner et al. 2014), but some came to the opposite conclusions (Roces-Diaz et al. 2018;Schirpke et al. 2019). Turner et al. (2014) believed cultural and regulating services had synergistic effects, and Qiu and Turner (2013) found forest recreation (related to the amount of forest habitat, recreational opportunities, proximity to population center, and accessibility) was synergistic with carbon storage and soil retention; however, Lyu et al. (2019) showed that recreational opportunity (relevant to accessibility, land use, scenic beauty, cultural sites, and potential activities like horseback riding, mountain climbing, recreational fishing, etc.) had trade-offs with carbon sequestration and nutrient retention. Such divergent viewpoints possibly are due to different ecosystem types and diverse ecosystem service indicators. In terms of supporting services, trade-offs were found between habitat provision and provisioning services , while there was a synergistic effect between habitat quality and crop production (Ai et al. 2015).
Some scholars attempted to assess ecosystem service aggregated patterns from two aspects. One is to identify hot and cold spots of ecosystem services. Wu et al. (2013), Li et al. (2016a), Liu et al. (2019), and Ran et al. (2019) respectively defined the values of top 10%, 20%, 50%, and above average as individual ecosystem service hotspots, and classified the research units according to the total number of hotspots inside. Qiu and Turner (2013) identified hot and cold spots as places containing 6 or more services in the highest and lowest 20%, respectively. Roces-Diaz et al. (2018) identified provisioning, regulating, and cultural service hotspots in Catalonia, Spain by Getis-Ord G � i method. Kim et al. (2018) determined the hot and cold spots of regulating services in Anshan, Korea based on perceptions of regional environmental stakeholders. Nevertheless, in the process of determining the overall hotspots distribution, the above research gave different ecosystem services the same weight, ignored the relationships between ecosystem services. Moreover, the results of hotspots analysis were discrete, which were not feasible for further analyses. These drawbacks to some extent increased the uncertainty and subjectivity of overall status of ecosystem services. The other main aspect of assessing aggregated ecosystem service patterns is to construct ecosystem service bundles. Ecosystem service bundles are sets of ecosystem services that repeatedly appear together across space or time (Raudsepp-Hearne, Peterson, and Bennett 2010). So far, several studies have evaluated the ecosystem service bundles, most of which used simple statistical methods such as principal component analysis (PCA), cluster analysis (CA) (Wu et al. 2015;Yang et al. 2015a), and self-organizing map (SOM) (Crouzat et al. 2015;Dittrich et al. 2017;Li et al. 2017;Mouchet et al. 2017). Lyu et al. (2019) identified three ecosystem service bundles using PCA, CA, and SOM in the urban belt along the Yellow River in Ningxia, China. Raudsepp-Hearne, Peterson, and Bennett (2010), Turner et al. (2014), Schirpke et al. (2019), and Yao et al. (2016) respectively evaluated the ecosystem service bundles in Quebec, Denmark, European Alps, and upper Hun River basin in China by PCA and CA. However, these studies can only reflect the overall situation of ecosystem services qualitatively or semiquantitatively.
There was little research on driving factors of change in aggregated ecosystem services. Ling et al. (2019) suggested that vegetation cover, human disturbance, temperature, and vegetation drought index had impacts on ecosystem service values, with contribution of 38.0%, 31.6%, and 30.4%, respectively.
Overall, existing studies on ecosystem service relationships are mainly qualitative or semi-quantitative assessments, but lack of quantitative exploration of multiple ecosystem services and related influencing factors. To fill some of previous research gaps and guide multiple ecosystem service management, we aim to: (1) display key ecosystem service patterns and demonstrate synergies and trade-offs between individual and among diverse categories of ecosystem services; (2) construct quantitative ecosystem service bundle index (ESBI) system and display its spatiotemporal dynamics; and (3) analyze factors associated with ESBI value changes.

Study area
Zhejiang Province (118°1′~122°26′E, 27°9′~31°11′N) is located in the southeastern coast of China ( Figure 1). As of 2015, it covers a total land area of 105,500 km 2 , including 11 cities, 90 counties, and 1350 towns. The terrain slopes down from southwestern mountains to northeastern plains, with 74.63% of hilly mountains, 20.32% of plains, and 5.05% of water bodies. Meanwhile, Zhejiang has 260,000 km 2 sea area with more than 3000 islands. Various landforms create diverse and abundant ecosystem services.
Zhejiang Province is one of the richest and most developed provinces in China (Zhao, Zhao, and Li 2019), and its economy has developed rapidly in recent years. From 2000 to 2018, GDP of Zhejiang increased by 13.09%, and exceeded 5,619.72 billion yuan in 2018. At the meantime, its proportion of urban population increased from 48.70% to 68.90%; the area of paved roads augmented from 90.50 to 459.46 million m 2 ; and the amount of compound chemical fertilizer used for agriculture nearly doubled (Zhejiang Statistical Yearbook, 2001. However, intense human activities also posed great threats to the supply of ecosystem services. Therefore, we chose Zhejiang as our study area, aiming to provide guidance for the government to better manage the ecosystem-economy balance, and also set a good example of multi-objective ecosystem management for other fast-developing areas in China and the rest of the world.

Data collection and analysis
We chose 12 types of ecosystem service that are indispensable for human well-being, covering provisioning, regulating, cultural, and supporting services (Table 1), specifically including crop production, meat production, aquatic production, cotton production, fruit production, carbon sequestration, water retention, soil retention, forest recreation, tourism infrastructure, tourist population, and biodiversity. We also collected data of potential socioeconomic and natural factors related to ESBI variations, including GDP in 2000, road density in 2000, precipitation in 2000, and so on (Table 2). We used administrative boundary (Ai et al. 2015;Li et al. 2017; Raudsepp-Hearne, Peterson, and Bennett 2010) of districts and counties to distinguish the research units. Due to data availability, we merged sub-administrative districts in the main urban areas of Hangzhou, Huzhou, Jinhua, Quzhou, and Jiaxing as separate larger research units, respectively. Meanwhile, we did not include island areas (i.e., Zhoushan City and Dongtou County of Wenzhou City) because of data scarcity for key indicators.
We normalized the values of ecosystem service indicators in 2000 and 2015 together to 0 ~ 100 by the minmax normalization method to ensure comparability across years. Using the "corrplot" package in RStudio (version 1.1.383), we conducted the Spearman's correlation analysis to demonstrate synergies and trade-offs among individual and diverse categories of ecosystem services. We constructed the ESBI system using Structural Equation Models (SEM) in Mplus (version 6.1) (Feng and Harring 2020). Then, we normalized the ESBI values and mapped their spatial distributions to show their spatiotemporal dynamics from 2000 to 2015.
SEM (Bollen & Noble, 2011) includes structural model and measurement model, and the general form of equations are as follows: The importance level of biodiversity protection where formula (1) is the structural model, formula (2) and (3) are the measurement models. � k is the overall index of unit k; η k is the vector of sub-indices of unit k; α η , α y , α x are intercept vectors; B is the correlation matrix among sub-indices; Γ is the weight matrix that sub-indices contributing to the overall index; y k and X k are the vectors of observing variables of η k and � k , respectively; Λ y and Λ x are the weight matrices of y k to η k and X k to � k , respectively; ς k , τ k and δ k are error terms.
In this study, the overall index is the ESBI; the subindices are ecosystem provisioning service set (EPS), ecosystem regulating service set (ERS), and ecosystem cultural service set (ECS); the observed variables are the 11 types of ecosystem service in Table 1 except biodiversity. We filtered observed variables for each sub-index based on the reliability test in Stata (version 14.1, StataCorp, College Station, Texas, USA), and modified the model by adding or removing paths between variables based on the modification indices and model fitness indices in Mplus. For the final model, the Comparative Fit Index (CFI) and Tucker-Lewis Index (TLI) often should be higher than 0.90; the Root Mean Square Error of Approximation (RMSEA) should be less than 0.08 (Wen and Liang 2015). Finally, we performed the multiple linear regression analyses to analyze factors associated with changes in ESBI values via R packages of "car" (Fox and Weisberg 2019) and "QuantPsyc" (Fletcher 2012). The dependent variable was the change in ESBI between 2000 and 2015, the independent variables are the initial value of ESBI and other factors shown in Table 2.

Spatial and temporal distributions of ecosystem services
Supporting and most provisioning services had higher values in the northeast, central, and coastal areas of Zhejiang, and the output of aquatic products in coastal regions was particularly high (Figure 2). Carbon sequestration in the north was relatively lower than that in the central and southern areas; water retention was higher in the southwest and lower in the northeast; soil retention was higher in the northwest, south, and southeast, and lower in the northeast. The distribution of forest recreation was similar to that of soil retention; tourism infrastructure scattered in Zhejiang, and tourist population in the south lagged behind that in the north.
From 2000 to 2015, most provisioning services and all regulating, cultural, and supporting services improved (Table 3 and Figure 2). For provisioning services, fruit production increased most significantly, which was 3.86 times higher in 2015 than that in 2000, with obvious increases observed in the central, northeast, and coastal areas; crop and cotton  production declined by 39.50% (4.85 million tons) and 26.08% (7,568.63 tons) respectively. Meanwhile, crop and cotton production were the only two with reduced overall outputs among the 12 ecosystem services. For regulating services, carbon sequestration increased the most (42.35%, 131.69 million tons), mainly distributed in the central, western, and southern regions. In contrast, water and soil retention increased only by 4.66% and 0.49%, and the latter changed the least among the 12 ecosystem services. For cultural services, tourism infrastructure and tourist population soared obviously. The average score of tourism infrastructure and the tourist population in 2015 were 15 times and 10 times of those in 2000, respectively. Tourist population in the middle and north increased faster than that in the south. For supporting service, biodiversity increased significantly in the whole province, with an average increase of 49.69% (0.80) in 76 districts and counties.

Trade-offs and synergies among ecosystem services
Among the 66 pairs of ecosystem services, 52 and 51 pairs were significantly correlated in 2000 and 2015, of which 27 and 30 pairs were highly correlated (Spearman's coefficient [r] ≥ 0.5), and 18 and 20 pairs were moderately correlated (0.3 < r < 0.5) (Figure 3). In both years, there were positive correlations among individual services in categories of provisioning (p < 0.01) and regulating (p < 0.001), respectively. Tourism infrastructure and tourist population in cultural services were positively correlated (p < 0.05), representing apparent synergies. In terms of different categories of ecosystem service, regulating services were negatively correlated with both provisioning (p < 0.05) and supporting (p < 0.001) services, indicating obvious trade-offs. Forest recreation in cultural services had trade-off with provisioning services (p < 0.01), but had synergies with regulating services (p < 0.001). There were synergies between provisioning and supporting services (p < 0.01). From 2000 to 2015, although the correlation direction was the same, the correlation degree changed, including 25 pairs of weakened ones and 22 pairs of enhanced ones. The synergies among crop, meat, and aquatic production weakened while the synergies between fruit production and meat, aquatic, cotton production enhanced. The trade-offs between crop, meat, aquatic production, and regulating services mostly became weaker, while the other two provisioning services' trade-offs with regulating services enhanced. The synergy between tourism infrastructure and tourist population strengthened. Except cotton, fruit production, and water retention, the relationships between biodiversity with other services weakened.

ESBI system and its spatiotemporal dynamics
As Figure 4 shows, ESBI overall index was significantly measured by three sub-indices: EPS, ERS, and ECS. The coefficients of EPS and ERS were both positive, while that of ECS was negative and the correlation between EPS and ERS was negative. The EPS had the greatest influence on ESBI, while the ECS had the least. The standardized coefficients of all indicators for subindices were positive (p < 0.01). EPS was measured by five provisioning services, of which the contributions of crop (standardized coefficient [SC] = 0.772) and meat production (SC = 0.722) were relatively high and that of aquatic production was the lowest. ERS was explained by three regulating services and one cultural service (forest recreation). Forest recreation contributed the most (SC = 1.010), while water retention the least. ECS was estimated by three cultural services and one regulating service (carbon sequestration), with the largest contribution of tourist population (SC = 0.948) and the smallest of forest recreation.
The spatial pattern of ESBI in 2000 was similar with that in 2015 (   (2) The subscripts of s and m refer to the sum value and the mean value, respectively; the numbers in the brackets are corresponding standard deviations; (3) CP, MP, AP, COP, FP, CS, WR, SR, FR, TI, TPOP and BD are abbreviations for crop production, meat production, aquatic production, cotton production, fruit production, carbon sequestration, water retention, soil retention, forest recreation, tourism infrastructure, tourist population, and biodiversity, respectively.
distribution of ESBI in 2015 was more uneven, with the standard deviation increased by 14.15% from 13.57 to 15.49, and the average decreased by 3.10% from 40.98 to 39.71. Table 4 shows the descriptive statistics of variables used for the construction of regression models. Table 5 shows that eight variables, covering initial value, socioeconomic, Blue and red circles represent synergies and trade-offs, respectively. Values in the lower triangles are correlation coefficients. The higher the absolute values of correlation coefficients, the darker the color, the bigger the circle, and the stronger the correlation. CP, MP, AP, COP, FP, CS, WR, SR, FR, TI, TPOP, and BD are abbreviations for crop production, meat production, aquatic production, cotton production, fruit production, carbon sequestration, water retention, soil retention, forest recreation, tourism infrastructure, tourist population, and biodiversity, respectively. * p < 0.05, ** p < 0.01, *** p < 0.001.  Table 6 suggests that eleven factors, including initial value, socioeconomic, and natural factors, collectively explained 74.64% of the total variance for the decrease in ESBI. Except road density in 2000 and natural factors, the other seven factors were statistically significant associated with the decrease in ESBI. According to standardized coefficients, household numbers in 2000 contributed the most across the model (SC = 1.302, p < 0.001), followed by GDP in 2000 (SC = −1.260, p < 0.001), gross product of agriculture, forestry, animal husbandry, and fishery (GPAF) in 2000 (2) *p < 0.05, **p < 0.01, ***p < 0.001; (3) the unstandardized coefficient of ERS was fixed as 1; (4) coefficients in the diagram are standardized with corresponding standard errors are in brackets; and (5) ESBI, EPS, ERS, and ECS represent the ecosystem service bundle index, ecosystem provisioning service set, ecosystem regulating service set, and ecosystem cultural service set, respectively; CP, MP, AP, COP, FP, CS, WR, SR, FR, TI, and TPOP are abbreviations for crop production, meat production, aquatic production, cotton production, fruit production, carbon sequestration, water retention, soil retention, forest recreation, tourism infrastructure, and tourist population, respectively. Please see more details in Appendix A.

Discussion
We selected 12 key ecosystem services in Zhejiang Province, displayed the spatiotemporal dynamics (between 2000 and 2015), and explored their potential synergies and trade-offs. Then, to explore the overall status of ecosystem services, we constructed the ESBI system and further analyzed the factors associated with changes in ESBI values. Below we compared our results with those of existing literature, summarized advantages and limitations of our study, and provided some insights for future research and management.
We found that at individual service level, provisioning services were synergistic with each other and so were regulating services; however, at the grouped service level, provisioning services had trade-offs with regulating services. These results were in line with those of Zhang et al. (2016), Xu et al. (2016), and Jia et al. (2014), but different from Lyu et al. (2019) and Ai et al. (2015). Such a discrepancy is probably due to the complexity of ecosystem service relationship, which is affected by ecosystem type, selected service indicator, ecosystem service assessment method, geographic and climate conditions, and other factors. Furthermore, in line with Yang et al. (2015a) and Lyu et al. (2019), our results suggest that there is no agreement on the relationships between cultural services and other services. Biodiversity in our study had tradeoff with soil retention, and this viewpoint is consistent with Sun and Li (2017), who estimated the relationship of soil conservation and biodiversity protection in Zengcheng District, Guangzhou City. However, Liu, Xu, and Jiang (2015) found that soil retention and biodiversity protection are synergistic in Wangqing County, Yanbian Autonomous Prefecture, Jilin Province. Biodiversity has synergies with provisioning services in our study may because plantations where most provisioning services are produced also provide   food and habitats for vulnerable and endangered animals and plants (Cox et al. 2006;Fang 2017;Schmid 2002). Our proposed ESBI system has several advantages in understanding relationships among multiple ecosystem services. The higher the ESBI values, the more synergies and the less trade-offs exist in the corresponding analysis units, the better aggregated ecosystem service provision, and vice versa. The increase in ESBI is beneficial while the decrease in ESBI is harmful. The advantages of our proposed ESBI system are as follows. First, as long as data are sufficient, our system can be applied to study areas at different scales. It provides possibility for other regions to study the aggregated pattern of ecosystem services at the scale of country, city, town, and even pixel. Second, it is a composite index system, including both the overall index and sub-indices. One can not only display the distribution of ESBI but also grasp the pattern of each sub-index. Third, the weighting of both overall index and sub-indices are objective and unequal. They are determined by the variation of data. If needed, one may also customize the weighting of any index by fixing its corresponding coefficient in the structural equation model (Yang, McKinnon, and Turner 2015c). Therefore, it is flexible and feasible for different management purposes such as providing most aggregated ecosystem services, providing most provisioning services, or developing tourism and minimizing the impacts on regulating services. Fourthly, the values of computed ESBI are continuous rather than categorical, and thus are easily available for further quantitative analyses. Finally, in contrast to monetizing all selected ecosystem services, our proposed ESBI system provides an alternative non-monetary solution to assess the aggregated pattern and relationships of multiple ecosystem services. Thus, our approach also avoids the disputes on discrepancy of different economic valuation methods (Yang et al. 2008).
There are also some limitations in our study. First, although we have taken data quality control measures, our results may be influenced to some extent due to the data scarcity and potential errors in some of our selected indicators obtained from government statistic yearbooks. Second, ecosystem service patterns may be dependent on the indicators used. For example, results may be different when taking forest land proportion and the number of hunting animals as the proxies of forest recreation service. Third, positive and negative correlations between ecosystem services suggest potential synergies and trade-offs, respectively. Fourth, we did not include policy factors in determining the drivers of changes in ESBI values.
We identified several issues for further research. First, it would be helpful to better understand the impacts of cultural services on overall ecosystem by dividing tourism attractions into different categories, such as natural spots, historic relics, and wildlife spots. Second, it would be interesting to understand how policy factors may influence the relationships among multiple ecosystem services as well as the ESBI values as a whole. Finally, future research may also use the overall index and sub-indices of ESBI for further analyses such as hotspot/coldspot identification, spatial planning, and priority settings for conservation and development.

Conclusions
We mapped the distributions of 12 ecosystem services in 2000 and 2015, analyzed the relationships among ecosystem services by Spearman's correlation analysis, constructed ecosystem service bundle index system via structural equation model, and explored driving forces of ESBI spatiotemporal variations using multiple linear regression. Our main findings include: (1) most ecosystem services were spatially clustered. There were synergies between individual ecosystem services in categories of provisioning and regulating services, respectively. Regulating services had trade-offs with both provisioning and supporting services, while there were synergies between the latter two categories; (2) Our proposed ESBI index system consists of overall index and sub-indices of provisioning, regulating, and cultural services. The higher the ESBI value, the more important the corresponding place for multiple aggregated ecosystem service provision. ESBI overall indices were higher in the north central, west, and southeast, and lower in the northwest corner and central areas. Compared to 2000, ESBI in 2015 distributed more unevenly, and the average dropped by 3.10%; and (3) the increase of ESBI was associated with the initial value, and four socioeconomic and natural factors; the decrease of ESBI was influenced by the initial value and six socioeconomic factors. Our study may give guidance to better realize multi-objective ecosystem service management and lay a solid foundation for further study. (2) * p < 0.05, ** p < 0.01, *** p < 0.001; (3) the coefficient of ERS was fixed as 1. The EPS and ERS was negatively correlated with a coefficient of −21.973 and standard error of 7.789 (p < 0.01); and (4) ESBI, EPS, ERS, and ECS represent the ecosystem service bundle index, ecosystem provisioning service set, ecosystem regulating service set, and ecosystem cultural service set, respectively; CP, MP, AP, COP, FP, CS, WR, SR, FR, TI, and TPOP are abbreviations for crop production, meat production, aquatic production, cotton production, fruit production, carbon sequestration, water retention, soil retention, forest recreation, tourism infrastructure, and tourist population, respectively.