Spatiotemporal assessment and trade-offs ofmultiple ecosystem services based on land use changes in Zengcheng, China

• We put forward the general structure of assessing comprehensive ecosystem service. • Trade-offs among various ecosystem services was analyzed. • Four alternative scenarios were designed to identify how to improve the comprehensive ecosystem service. • More sustainable intensification of agriculture should be adopted to increase the food production. ⁎ Corresponding author. E-mail address: lifeng@rcees.ac.cn (F. Li). http://dx.doi.org/10.1016/j.scitotenv.2017.07.221 0048-9697/© 2017 Elsevier B.V. All rights reserved. a b s t r a c t a r t i c l e i n f o


Introduction
In recent years, ecosystem services (ESs) have increasingly gained attention worldwide (Fisher et al., 2009;Kragt and Robertson, 2014). ESs represent the human benefits that have been formed and sustained by ecosystems and ecological processes. The decrease of ESs seriously threatened human health and survival (MA, 2005). Conversely, human activities also affect ecosystem processes and components (Burkhard et al., 2012). Nearly two thirds of the ecosystems around the world have been destroyed by humans (MA, 2005). Urban areas represent the densest concentrations of human activities (Wang et al., 2011). In the process of urbanization, human activities have a great influence on land use changes, which would affected the ESs in reverse (de Groot et al., 2010;Liu et al., 2014). The changes of urban structure and function have led to a series of problems such as decreased carbon storage (Jiang et al., 2017), degraded water quality (Gómez-Baggethun and Barton, 2013), serious soil erosion (Yao et al., 2016) and decreased biodiversity (Maes et al., 2012).
The assessments of ESs can provide us quantitative information about ESs changes in the urbanization (Robinson et al., 2013). Most of studies mainly examined the natural ecosystems (Delgado-Aguilar et al., 2017;Tian et al., 2016;B. Zhang et al., 2017;X. Zhang et al., 2017) and only focused on one or a few ESs (Guo et al., 2013;Schirpke et al., 2016;Yao et al., 2016). However, few studies examined the temporal and spatial changes of overall ES in urban areas . Until now, land use change analysis has become a fundamental tool for assessing the impacts of human activities on ESs in urban areas (Du et al., 2009). Study demonstrated that woodlands, wetlands, water bodies and orchards in the urban areas play the most important role in ESs. The combined ES value of these land use categories accounts for over 90% of the total value (Li et al., 2010). Moreover, the rapid expansion of urban and industrial land in urban areas has led to lower food supply, carbon sequestration, soil water storage, and habitat suitability . To determine the importance of land use changes on ESs, many studies have not only evaluated the impact of actual land use changes, but also compared a suite of alternative land use scenarios (Geneletti, 2013;Zhu et al., 2011). Scenario analysis was used to predict the future impact of land use change on ESs by establishing different alternatives associated with the future implementation of plans and policies; this can help achieve a variety of goals and can be used to inform policy-maker and stakeholder decisions about the most optimal ES provisions (Polasky et al., 2011;Geneletti, 2012). These scenarios have been employed in both theoretical frameworks (Seppelt et al., 2013) and empirical applications (Butler et al., 2013). The alternative scenarios in most studies were designed according to local land use polices (Nelson et al., 2009;Polasky et al., 2011) or only focused on improving one specific ES (Chuai et al., 2015). A few studies designed the scenarios by considering both local land use polices and the trade-off analysis of multiple ESs. Trade-off analysis could help us explore the optimal land use pattern which can improve the specific ESs or overall ES (Feng et al., 2017).
Quantitative modeling also plays a vital role in assessing ESs. Current evaluation methods range from simple spreadsheet models to complex software packages (Bagstad et al., 2013), including biophysical, empirical, GIS-based, decision support and survey-based models and less widely used approaches such as lookup tables (Haase et al., 2014). Biophysical models have mainly focused on the value of forests in reducing air pollution (Jim and Chen, 2009). Most of the empirical methods were used to analyze services regulating green infrastructure (Haase et al., 2014). GIS-based models, such as the Social Values for Ecosystem Services (SolVES) model and the Recreation Opportunity Spectrum (ROS) model, have been widely applied to assess ES provisions and can reveal the heterogeneity of spatial distributions (Angold et al., 2006). SolVES is an ArcGIS toolbar for mapping the social values of ESs based on survey data or transfer values. However, this model is limited when the study site is drastically different from past study sites (Sherrouse et al., 2011). ROS model was designed to provide a framework for outdoor recreation management and can map recreation opportunity indices in ArcGIS. It was developed in the United States and has been applied in different countries, such as the European Union (Paracchini et al., 2014) and China . Many decision support tools have been developed for integrating ESs into both public and private sector decision-making processes (Bagstad et al., 2013), such as the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) and the Artificial Intelligence for Ecosystem Services (ARIES) models. Among them, the InVEST model has been widely used because it can be independently tested and the software is open sourced. Both ES spatial distribution and trade-off maps can be presented by this tool (Grafius et al., 2016). The ARIES model is a web-accessible application that can build and run ES spatial flow analysis in a given area based on user goals; however, this new tool is currently time consuming and cannot be generally used until global models are completed (Villa et al., 2011). In addition, survey-based models were designed to improve the understanding of both how human behavior affects ES provisions and how people respond to ESs; they have been mainly used to estimate the recreation potential of green areas (Maas et al., 2006). Look-up tables have been generally used to transfer results from previous studies to current studies of interest (Kreuter et al., 2001).
Because of the rapid economic development, urbanization has occurred at an unprecedented rate in China since the beginning of the 21st century (Tao et al., 2016), especially in a few large, highly developed cities, such as Beijing, Guangzhou and Shanghai. Compared with the ESs of other cities in China, the ESs of these megacities were more negatively affected by land use change during urbanization. Therefore, assessing the spatiotemporal effects of urbanization on multiple ESs in these megacities will help guide urban planning and formulate land policies in China (He et al., 2016). Zengcheng, Guangzhou was chosen as a typical study area to quantitatively assess the ESs, including climate regulation, water supply, water purification, soil conservation, biodiversity protection, recreation opportunity and food production, in this study. In addition, this study reports a new index for measuring the comprehensive ecosystem service (CES). Different alternative scenarios were designed according to local policies and trade-off analysis of the ESs. The main objectives of this study were to: (1) integrate different models to assess the CES in Zengcheng; (2) put forward some ecological engineering and land management measures for city planners and policy makers to improve the corresponding ESs and achieve sustainable urbanization in Zengcheng; (3) explore the optimal alternative land use pattern, which can provide a higher CES and simultaneously minimize the trade-offs of different ESs.
Zengcheng has abundant natural resources and a mature industrial system. Local government has been actively exploring sustainable land use planning policies in recent years. The northern region including PT, XL, ZG and northern ZX, is mainly covered by mountains and forest land. This region has focused on the development of agriculture and ecological tourism. However, the environmental problems, such as biodiversity threats, river pollution and a rapid increase in energy consumption, have emerged because of rapid tourism development. Therefore, economic activities that potentially affect the environment in this region are now strictly restricted by the local government. The central region is mainly hilly and includes ZC, LC and ZJ. This region was designed to promote urban culture services and residents' quality of life. Educational institutions and leisure centers in Zengcheng are mainly distributed in this region. The southern region, which is near downtown Guangzhou, was designed to be the major urbanization development area. Industrial lands expanded rapidly in this region, especially in XT. Local governments devoted their efforts to developing denim jeans and the automobile industry. Available land resources have become very scarce in this region because the population density is increasing. Additionally, industrial activities have created substantial pressure on the ESs in the southern region (ZMSB, 2013).

Land use classification in Zengcheng
In this study, we used the raster maps generated in 2003, 2008 and 2013 at a 30-m spatial resolution to detect past land use changes in Zengcheng during urbanization. The raster maps were converted by shapefile data, which were drawn by the local planning bureau according to investigational information gained from actual land use in Zengcheng. Land use types in this study were classified into seven categories according to the national standard for classification of land use (http://www.mlr.gov.cn/): woodland, grassland, cropland, water bodies, built-up land, traffic land and undeveloped land (Table A.1). Fig. 2 shows the general structure for assessing the CES in Zengcheng. We used a subset of the InVEST model and focused on reporting the services for climate regulation, water supply, water purification, soil conservation and biodiversity protection in biophysical terms. The InVEST model estimates the ESs according to the given land use maps and biophysical data for the region ). Additionally, the Recreation Opportunity Spectrum (ROS) model was used to map the recreation opportunity service, and food production was predicted by the yield model. The data sources for different parameters in the models were shown in Table A.2.

Climate regulation
The InVEST model maps and quantifies the amount of carbon stored and sequestered based on four carbon pools: above ground biomass, below ground biomass, soil and dead organic matter. The biomass carbon stocks per unit area of each land use type assessed in this study were derived from the results of local studies (Fang et al., 2010;Guo et al., 2013;Ni, 2004) (Table A.3).

Water supply
Water supply in the InVEST model is defined as the amount of water runoff from the landscape (Tallis et al., 2011). Water yield is based on the principle of water balance, and average annual precipitation (P x ) and actual annual evapotranspiration (AET x ) are used to calculate the annual water yield (Y x ) in each 30 m × 30 m grid cell according to the following equation: where AET x is the actual annual evapotranspiration for pixel x and P x is the average annual precipitation in pixel x.

Water purification
Water purification service refers to the capacity of ecosystems to mitigate water pollution by retaining some non-point source pollutants in vegetation and soil (Tallis et al., 2011). Higher nitrogen and phosphorus export values for each LULC reflect lower water purification capacity. This module estimates the quantity of pollutants retained for water purification through a three-step process. Average amount of nutrients exported annually from each grid cell can be calculated using Eqs.
(2)- (4): where ALV x is the adjusted loading value at pixel x, pol x is the export coefficient at pixel x, and HSS x is the hydrological sensitivity score at pixel x. Next, HSS x can be calculated as: where λ x is the runoff index at pixel x, λ w is the mean runoff index in the watershed, and λ x can be calculated using the following equation: where ∑ U Y u is the sum of water pixel yield along the flow path above pixel x. We modified literature reports (Bai et al., 2013;Fu et al., 2014;Huang, 2014) and consulted with experts to represent the export coefficients of this region (Table A.4).

Soil conservation
Sediment retention module is used to calculate the average annual soil loss from each parcel of land, and reduced soil loss indicates improved soil conservation. This module employs the Universal Soil Loss Equation (USLE) (Wischmeier and Smith, 1978) on a grid cell scale to predict the average annual soil loss rate: where R x is the rainfall erosivity for pixel x, K x is the soil erodibility factor for pixel x, LSx is the slope length-gradient factor, C x is the crop management factor for pixel x and P x is the support practice factor for pixel x. Soil erodibility factor was estimated by using the EPIC model (Williams et al., 1983). Crop management and support practice factors for each land use type were modified according to the relevant literature (Hamel et al., 2015;Yang et al., 2003;Yao et al., 2016) (Table A.5).

Biodiversity protection
Habitat quality refers to the ability of an ecosystem to provide suitable conditions for the persistence of individuals and populations (Mckinney, 2002). Habitat quality depends on the relative impact of threats, the sensitivity of habitats to threats, the distance between habitats and the threat sources and the location of the protected areas. The impact i rxy of threat r from grid cell y on the habitat in grid cell x can be represented using the following equations: where d xy is the linear distance between grid cells x and y and d rmax is the maximum effective distance of the threat. In this study, the threats to croplands and built-up land followed exponential decay and the threat to traffic land followed linear decay. The total threat level D xj in a grid cell x with land use/land cover j is calculated as The habitat quality Q xj of land use/land cover j is finally calculated as In this model, the sensitivity of different threat sources for land use/ land cover was based on previous studies (Leh et al., 2013;Terrado et al., 2016).

Recreation opportunity
Outdoor recreation service in this study was provided by green spaces and other natural landscapes. A flowchart of the procedure for the ROS model is shown in Fig. 3. This model was associated with both the recreation potential index and the zoning of Zengcheng in terms of the remoteness and accessibility of reaching recreation sites. The remoteness/accessibility index was obtained by calculating the distance from urban areas and roads in grid cells using the ArcGIS toolbox.

Food production
Zengcheng is considered to mainly produce three types of food: grains, vegetables and fruits (http://nyj.zengcheng.gov.cn/). Food production can be calculated using the yield model: where PRO G , PRO V and PRO F represent the production of grains, vegetables and fruits, respectively. A i is the area of town i in Zengcheng; R Gi , R Vi and R Fi represent the area proportions of grains, vegetables and fruits in town i, respectively; and P Gi , P Vi and P Fi represent the per unit area yield of grains, vegetables and fruits for each town, respectively.

CES index
A comprehensive index for ESs was developed in this study to compare the CES level. This new model weighs various values that constitute the total index value for the CES. The model is expressed as follows: where ES j is the total index value for the CES in the year or scenario j, w i is the weight assigned to the ith ES by stakeholders and experts, S ij is the standardized value for the ith ES in the year or scenario j, and n is the number of ESs being evaluated. In addition, w i is the relative importance of each individual ES to various stakeholders. In this study, the Analytical Hierarchical Process (AHP) technique (Zhang and Lu, 2010) was used to derive the weight for each ES. AHP is always used to make complex decisions that involve multiple criteria; its structure comprises the goal, criteria and alternative levels. The weights for the different ESs in Zengcheng (Table 1) shows that nutrient export, sediment export and recreation opportunity account for higher proportions of weights. For Zengcheng, the main environment problems were water pollution in the southern region and soil loss in the northern region. Compared with other ESs, these two services played more important role in urban sustainable development. In addition, with the increase of people's quality of life, tourism and leisure has become a more popular lifestyle. Therefore, recreation opportunity is also a very important service in Zengcheng. To ensure that different ES indices could be added, we standardized the values of these services (see Supplementary data). Fig. 4 shows the spatial distribution of land use in Zengcheng from 2003 to 2013. Woodland and cropland were the two largest land use types, while grasslands and undeveloped land were accounted for the smallest proportion in Zengcheng from 2003 to 2013. The northern area was mainly covered by cropland and woodland. The western central district was mainly covered by cropland. Downtown was in the eastern central part, which was mainly covered by built-up land. The southern area, especially the southwest district, was highly urbanized. Built-up land and traffic land accounted for the largest proportion in this region.

Land-use changes from 2003 to 2013 in Zengcheng
For Zengcheng, the woodland, grassland, built-up land and traffic land areas all increased from 2003 to 2013. Among them, built-up land and traffic land expanded the most, as they expanded by 106% and 129%, respectively. The water bodies, cropland and undeveloped land areas all decreased during this same period, especially cropland, which decreased the most. Table 2 shows the land use changes in all the towns from 2003 to 2013. Built-up land and traffic land in all towns increased a lot, while the cropland in most of towns, except in ZC, all decreased during the past ten years. Among the towns, XT exhibited the greatest decrease in cropland and the greatest increase in builtup land and traffic land. For the woodland area, LC and ZX experienced the largest increase, while XL exhibited the largest decrease compared with the areas in the other towns. In addition, except for in XT and PT, the grassland and undeveloped land areas in most of the towns had no obvious changes. Regarding the water bodies, LC experienced the largest increase, while XT exhibited the largest decrease compared with the area of water bodies in other towns.   Table 3 shows the indicator values for various ESs in Zengcheng from 2003 to 2013. Carbon storage was decreased by 2% during the past ten years, while water supply increased from 2003 to 2008 and then decreased from 2008 to 2013, mainly as a result of changes in precipitation. Water purification was improved because of decreases in nitrogen export (−6%) and phosphorus export (−11%). Soil conservation became worse as sediment export increased by 25%. Biodiversity protection and recreation opportunity were decreased by 4% and 6%, respectively. Regarding food production, both vegetable and fruit production increased by 33% and 68%, respectively, while grain production decreased by 24% from 2003 to 2013. Fig. 5 shows the ESs in each town as a percentage of the entire area from 2003 to 2013. Carbon storage, water supply, biodiversity and recreation services were mainly provided by the northern region, especially PT and ZG. XT in the southern region exhibited the largest proportion of nitrogen and phosphorus exports because of the heavy pollutant loads on impervious surfaces and weak purification ability caused by low vegetation coverage. Largest proportion of soil loss occurred in the northern region, especially in PT. In addition, food production showed different distribution patterns. Grain production was mainly provided by the northern region, especially PT, and vegetable production was mainly provided by the southern region, especially ST. However, fruit production was mainly provided by ZX in the central area.
The supply of multiple ESs per unit area in each town from 2003 to 2013 is shown in Fig. 6. Carbon storage per unit area followed a pattern similar to that of vegetation coverage, as its value decreased from north to south. From 2003 to 2013, the water yield change was consistent with precipitation, which increased at first and then decreased on the whole. Although evapotranspiration in the north was higher because of its higher vegetation coverage, water supply in this area was still higher because of the larger amounts of precipitation. Water purification was higher in the north than in the south, as nitrogen and phosphorus exports per unit area were lower in the north. This is because the proportion of cropland was higher in the southern region and constant agricultural activity led to increased pollution exports. Additionally, the higher proportion of built-up land and traffic land produced higher nitrogen and phosphorus loads than the other land use types in the southern region. However, the northern region suffered a greater soil loss, as sediment export per unit area was higher in the north than in the other regions. This is mainly because the steeper slope and greater precipitation in the north had resulted in higher rainfall erosion and more soil loss. Soil erosion was most serious in 2008 because of the heavy rainfall in this year. Biodiversity index in the north was higher than that in the south because the northern region had higher vegetation coverage and lower degree of habitat fragmentation. Recreation service was also decreased from north to south because natural landscape was mainly distributed in the northern area. Grain, vegetable and especially fruit production per unit area were all improved during the past ten years. The grain production per unit area was not obviously different across the various towns. The food productions were affected by many factors, such as climate, soil properties, agricultural technology and management practice. XL exhibited the largest vegetable yield per unit area, while PT and ZX exhibited a larger fruit yield than the other towns.  In this study, we standardized the ESs by first using Eqs. (1) and (2) Fig. 7 shows the spatial distributions of CES in Zengcheng. Higher CES values were mainly distributed in the northern region, while lower CES values were observed in XT in the southwest and LC in the central region. During the past ten years, the CES value increased in XT mainly because the cropland area decreased and the vegetation coverage increased, which improved the water purification. The region north of ZG exhibited an increase in CES value mainly because the undeveloped land area decreased and the vegetation coverage increased, which reduced soil loss and improved the recreation service in this area. However, increased intensive traffic network in all the towns resulted in a decreased CES value. Moreover, the   central towns, especially LC and ZC, exhibited significant declines in CES because rapid expansion of built-up land had led to the decreases in carbon storage, biodiversity and recreation opportunity.

Trade-offs among multiple ESs in Zengcheng from 2003 to 2013
Correlation analysis can be used to help determine the relationships between different ESs (Wu et al., 2013). We visualized interactions by plotting the correlation between pairs of ESs in 2013 at the grid level. This produced both a correlation coefficient (r-value) and a significant value (P-value) ( Table 4). The results showed significant (P b 0.01) positive correlations among carbon storage, water yield, water purification, biodiversity protection, and recreation opportunity. The areas that had higher vegetation coverage and were less urbanized exhibited higher carbon storage, water purification, biodiversity protection, and recreation opportunity values. Conversely, soil conservation had significant (P b 0.01) trade-off relationships with other ESs except food production. This is because most high ES values in Zengcheng were mainly distributed in the northern mountains. However, the steeper slope in these areas resulted in greater soil loss. Additionally, food production had trade-off relationships with carbon storage, water purification, biodiversity protection, and recreation opportunity, indicating that that areas with higher provisioning service (food production) cannot simultaneously provide higher regulating services (carbon storage and water purification), supporting service (biodiversity protection) and cultural service (recreation opportunity).

Trade-offs among various alternative scenarios in Zengcheng
In this study, land use in 2013 was chosen as the base case scenario (baseline). According to correlation analysis, climate regulation, water supply, water purification, biodiversity protection, and recreation opportunity have synergetic relationships. To improve these services, we designed a woodland buffer (WB) scenario because the government in Zengcheng has been working on an "Ecological corridor construction" project, which includes developing riparian vegetation zones and highway green belts (OWZM, 2015). Correlation analysis also showed that the high ES values in Zengcheng were most often accompanied by great soil loss. To reduce soil loss, we designed a soil conservation (SC) scenario modeling that of the "Grain for Green" policy in China. This policy mainly aims to restore vegetation on previously steep croplands (Wang et al., 2011). In addition, we need to increase food production when improving other ESs. An agricultural expansion (AE) scenario then was developed to increase the food production. Finally, we designed a combined development (CD) scenario to balance various ESs according to above measures. In total, we designed four alternative scenarios to explore ES improvement. Table 5 shows the descriptions of four alternative scenarios and their associated land cover changes.
Our assessment results showed that average CES values under both the WB (0.49) and SC (0.49) scenarios perform better than that under the base case scenario. However, the AE (0.43) and CD (0.46) scenario CES values perform worse than that under the base case scenario. Fig. 8 shows the change rates of the various ESs under alternative scenarios relative to the base case scenario. WB and SC scenarios perform better than the base case scenario in terms of carbon storage, water purification, soil conservation, biodiversity protection, and recreation opportunity. Among them, greater improvements in carbon storage, water purification, biodiversity protection, and recreation opportunity are mainly located in the central and southern regions, while greater  Cropland and undeveloped land within 100 m surrounding large rivers and reservoirs were converted to woodland buffer strips. Simultaneously, the cropland and undeveloped land within 50 m of each side of highways and railways were converted to woodland buffer strips. The width of the woodland buffer strips was set based on the actual land usage in Zengcheng and practical experiences in other Chinese districts (Yang et al., 2003;Zhang et al., 2002).

Soil conservation (SC) scenario
The croplands and undeveloped lands with a slope greater than 15°were converted to woodlands. Previous research has shown that steeper slopes equate to easier, more intense soil erosion, especially when the slope is greater than 15° (Bai et al., 2013;Sun et al., 2015).

Agricultural expansion (AE) scenario
The woodlands and undeveloped lands with a slope less than 6°were transformed to cropland because studies have indicated that these flatlands are suitable for farming (Bai et al., 2013;Ge et al., 2012). Under this scenario, we did not change the woodlands and undeveloped lands in urban centers, squares, parks, scenic areas or nature reserve.

Combined development (CD) scenario
The woodland buffer strip setting was the same as in the woodland buffer scenario. As in the soil conservation scenario, the croplands and undeveloped lands with a slope greater than 15°were converted to woodlands. Moreover, as in the agricultural expansion scenario, the woodlands and undeveloped lands were transformed to croplands if the slope was b6°.
improvements in soil conservation are mainly distributed in the northern region. However, these two alternative scenarios perform worse than the base case scenario in terms of food production. AE and CD scenarios perform better than the base case scenario in terms of food production but perform worse in terms of the other ESs. Most of the services showed larger decreases in the north than in the south under both the AE and CD scenarios. Considering all the alternative scenarios, WB scenario increased the most in recreation opportunity (16%), while SC scenario exhibited the greatest improvement in soil conservation (8%). AE scenario exhibited the largest increases in grain (22%), vegetable (22%) and fruit production (37%) but performs worst in terms of the other ESs, except water supply. In addition, although various ESs under the CD scenario were more balanced than those under the other scenarios, the CES value under this scenario was reduced compared with that under the base case scenario.

Linking land use strategies and ecological engineering measures with multiple ESs
Optimal land management and urban planning require joint consideration of values of all the objectives (Polasky et al., 2011). Integrating ESs into urban planning and management remains a major challenge and has not been applied in practice. During the past ten years, artificial land, including built-up land and traffic land, increased by 100% and 124%, respectively. Rapid expansion of artificial land had led to the decreases in carbon storage, biodiversity protection and recreation opportunity in Zengcheng. By 2013, artificial land accounted for 15% of the total land area but only provided 8%, 7%, 3% and 3% of carbon storage, water purification, biodiversity protection, and recreation opportunity services, respectively. Ecological land, which includes woodland, grassland and water bodies, can provide more of these services per unit area than artificial land. Thus, the artificial land expansion should be effectively controlled, especially in the southern region. In addition, the per capita area of built-up land in Zengcheng increased from 112 m 2 to 223 m 2 during the past ten years. Improving the efficiency of land use would be a better choice for future land use planning. New urbanization development zones should maintain a certain proportion of ecological land, which can facilitate the functionality, sustainability and vitality of urban ecosystems . To control the proportion of each land use type, ecological land protection policies in China, such as setting the "ecological control line" (Hong et al., 2017) should also be implemented in Zengcheng.
City planners and policy makers in Zengcheng should take corresponding ecological engineering and land management measures to improve ESs in future urban planning.
(1) CES value was higher in the northern region than that in other regions. However, the northern region suffered more severe soil erosion due to its steep terrain and heavy rainfall. This study proved that the SC scenario can significantly improve the soil conservation service, especially in the northern region. Thus, the "Grain for Green" policy should also be carried out in Zengcheng, especially in the northern region. Afforestation should be implemented on steep slopes and badly drained soil, whereas deforestation should be implemented in relatively flat, well-drained areas with loam soil (Rompaey et al., 2002). In addition, farmers and land managers should use mulch with vegetative residues to reduce soil and water loss (Prosdocimi et al., 2016). (2) In the central regions, the CES value decreased because of the rapid expansion of artificial land. To improve the services for climate regulation, biodiversity and recreation opportunity in this region, promoting built-up land utilizing more efficiently and intensively will be a good strategy in the future. In addition, this study showed that WB scenario can effectively improve the CES value, especially in the recreation service. As the cultural and recreation center of Zengcheng (OWZM, 2015), the central region should continue implementing the "Ecological corridor construction" project by developing greenways in the riparian and highway zones. This measure can reconcile political objectives and urban development while not only improving landscape quality but also providing new opportunities for public recreation and education (Ribeiro and Barão, 2006). (3) Southern region performed the worst in terms of climate regulation, water purification, biodiversity protection, and recreation opportunity. Industries of Zengcheng are mainly distributed in this region. By 2013, artificial land accounted for the largest proportion (40%) in this area. Therefore, this region had lower carbon storage and higher carbon emission than other regions. To balance carbon storage and emission, planners should not only control the spreading of built-up land and maintain a certain proportion of ecological land but also develop low carbon industrial methods, such as using renewable energy, enhancing energy efficiency, reducing the intensity of carbon emissions . Regarding water purification, the vegetation restoration in riparian zones in this study proved to be an effective way to improve water quality because these zones are important ecological barriers for water bodies. Additionally, developing green industry with low pollution and low consumption, adopting new water pollution mitigation technologies (Hasanbeigi and Price, 2015), and using pesticides that produce less pollution are all useful strategies for reducing water pollution. Biodiversity protection can be improved by enriching plant species, reducing landscape fragmentation and enhancing the connectivity of ecological corridors (Maes et al., 2012). Many provinces in China have implemented an "Ecological function zoning" policy (Liu et al., 2017) by establishing nature reserves to protect biodiversity. This policy can also be applied to future urban planning in Zengcheng. In terms of recreation opportunity, people prefer to travel to more natural habitats that surround areas of water (Paracchini et al., 2014). According to the population density and geographical characteristics, parks and green spaces should be reasonably and evenly distributed in Zengcheng to guarantee that these spaces are within walking distance. In addition, the connectivity, liquidity and quality of water bodies also should be improved to increase recreation opportunity.
Different scenarios in this study did not produce an obvious change in water supply. The water supply project in Zengcheng processed and provided 1.12 × 10 9 m 3 of water in 2013, and 98.8% of this water came from surface water. The total water consumption was 2.54 × 10 8 m 3 in 2013 (ZMSB, 2013). Therefore, the water resources in Zengcheng are relatively abundant, as we demonstrated that the surface water yield was 2.97 × 10 9 m 3 in 2013. For food production, the total consumption of grains, vegetables and fruits in 2013 were 1.19 × 10 5 , 8.37 × 10 4 , and 4.10 × 10 4 tons, respectively (http://tongji.cnki.net/ kns55/). Thus, food production under base case, agricultural and combined development scenarios can meet the local demand; however, grain production under woodland buffer and soil conservation scenarios are not sufficient to support local population. Although the woodland buffer and soil conservation scenarios perform the best in CES value, they improve most of ESs at the expense of reducing cropland and food production. To maintain a higher CES value and simultaneously meet the local food demands, farmers should increase the proportion of planted grain in Zengcheng. In addition, the agricultural sustainable intensification strategy (B. Zhang et al., 2017, X. Zhang et al., 2017 should be adopted to increase food yield without decreasing other ESs or using additional land. Such strategies include improving crop variety, optimizing fertilization, incorporating straw to cropland, developing advanced agricultural machinery, and developing effective irrigation technology.

Limitations and uncertainties
Some uncertainties and limitations exist in this study. First, the limitations in InVEST model and ROS model could lead to the inaccuracies of ES assessment. The climate regulation module assumes that none of land use types gain or lose carbon over time. However, many regions are actually recovering from past land use or are undergoing natural succession. This problem can be addressed by dividing the types into different age classes. Water supply module is based on annual averages, which neglects the temporal dimensions of water supply and hydropower production (Leh et al., 2013). Compared with other parameters in the module, water yield is highly sensitive to the changes of precipitation (Redhead et al., 2016). Future studies should consider the seasonal variability of precipitation in particular. Water purification module does not consider any chemical or biological interactions that may degrade the pollutants except filtration by terrestrial vegetation. Thus, the nutrient retention would be underestimated in the delivery process (Sharp et al., 2016). For this module, we can further assess the accuracy of ES simulations by field survey validation. The limitation in soil conservation module is that it only calculates the rill/inter-rill erosion processes in the USLE. The equation does not consider other sediment sources, including gully erosion, stream bank erosion, and mass erosion. R (climate) and K (soil) factors in the USLE have proven to be the most influential parameters for sediment export (Hamel et al., 2015). Local empirical data will be helpful for decreasing the uncertainty in these parameters. There also exists limitation in habitat quality module. This module assumes all threats are additive. However, in most of cases, the collective impact of multiple threats is much greater than the sum of individuals (Tallis et al., 2011). In addition, the ROS model mainly focuses on the recreational function of natural infrastructure, especially the water component. However, it ignores the recreational value of artificial infrastructure in built-up area. Therefore, this model underestimated the recreation service in Zengcheng. Second, four alternative scenarios in this study cannot represent all possible land use planning measures. Future studies should design various detailed scenarios to improve ESs according to the actual local situations and planning policies.

Conclusions
In this study, we developed a general structure to assess the spatiotemporal changes of multiple ESs and CES. Alternative scenarios were also designed to explore the optimal land use strategies that can provide higher CES values and minimize the trade-offs among various ESs. During the past ten years, rapid expansion of artificial land in Zengcheng had resulted in the decrease of CES values. We have suggested some ecological engineering and land management measures to improve the ESs in different regions according to the scenario analysis and local situations. City planners and policy makers should focus on improving the soil conservation in the northern region; recreation opportunity in the central region; and climate regulation, water purification, biodiversity protection, and recreation opportunity in the southern region. Scenario analysis showed that WB and SC scenarios perform best in CES and improve most of ESs at the expense of food production. Thus, the "Ecological corridor construction" project and the "Grain for Green" policy can be implemented on the premise of increasing food production. To minimize the trade-offs, sustainable intensification of agriculture should be adopted to increase the food supply without threatening other ESs or occupying additional land.