Study on the Comprehensive Improvement of Ecosystem Services in a China’s Bay City for Spatial Optimization

: Ecosystem services are characterized by region and scale, and contribute to human welfare. Taking Yantai city, a typical bay city in China, as the example, its three representative ecosystem services: food supply (FS), carbon sequestration (CS) and water yield (WY) were chosen as study targets. Based on analyzation of six different aspects of the supply and variation characteristic of demand, this study tried to propose advices for comprehensive improvement of ecosystem services for spatial optimization. The results showed that: (1) ecosystem services supply was strong in central and southern areas of Yantai, while the northern coastal areas were relatively weak; (2) synergistic relationships were found of FS-CS, FS-WY and CS-WY both in 2009 and 2015, with the strongest one for FS-WY. Additionally, in the synergistic relationships, each pair of ecosystem services was dominated by one ecosystem service; (3) most of the three pairs of synergistic relationships had the tendency to strengthen with larger scales; (4) four ecosystem demands changing areas were observed and comprehensive improvement suggestions for them were proposed. This work provides a new attempt to improve ecosystem services based on its supply-demand relationship, which will give a baseline reference for related studies in Yantai city, as well as other similar bay cities.


Introduction
Ecosystem services are the various direct or indirect benefits that people derive from ecosystems, including provisioning, regulating, cultural and supporting services, which are important links between ecosystems and human well-beings [1][2][3]. Influenced by the diversity of ecosystem services, the complexity of ecological processes and human preferences, different ecosystem services exhibit complex interactions, which can be described as trade-offs and synergies [4,5]. Trade-offs between multiple ecosystem services occur when the improvement of one ecosystem service is at the expense of the reduction of another ecosystem service. The trade-offs are common between provisioning services, and between provisioning and regulating services [6,7]. Conversely, the relationship is synergistic when an improvement in one ecosystem service leads to an increase in another [8]. The synergy is relatively common for regulating, supporting and cultural services [9,10]. In recent years, the contradiction between the excessive human demand for material and resources and the limited natural resources has become increasingly evident. Meeting the needs for food, water, materials and shelter for more than 7 billion people on the planet comes at the expense of natural ecosystems [11,12]. Rapid urbanization changes the land use pattern, services and their changes [62], ecological risks [63] and ecological efficiency in Bohai Rim [64]. While using authoritative land-change survey data to analyze the spatiotemporal characteristics of trade-offs/synergies, scale effects, driving mechanisms and improvement of ecosystem services is lacked. Therefore, it is necessary to explore the spatial patterns of ecosystem services trade-offs and synergies, and propose some improvement suggestions in the context of rapid socio-economic development in Bohai Rim. Yantai City, a typical city in the Bohai Rim, was selected as the study area. Even though as a bay city, water shortage is always a big problem for it. Per capita water resource is lower than 500 m 3 in over a half of the city. Meanwhile, carbon emission is large in the city because of large proportion of industrial production and mining. However, the city has relatively developed agricultural production, especially in fruits planting [65]. Thus, food supply, carbon sequestration and water yield were chosen as the representative ecosystem services in this study. Based on the land-change survey data, the spatiotemporal distribution characteristics of the three types of ecosystem services were analyzed. Through correlation analysis and the introduction of an ecosystem service trade-offs/synergies index, the spatiotemporal characteristics, scale effects and driving mechanisms of the trade-offs/synergies among the three types of ecosystem services were investigated. Meanwhile, on the basis of tradeoffs/synergies and changes in ecosystem services demand, improvement suggestions of ecosystem services were proposed to provide decision support for ecosystem conservation and optimal management in Yantai and other bay city.

Study Area
Yantai is located in the eastern part of the Shandong Peninsula, China, and is an important port city in the Bohai Rim. The city governs 12 county-level administrative regions, with complex terrain conditions, as shown in Figure 1. The elevation is high in the middle of the territory, and low in the south and north, with the range of 0~922.8 m. The climate belongs to the warm temperate continental monsoon climate. It has an annual average temperature of 13.0 • C, a sunshine hour of 2611.6 h and an average annual precipitation of 680 mm. The precipitation is concentrated in summer (from June to September), accounting for more than 70% of the annual total precipitation. However, the spatial distribution of rainfall is uneven, which is significantly more in the central and southern hilly than the northern coastal areas. The water volume of rivers is greatly affected by rainfall, and drying-out of rivers are common to see during non-flood season. Water resources are extremely scarce in the city, with an average total water resources of only about 3.2 billion m 3 over the years. Additionally, per-capita water resource is less than 500 m 3 , which is only one-fifth of the national average level. The lack of water resources has become an important limiting factor for future development. As the acceleration of urbanization, there are a large number of industrial zones, ports and docks under construction. Compared with other coastal cities, the proportion of industrial land in Yantai is relatively large, which has a greater negative impact on the environment. The urban ecological space has been drastically reduced, and environmental pollution has become more and more serious.

Data Sources
The land-change survey data in 2009 and 2015 was from Yantai Bureau of Land and Resources, including 8 primary types (arable land, garden plot, woodland, water area, construction land, grassland, traffic land and other land), and 33 sub-class types. The annual rainfall and Normalized Difference Vegetation Index (NDVI) data were obtained from the data center of resources and environment science, Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 1 August 2020)), with resolution of 1000 m× 1000 m. Evapotranspiration data was from MOD16 (http://www.ntsg.umt.edu/project/modis/ mod16.php (accessed on 1 August 2020)), with resolution of 1000 m × 1000 m. Soil data were obtained from national soil survey data with resolution of 1000 m × 1000 m. Total water consumption and total population in 2009 and 2015 came from Shandong Yantai statistical yearbook (http://www.shujuku.org/shandong-statistical-yearbook.html (accessed on 1 August 2020)). The total output of grain, oil crops and vegetables, per capita food demand came from Yantai statistical yearbook (http://tjj.yantai.gov.cn/col/col118 /index.html (accessed on 1 August 2020)). Population density data is from WorldPop (https://www.worldpop.org/ (accessed on 1 August 2020)), with resolution of 100 m × 100 m. China's per capita carbon emissions are adopted for per capita carbon emission, coming from Emissions Database for Global Atmospheric Research (EDGAR). It should be noted that land-change survey data was from the Second National Land-Change Survey of China. The data is obtained through field investigation on the land use status, land ownership and administrative division changes according to the land use status database and remote sensing image data issued by the superior department. The survey is a remarkably comprehensive survey, as a statutory survey in China and it employs highresolution remote sensing images (e.g., <1-m spatial resolution sensing images in the eastern regions, or <5 m in the western regions), aided by the field surveys [13]. The software used was ArcGIS 10.6 (Esri, Redlands, CA, USA) and Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) 3.7.0 (Natural Capital Project, Stanford, CA, USA).

Estimation of Ecosystem Services Supply and Demand
Ecosystem services supply and demand was estimated by the following formulas using the land-change survey data in 2009 and 2015. Based on this, changes in the supply and demand can be obtained.

Water Yield Service
Freshwater supply is an important provisioning service, which is of great significance for the development of irrigated agriculture, ensuring population growth, the development of industry and tourism and the improvement of living standards. In this study, the water yield (WY) model embedded in InVEST was used to estimate the annual average water yield. The model is mainly based on the Budyko hydrothermal coupling equilibrium hypothesis, which assumes that all water except evapotranspiration reaches the outlet of the basin [66,67]. Using the land-use, climate and soil data, the annual water yield of the study area can be obtained. Based on total water consumption and population data, we got the per capita water consumption. Combined with population density, water demand (WD) was obtained. The formula is as follows: where WY refers to annual water yield (mm); AET(x) is the annual actual evapotranspiration of grid unit x (mm); P(x) is the annual precipitation amount of grid unit x (mm); PET(x) is the potential evapotranspiration of grid x (mm); Kc(x) is the plant (vegetation) evapotranspiration coefficient; ETo(x) is based on the evapotranspiration of a reference vegetation, which reflects the local climatic conditions; AWC(x) is the volumetric (mm) plant available water content; W(x) is a non-physical empirical fitting parameter that characterizes the land surface properties of catchments; Z is an empirical constant, which captures the local precipitation pattern and additional hydrogeological characteristics [68]; WD is water resource demand; Dpcwc is per capita water consumption; Ppop is population density. The biophysical coefficients of each land use class used in InVEST water yield model can be found in Table 1, including the information about vegetation, plant evapotranspiration coefficient Kc and root depth for each land use class. Land use_veg shows the information used in AET equation. The value should be 1 for vegetated land, and 0 for all other land. The Kc is the plant evapotranspiration coefficient for each land use class, and it estimates plant evapotranspiration based on FAO 56 guidelines [69,70]. The root depth was obtained following related researches [69][70][71].

Carbon Sequestration Service
Carbon sequestration (CS) is an important regulating service, which is closely related to the productivity and climate regulation capacity of terrestrial ecosystem [72]. The carbon sequestration module of InVEST model was used to calculate the CS of Yantai in 2009 and 2015. The InVEST carbon sequestration module uses land-use data and stocks in four carbon pools (aboveground biomass, belowground biomass, soil and dead organic matter) to estimate the current carbon storage in the landscape. The aboveground biomass includes all plant matter on the soil (such as bark, trunk, branches and leaves). The belowground biomass mainly includes the living root system of the aboveground biomass. The soil organic matter is an organic component of soil and the largest carbon pool on land. The dead organic matter includes vegetation litter and dead wood. In addition, there is a fifth carbon pool, that is, the carbon pool of wood products or forest by-products. As the fifth carbon pool is closely related to the rotation frequency and type of forest products, it is difficult to obtain relevant data in China. Calculation of CS in this study did not involve the fifth carbon pool. The carbon density data used in the model was obtained by referring to the literature related to the study area and providing example data with the model [73,74]. In order to ensure the comparability of the calculation results of CS in 2009 and 2015, this study did not consider changes in carbon density due to changes in vegetation cover. Carbon demand (CD) was obtained by multiplying per capita carbon emission data by population density data. The formula for calculating the CS and CD is as follows: where CS is the sum of carbon storage of all land-use types, Cabove, Cbelow, Csoil and Cdead is the carbon densities of aboveground biomass, belowground biomass, soil organic matter and dead organic matter [75]. Biomass and soil carbon density (Mg/ha) in different land use class are shown in Table 2. S is the area of different land-use types; CD is carbon demand; Dpccc is per capita carbon emission.

Food Supply
Food supply (FS) is an important ecosystem service obtained by human beings from the ecosystem, which plays an important role in ensuring the normal production and life of human beings. Some studies have shown that there is a linear correlation between total grain output and NDVI [76,77]. The total grain output can be allocated to the agricultural grid according to the NDVI value [77]. Specifically, crop yield can be estimated according to the ratio of NDVI value of a grid to the total NDVI value of land type. For example, the yield of grain, oil crops and vegetables in each grid can be estimated for cultivated land, and the yield of melon and fruit in each grid can be estimated for garden land. On this basis, this study can show the FS capacity of different grids, and finally acquire the FS map of Yantai in 2009 and 2015. Food demand (FD) is obtained by multiplying per capita food demand by population density. The formula for calculating the FS and FD is as follows: where FSi is the yield of grain, oil crops, vegetables, melon and fruit assigned by grid i. Gsum is the total production of grain, oil crops, vegetables, melons and fruits in Yantai. NDVIi is the NDVI value of grid i. NDVIsum is the sum of NDVI values for arable land and garden plot in Yantai; FD is food demand; Dpcfc is per capita food demand.

Changes in Ecosystem Services and Land Use
In order to better analyze the spatial changes in the supply and demand of ecosystem services, the difference method was used to explore the changes of various services from 2009 to 2015. The calculation formula is as follows: where ESCi is the change value of the ecosystem services of type i, ESit is the value of the ecosystem services of type i at t time and ESip is the value of the ecosystem services of type i at time p. When ESC > 0, ecosystem services are in a state of increase; when ESC = 0, there is no change in ecosystem services; when ESC < 0, ecosystem services are in a state of decrease. ESCi here only represents the relative increase and decease of services provided by each ecosystem. In the study, the ecosystem services in 2015 were called ESit, and the ecosystem services in 2009 were called ESip.
To clarify how different land-use types were transformed, the study calculated the land-use transfer matrix of Yantai from 2009-2015. Land use transfer matrix can be obtained by overlaying land use maps in 2009 and 2015.

Ecosystem Service Trade-Offs/Synergies Analysis
This study used correlation analysis method to quantify the trade-offs/synergies relationship between different ecosystem services at different times in the study area by referring to the study of Raudsepp-Hearne [78] and Turner [7]. When the correlation coefficient between the two ecosystem services was negative and passed the significance test of 0.05 level, it was considered that there was a significant trade-off relationship between the two ecosystem services. Conversely, if the correlation coefficient passed the significance test and the value of the correlation coefficient was positive, a significant synergistic relationship was considered to exist.
To better spatially analyze the trade-offs/synergies relationship of ecosystem services, this study proposed an ecosystem service trade-offs/synergies index T (−∞< T < +∞), drawing on the research of Pan [25]. The closer T is to 0 indicates the degree of the trade-offs/synergies is weaker, and vice versa. The calculation formula is as follows: where Tij is the trade-offs/synergies index T of the ecosystem services type i and j, and index T was obtained by correlation analysis. ESpi and ESpj are the ecosystem service supply indexes of type i and type j, respectively. A positive value for Tij means the ecosystem services of type i are dominant in the trade-offs/synergies, while a negative value indicates the type j is dominant.

Classification Rule of Changes in Ecosystem Services Demands and Synergies Relationship
Based on Formula (10), this study tried to establish a classification rule to show the changes in three ecosystem services demands, shown as Table 3. During 2009 to 2015, the three demands increased or decreased. Based on the two kinds of variations, 8 types of changing combinations was classified. For example, I-I-I means all of WD, CD and FD increased from 2009 to 2015. Similarly, each synergy was divided into two categories: the ecosystem services of type i or j are dominant, according to formula 11. Then, eight combination types of dominant ecosystem service in the synergies were classified by the rule in Table 4.

Analysis of Land-Use Change
Yantai was mainly composed of arable land, garden plot, woodland and construction land, with the four major categories accounting for 85.5% of the total land area in 2009. In 2015, the spatial distribution of the various land-use types was likely to that of 2009, as shown in Figure 2. During 2009-2015, the total area of arable land only changed 50.85 km 2 , with a change percentage lower than 1%. Garden plot, woodland, grassland and water area also remained relatively stable, all of which have the change percentages lower than 6%. The area of traffic land and construction land increased considerably in 2015, with the main sources of arable land and garden plot. Other land also changed dramatically, with up to 13.31 km 2 and 18.9 km 2 of arable land and water area being converted to it. In terms of the direction of outflow, there was a lot of agricultural land transformed into construction land, and ecological land (including garden plot, woodland, water area and grassland) was mainly converted into traffic and construction land. From the transfer matrix (showing in Table 5), it can also be seen that most of the conversions between different land classes were concentrated in agricultural land-construction land conversion, and ecological landconstruction land conversion, showing an overall pattern of increasing construction land and decreasing agricultural land and ecological land. There was relatively few conversion between agricultural land and ecological land.

Spatial Distribution and Changes of Ecosystem Services
From the spatial distribution, the pattern of WY in 2009 and 2015 was quite different, as shown in Figure 3. In 2009, the areas with high WY were mainly concentrated in the central-eastern part of Yantai. This was mainly related to the high rainfall and elevation in this hilly area. Additionally, there was a weaker land surface evapotranspiration in here than in urban core and developed coastal areas. WY was significantly weaker in the urban core (Zhifu, Fushan, Laishan), western Laizhou and northern Longkou. In 2015, Haiyang and Muping were the areas with the strongest WY in the city. This had a strong correlation with the distribution of precipitation in 2015, which was obviously high in the south and low in the north. WY capacity in the urban core and coastal areas of Yantai was still weak while it in Laizhou was significantly enhanced in 2015. Even though the precipitation of Qixia was low, its WY capacity was still at a medium level in the city with its hilly terrain. The analysis indicated that precipitation and topography were the key factors affecting WY.  The high-value areas of CS in 2009 were mainly concentrated in the inner area of the city. These areas were mostly hilly and had a dense distribution of woodland and garden plot, which can store more carbon compared to other land-use types. There was low CS value in coastal areas and the urban core (Zhifu, Fushan and Laishan). This was mainly due to the fact that these areas were mostly plain, with high population density, developed economy, less ecological land and weak CS capacity. Compared with 2009, CS in 2015 changed little, and most of them changed in the range of [−11 t/hm 2 , 0 t/hm 2 ], as shown in Figure 4. The CS capacity has a slight decrease in most areas of Yantai. Due to the decrease of ecological land, the CS capacity of the region reduced [79][80][81]. In 2009, the high-value areas of FS were mainly concentrated in the south-central part of Yantai. Laizhou and the urban core were the areas with the weakest FS capacity. In 2015, the overall supply capacity was enhanced, but the distribution changed little, as shown in Figure 5. The results of changes in FS from 2009 to 2015 showed that the FS capacity of Yantai was enhanced. The areas with the largest increase were mainly concentrated in Qixia and Laizhou, with a maximum of 10.3 t/hm 2 . Haiyang, Longkou and Penglai also strengthened. There were also some areas in Yantai where the FS capacity was relatively low, and most of the reduced areas were the construction land. Due to economic development, economic activities had to be carried out by encroaching on agricultural land [44,82]. This was consistent with the land-use change that there was a trend of decrease of agricultural land and increase of construction land in most areas of the city.

Analysis of Trade-Offs/Synergies in Ecosystem Services
At the grid, town and county scale, the three pairs of ecosystem services in the Yantai in 2009 and 2015 mostly showed synergistic relationships, as shown in Figure 6. The strongest synergistic relationship was found for WY-CS. Only the FS-CS at the county scale and the WY-CS in 2015 had no trade-off/synergistic relationship. The study found that the aforementioned synergistic relationships had different trends in different scales and years. At the grid and town scale, the synergistic relationship of WY-CS and FS-CS tended to weaken during 2009-2015. However, the synergy of WY-FS was enhanced at grid and town scale. The study also found that most of the three pairs of synergistic relationships had the tendency to strengthen with larger scale. The urban core of Yantai, Laizhou and the northern and southern coastal zone were the areas with the strongest synergy. Except for the above areas, most areas of Yantai had a synergy index between [−2, −0.5], and the synergy between CS and WY was slightly weaker, as shown in Figure 7. These areas were rich of agricultural land and ecological land, with relatively high vegetation coverage. Since surface vegetation had a direct effect on CS and WY, strengthening agricultural land and forest protection could also help increase CS and WY services. In 2015, the synergetic relationship between CS and WY changed greatly. Areas dominated by WY decreased while areas dominated by CS increased. In areas dominated by WY, the intensity was significantly weaker than in 2009. The WY-oriented areas in the central and northern part of Yantai were significantly reduced. The urban core and Penglai with a synergy index less than −2 in 2009 changed into the range of [−0.5, 0], indicating the increase of synergy between WY and CS in these areas.
In 2009, the synergy between FS and WY was highly dominated by WY, with the synergy index in the range of [−4, −3] in most regions. In 2015, the synergy between FS and WY increased but the leading role of WY was weakened, as shown in Figure 8.    The synergy between FS and CS was likely to the previous two. The three pairs of ecosystem services were all dominated by one of the ecosystem services. In 2009, the synergy was dominated by CS. The areas with the strongest leading role of CS were mostly concentrated in Laizhou, Penglai and Haiyang, with a synergy index less than −5. In the vast inland areas of Yantai, including Qixia and Mouping, CS dominated as well. In 2015, the dominant ecosystem service was still CS, but the intensity decreased, as shown in Figure 9. The area where the synergy index was between [−3, −2] was reduced, while [−2, −1] increased significantly.

Adjustment and Improvement of Ecosystem Services Based on Spatial Imbalance in Supply and Demand of Ecosystem Services
As shown in Figure 10, WD, CD and FD changed differently. Increased WD and CD took up more areas than decreased one, but FD changed in the opposite. As shown in Figure 11, the city was divided into four demands changing areas: I-I-I, I-I-D, D-I-D and D-D-D, according to the classification rule in Table 3. I-I-I, which represents the increase of three ecosystem services demands, accounts for the largest proportion of 35.99% in the whole city. I-I-D means increase in WD and CD but decrease in FD, was least in the four. There only was an increase in CD in D-I-D. Additionally, in D-D-D, all the three demands decreased. It was urgent to reinforce the ecosystem services in I-I-I area. Additionally, in this area, the combination types of dominant ecosystem service in the synergies were shown in Figure 12 and Table 6. Type I, in which dominant relationship in three synergies was WY > CS, WY > FS, CS > FS, constituted the largest proportion of 83.58%. In this type, FS was in the most vulnerable position. Thus, FS improvement would be given the priority, followed by CS, and finally WY. In type II, FS was relatively weak as well. With the continuously upgrading grain consumption structure, as seen in Figure 10, FD increased in a small part of the city, and it was all in I-I-I area. That meant over 98% of I-I-I area was to take FS enhancement as primary task, to meet the increased FD.  In I-I-D area, WD and CD increased, FD decreased. Thus, more attention should be paid to the synergy of WY and CS. As shown in Figure 13, 85.74% of the I-I-D region was dominated by WY. Priority should be given to develop CS in WY-dominated area. In D-I-D, only CD increased. Then CS improvement was in the first place. In D-D-D, three ecosystem services demands reduced, and the improvement priority of different types was consistent with in I-I-I. However, D-D-D region provided ecosystem services to other areas in need, as shown in Table 7. Specifically, it provided: WY in type I and IV to type III and VI in I-I-I region, and CS-dominated area in I-I-D region; CS in type II and III to type IV and V in I-I-I region, WY-dominated area in I-I-D region and the whole D-I-D region; FS in type V and VI to type I and II in I-I-I region.
Based on above, priority in ecosystem service improvement in Yantai city was classified into three zones, as shown in Figure 14. FS and CS priority improvement zone composed 98.66% of the whole land, and FS priority improvement zone covered over 60%. WY-CS-FS showed an interactive and supportive relationship in the city [83], and WY has always been dominated. To some extent, development of CS and FS were based on WY. In guaranteeing the water balance between supply and demand, it was FS and CS that in the first place to improve. However, this did not mean stop of WY development. For the government, with the exception of development of weak ecosystem service, policy tradeoff in different ecosystem services and resource scheduling were key issues. It is important to dispatch ecosystem service from oversupply area to demand dissatisfaction area.

Further Analysis of Ecosystem Services Caused by Land-Use Change
Land-use change affects the structure and process of ecosystem services, which in turn affects the ecosystem services trade-offs/synergies [8,10,13,39]. Land-use change in Yantai was relatively drastic, with an overall trend of decreasing agricultural land and ecological land and increasing construction land. A large amount of ecological land would inevitably be converted to construction land in the future development process, destroying the environment and thus weakening the ecosystem services supply [44,45]. For example, the expansion of construction land can lead to a reduction in the vegetation cover, thereby reducing carbon sequestration services. It can also be seen from Figure 3 that the reduction of forest in the urban core (Zhifu, Fushan and Laishan) reduced the carbon sequestration capacity more significantly. In the following-up study, we should focus on the change of ecosystem services interrelationships due to land-use change, and use scenario analysis to explore the change of ecosystem services trade-offs/synergies in different development scenarios to achieve sustainable development of the region [31].

Driver Exploration of the Synergies in Ecosystem Services
This study found that there was a clear scale effect in ecosystem services synergies. The characteristics of synergies differed in regions and scales. Synergistic relationships had the tendency to strengthen with larger scale. This may be due to the fact that the grid scale was only representative of the ecosystem service provisioning capacity in a small area, and the differences between the grid scales were relatively large. As the scale increased, the ecosystem service values in the town and county were more comprehensive, which can better reflect the regional ecosystem service provisioning capacity. The results of this study were similar to those of Zhang [84] on Beijing-Tianjin-Hebei region, which showed that the trade-offs/synergies between multiple pairs of ecosystem services in Beijing-Tianjin-Hebei region also tend to strengthen as the scale increases (from 1 km grid scale to county and city scale). However, the study also noted that this trend was not very clear in Qin-Ba Mountains [58] and the farming-pastoral ecotone [85]. This may be due to the spatial heterogeneity of natural elements and the differences in the demand of different stakeholders. Therefore, there is a need to explore the underlying causes of changes in the relationship between ecosystem services in different areas in order to better formulate locally appropriate policies for ecosystem management.
In this study, the relationship between CS and WY in Yantai was synergistic, which is same to in China's national ecological barrier zone [86] but opposed to in Northwest China [87]. Additionally, from 2009 to 2015, the synergy changed significantly, which was related to the decrease in WY caused by the lack of precipitation in Yantai [88]. In addition, the weakening of the WY may also be related to the large-scale implementation of ecological projects such as the greening of barren hills, the greening of villages and the Grain-to-Green. According to the government work report of Yantai, the area of woodland increased from 534,600 hectares to 546,500 hectares from 2009 to 2015. Studies have shown that forest cover can be improved to some extent through afforestation, while also playing a positive role in ecosystem services such as CS [89][90][91]. On the other hand, the evapotranspiration of the region was enhanced to some extent and WY was weakened through afforestation [92][93][94]. This was an important reason for the increase in the area with CS as the dominant ecological service in 2015.
Synergy between FS and WY was dominated by WY in 2009 and 2015. This may related with photosynthesis. WY enables plant tissues to absorb water for photosynthesis, and the accumulation of crop yields mainly depends on photosynthesis, which converts solar energy into chemical energy by green plants. Meanwhile, soil moisture is directly related to plant root water absorption and leaf transpiration, and affects the dry matter accumulation, and finally the grain yield [95,96]. Luo [97] and Zhang [98] also found that increasing soil moisture content improved the photosynthetic properties of leaves and increased the biomass and grain yield of wheat. The reduction in rainfall was one of the important factors [88]. The synergy between FS and WY in the southern area of Yantai changed little in 2015, and remained within the range of [−4, −3], which was related to the fact that the region had more precipitation.
There was also a synergistic relationship between FS and CS, and the synergy was dominated by CS. This was again similar to the findings at other scales [32,97]. In the growth and development, the plants rely primarily on photosynthesis to absorb carbon dioxide, which undergoes a series of complex chemical reactions to produce carbohydrates that are stored by vegetation and ultimately provided directly to humans in the form of food [98]. This is an important reason why the synergistic was dominated by CS. From 2009 to 2015, reduced precipitation was unfavorable to vegetation growth, resulting in weakened CS services. This may account for the weaken of CS's leading role.

Analysis of the Changes in Ecosystem Services Demand and Improvement of Ecosystem Services
The study also found that from 2009 to 2015, WD and CD increased in most areas but FD was on the opposite. It was closely related to per capita demands of each ecosystem service. With the improving living standards, there was a decrease in FD but an increase in meat, eggs, milk and other food, calling for a more reasonable dietary structure. For those little areas of FD increased, on the one hand local FS should be strengthened first by improving food production technology and scientific planting and breeding ability [90]. It was also important to give more policy support and subsidy for food production in this area was. On the other hand, it was a good way of making better use of surplus FS in most areas of the city. Additionally, it required a construction of food transportation corridors and distribution stations [51,56,90]. However, in this dietary structure, CD increased a lot, especially in more developed area [99]. It was because that consumption of meat, eggs and dairy products produced higher carbon emissions than food consumption did. Conservation and restoration of forest was the key to increase CS. Furthermore, it should strengthen mine recovery and urban greening. To reach carbon equilibrium, it also mattered that people should try to reduce carbon emission to decrease the CD, by energy structure adjustment and improvement [90]. As one of the typical water shortage cities, WD increased a lot in Yantai, owing to industry production, coal mining, vegetable and fruits planting and living [57], however. A greater imbalance between water supply and demand was seen. It was urgent to develop WY. It was hard to regulate precipitation so that water reuse and construction of sponge city may be good solution for it.
Based on changes in ecosystem services demand and ecosystem services synergies, three prior zones of ecosystem service improvement were divided: FS-prior, CS-prior and WY-prior. Rapid urbanization makes the neglect of agriculture planting and ecological land protection, policy support and subsidy from government is the key for FS-prior area. The government should ensure that the red line for protecting cultivated land is not crossed and strictly control land uses. On this basis, it is important to improve land quality and grain yield. In CS-prior area, the most important thing is to strengthen the protection and restoration of forestry resources. For Yantai, the government should try to strengthen afforestation projects in barren mountains and the construction of key coastal shelterbelts, and implement the "Grain for Green Programs" in sloping cultivated land above 25 degrees, key sandy land and important water sources. In addition, the management of newly built forest land and the tending of young and middle-aged forests are important as well. For WYprior area, the first one is to ensure the fully protection of water source areas. Furthermore, water purification projects can be constructed at river confluence reaches, river estuaries and other suitable places. Additionally, within the flood control embankments of rivers and lakes, the programs of conversion of farmland to wetland and fishing to lakes can be carried out.

Limitation and Further Research
It should be mentioned that due to the time-cumulative nature of ecosystem services, the variability of different types of services varied greatly, which affected the spatiotemporal expression of ecosystem service relationships [100]. The relationships between ecosystem services showed stages and differences over time. Due to the difficulty of acquiring landchange survey data, this study was limited to a short spatiotemporal scale and may not be able to fully clarify the interrelationships among ecosystem services. Moreover, this study focused on the relationship between provisioning and regulating services, with only three ecosystem services selected. Quantitative analysis of the mechanisms driving trade-offs/synergies between ecosystem services was still lacking. There was still a need to conduct more detailed and long time series of ecosystem service assessments, form ecosystem service bundles, study the trade-offs and synergies in the bundles and quantitatively investigate the underlying causes of changes in ecosystem services at different scales on the basis of macroscopic studies, to propose more precise land management policies according to local conditions [77,101].
Meanwhile, there was a limit in the characterization of supply and demand. In the calculation of WY, CS and FS, most of the parameters used were based on relevant research results and manual recommendations. InVEST model was applied in the estimation of CS. InVEST has the advantages of dynamic, spatial, multi-level and multi-module. It is highly scientific, intuitive and operable and widely used in different types of ecosystem service assessment at home and abroad [69,[72][73][74][102][103][104][105]. CASA, i-Tree and some other model are also used for CS estimation. However, in China, few studies use i-Tree model [106], owing to the difficulty of acquiring the data of species composition, diameter/size distribution, species diversity, tree cover and so on. In further research, there was a need to collect the data for i-Tree model and make a comparison between the results, through field research and environmental monitoring. The demands in the study were calculated by per capita data demands. However, the quantification of ecosystem services demand is influenced by cultural and socio-environmental factors such as gender, age, marital status, education, income and ethnicity, and may change over time and spatial scales [107][108][109]. The trade-offs between different ecosystem services can then also lead to trade-offs between the wellbeing of different populations [110]. Meanwhile, the impact of external forces and global markets (including corruption and governance) on ecosystem services tradeoffs/synergies and human well-being cannot be ignored [111]. For example, agricultural management is influenced primarily by ecosystems related to food supply, but also by biophysical and socio-economic changes and management practices, as well as by access to markets and patterns of trade [112]. The complexity of ecological, social, natural and economic factors also exacerbates the difficulty of studying the trade-offs/synergies between ecosystem services and human well-being. Future study should further enhance interdisciplinary synergy and collaboration from ecology, sociology, geography and economics. Quantitative adjustment and cross-regional flow need further consideration in improvement, though with different regional interests, cross-regional ecosystem conservation and management was difficult to realize [113]. To gradually overcome the problem and made an overall scheduling of cross-regional ecosystem services, the government need to play its leading role.

Conclusions
Human beings meet and improve well-being by consuming a variety of ecosystem services and products. However, the types and quantity of ecosystem services required by stakeholders differ in time and space. The spatial mismatch between supply and demand of ecosystem services directly leads to environmental injustice and affects human well-being seriously. This study analyzed temporal and spatial changes in supply and demand of three ecosystem services (food supply, carbon sequestration and water production) in Yantai, a typical city in Bohai Bay, China. Trade-off/synergy relationship and driving mechanism of the three types of ecosystem services were revealed. Classification rule of changes in ecosystem services demands and dominant ecosystem service in the synergies were put forward to acquire supply-demand imbalance areas of ecosystem service and provide improvement suggestions.
The results showed that from 2009 to 2015, the conversion of different land-use types in Yantai was mainly concentrated in agricultural land-construction land and ecological land-construction land, showing a pattern of increasing construction land and decreasing agricultural land and ecological land. Ecosystem services differed significantly in Yantai. The high-value areas of WY, FS and CS were mainly concentrated in the hilly area of south-central Yantai, while the urban and industrial areas along the north coast were obviously weak. From 2009-2015, the capacity of FS increased in most areas while WY decreased in most areas, and CS changed little. In 2009 and 2015, Yantai showed synergistic relationships among the three pairs of ecosystem services, with the strongest synergistic relationship for the FS-WY. Most of the three pairs of synergistic relationships had the tendency to strengthen with larger scales. The synergistic relationships of CS-WY and FS-WY were dominated by WY, and FS-CS was dominated by CS. According to the four regions of changes in ecosystem services demands and three trades-off/synergies in each region, improvement suggestions of ecosystem services were proposed.