Spatial Heterogeneity of Typical Ecosystem Services and Their Relationships in Different Ecological– Functional Zones in Beijing–Tianjin–Hebei Region, China

: Recognizing changes in ecosystem services (ES) and their relationships is the basis of achieving sustainable regional development. Regional collaborative development has become the core strategy of the development of the Beijing–Tianjin–Hebei (BTH) region. However, sub regions have different ecological changes and relationships. Here, we quantify and map ES, including water yield, sediment retention, carbon sequestration and grain productive capacity in 2000, 2005, 2010 and 2015, using several biophysical models and explore the relationships of spatial correction, trade-offs and synergies among multiple ES in different spatial scales. Results across the four years show that the quality and variation tendency of ES from each region are spatially heterogeneous. The relationship between ES that are not signiﬁcant in the entire region shows different correlations in individual ecological–functional zones. From the perspective of regional disparity, the effect of land use factor and correlative mechanisms among ES are analyzed. To observe the spatiotemporal variations and relationships of ES in individual regions, land use management policies are proposed on the basis of the results of the relationships among ES.


Introduction
Ecosystems provide human beings with a wide range of material goods and services to survive and develop including a great series of material services (such as agricultural products, timber, water and other raw materials) and immaterial services (such as water and air purification, climate regulation and aesthetic benefits) [1,2]. All these material and immaterial benefits people obtain from ecosystems are generally called ecosystem services (ES) [3]. ES are generally classified into provision, regulation, support and culture services [2]. In the past 50 years, ecosystems have been destroyed in rapid urbanization. A total of 60% of ES have degraded worldwide, especially in areas with considerable human influence and dramatic land use changes [2]. As the second largest economy in the world, China has demonstrated remarkable achievements in economic growth since 1979. However, this rapid growth has come at a price and the environment is quickly deteriorating because of pollution and irrational utilization [4,5]. The Organization for Economic Cooperation and Development (OECD) pointed out that the economy of China is close to that of the developed world but its environment is similar to that of the poorest countries [6].
The Beijing-Tianjin-Hebei (BTH) region, which is the largest plain in Northern China, has also experienced rapid population growth and economic development. However, human activities, such as deforestation, mountain reclamation and excessive resource exploitation, have exacerbated the environmental problems in the region [7,8], thereby causing the decay of ES. The BTH region has a long history of agricultural reclamation. Considerable high-quality farmland has been lost due to the urban sprawl in the region and food production is projected to decrease from 109.61 × 10 6 t in 2013 to 106.14 × 10 6 -108.14 × 10 6 t in 2040 [9]. Furthermore, the entire BTH region is classified to be under absolute scarcity level by UN standards (2013) [10]. Groundwater level has declined continuously at a rate of 0.36 m/a in the past 50 years [11]. Meanwhile, the degradation and destruction of forests and grasses has led to massive soil erosion and desertification [12]. In the past few decades, the mountains Taihang and Yan-shan have been subject to serious soil erosion because of vegetation deterioration caused by overgrazing and overfarming [13]. Photosynthesis of vegetation can improve air quality and regulate climate. However, the urbanization and industrialization in the BTH region does not only reduce urban green space, thereby lessening the capacity of carbon sequestration but also produces considerable carbon emissions, thus leading to the heat island effect and air pollution [14,15].
With respect to the dilemma between economic development and environmental deterioration, many studies have concentrated on revealing the relationship among various ES and seeking an improved solution to achieve a mutually beneficial objective. Nelson et al. [3] modeled the changes in ecosystem services, biodiversity conservation, commodity production in a spatially explicit manner, found tradeoffs between them in the Willamette Basin, Oregon and made a deal with more effective, efficient and defensible natural resource decisions. Goldstein et al. [16] evaluated the environmental and financial implications of seven planning scenarios in Kamehameha Schools and found that tradeoffs existed between carbon storage and water quality as well as between environmental improvement and financial return, which helps to guide the local land-use decisions by involving tradeoffs between private and public interests. Relevant researches were also conducted in the BTH region. By developing the relationship between net primary productivity (NPP) and grain yield, Peng et al. [17] quantified and mapped the grain production of farmland and found that the main high-production areas were mainly located in the central plain of the BTH region and the eastern piedmont plain. The plain had low soil retention because of the low potential soil loss in farmland. Li et al. [11] evaluated regional water security in the BTH region via a freshwater ES flow model and concluded that freshwater ecosystem service provision was high in mountainous areas, especially in the Yan-shan and Taihang mountains and low in the Northern China Plain and urban areas. Nian et al. [18] analyzed the net capacity of ecosystems that provide carbon fixation and oxygen release (CFOR) in the BTH region using the "provision-reception" framework and the net CFOR exhibited a declining trend from the northwest to the southeast. To fulfill people's demand for ES in the BTH region, Zhang et al. [19] proposed an evaluation framework of landscape ecological security pattern integrating ES assessment and landscape connectivity analysis with human ecological demand to identify the local ecological characteristics. Although these studies have provided good approaches and frameworks to quantify ES in the BTH region, no unified suggestions or environmental policies can be applied in the whole region because of the significant spatial heterogeneity within BTH.
In 2015, the outline of collaborative development of Beijing, Tianjin and Hebei Province was released by the Chinese government to strengthen the ecosystem management in the region. Formulating effective policies that target regional characteristics is a major challenge in improving environmental performance, especially at the sub region level of the BTH area. Therefore, a clear exploration is needed to deepen the correlations among ES and the spatiotemporal changes over time in different physical-geographical regions [7,20].
The scientific interests of this study focused on the following: (1) determining the correlation among different ES in the BTH region and their change over time, (2) identifying variance in these correlations in the different ecological-functional zones of the BTH region and (3) learning the suggestion of the sub regional variance in these correlations to local management. Four types of ES, namely, water yield (WY), soil retention (SR), carbon sequestration (CS) and grain production (GP), in 2000, 2005, 2010 and 2015, were quantified and mapped. The correlation among the different ES in various ecological-functional zones was analyzed by Spearman nonparametric tests to reveal the spatial relationship. Finally, on the basis of the results of the correlations among ES, we proposed sustainable development advice for the different zones.

Background of the Study
Changes in one ES can cause changes in others through three main relationships as follows.
(1) The synergy describes a positive reaction of ES to a change in the driver [21]; (2) trade-off is an antagonistic situation that involves losing one quality of an ES to gain another [22]; and (3) compatibility means that a change in ES does not significantly influence others. ES relationships have been discussed previously. Several scholars found a synergy between supply service and biodiversity [3,23] and others found a trade-off correlation between provision and regulation services [24].
At present, common research methods of analyzing trade-off or synergy correlation mainly include map comparison, scenario analysis and trade-offs using optimized landscapes [25]. The first method maps each ES and analyzes their correlation cell by cell via spatial overlay [26]. As the most commonly used method, scenario analysis sets a series of situations and analyzes trade-off or synergy correlation among ES under various conditions [27]. The third method quantifies ES with biophysical or statistical models and then compares trade-off or synergy correlation on the basis of optimized land use patterns [28].
The trade-off or synergy correlation among ES is not uniform because of the varying geographical characteristics. This correlation depends on the ecosystem compositions, structures, conditions [3,7], complex interactions between physical (e.g., topography, geology) and biological (e.g., vegetation and micro-organism) factors and land use and management [8,29]. Jiang et al. [4] and Qin et al. [30] analyzed the trade-off and synergy between ES in the Three-river Headwater area and the Guanzhong-Tianshui economic region, respectively. The studies found a trade-off between water yield provision and carbon sequestration in the former region, which was in the hinterland of the Qinghai-Tibet Plateau and had a plateau continental climate and a strong synergy in the latter region, which was in the monsoon area of transition from semi-humid to semi-arid. REDD+ policies have attempted to achieve mutually beneficial solutions to enhance the carbon sequestration service and the conservation of tropical biodiversity; however, the results demonstrated an apparent trade-off between them [31]. Therefore, economic and ecological criteria should be regionally balanced in the application of land use policies because various regions respond differently to ES.

Study Area
The BTH region (36 • 01 -42 • 37 N, 113 • 04 -119 • 53 E) is in Northern China and includes Beijing City, Tianjin City and Hebei Province. The region has the typical temperate semi-humid to semi-arid monsoon climate. Rainfall mainly occurs in summer and presents a unimodal distribution. Average annual rainfall is roughly 538 mm and presents a decreasing spatial pattern from eastern coastal to western inland. The northwest region has an undulating and mountainous topography with steep slopes and the southeast plain is characterized as an agricultural production base. The region has an elevation range of 0-2800 m and includes highland, mountain land, piedmont plain, low plain and littoral plain. On the basis of the physiographic conditions, the BTH region can be divided into four ecological-functional zones (Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, 2000), namely, deciduous broad-leaved forest ecoregion of

Materials
Land use datasets of the BTH region were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn), including farmland, forestland, grassland, water and wetland, settlement, desert and others. The boundaries of the ecological-functional zone were provided by the Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences (http://www.ecosystem.csdb.cn/). Raster-formatted soil data, that is, soil texture, soil types, organic matter content, reference soil depth, root-restricting layer depth and plant available water content (PAWC), were acquired from the Harmonized World Soil Database. Global digital elevation model data were downloaded from the United States National Aeronautics and Space Administration. Climate data were obtained from China National Meteorological Information Center. The boundaries of the BTH region were obtained from China Data Sharing Infrastructure of Earth System Science (2013). Vector data of watershed were required in water yield assessment; the entire study area was delineated into 122 watersheds using the Integrated Valuation of Ecosystem Services and Tradeoff model (InVEST model) [32]. All the data were uniformly resampled at a resolution of 1000 m. Average annual reference evapotranspiration was provided by the MOD16A3 product from the Numerical Terradynamic Simulation Group (http://www.ntsg.umt.edu/) but 2015 data was unavailable and thus replaced by 2014 data. Annual NPP was provided by the MYD17A3H Version 6 product from the Land Processes Distribution Active Archive Center (https://lpdaac.usgs.gov/).

Methodology
In this study, water yield and soil retention were evaluated by InVEST model, which was a model designed to map and quantify the distribution of ES through a set of standalone but linkable models and a tool for exploring how changes in ecosystems [32]. According to photosynthesis equation, NPP can be convert to the amount of carbon sequestration, therefore, NPP was used to represent the carbon sequestration, the capability of vegetation to absorb carbon dioxide and was calculated by the Carnegie-Ames-Stanford approach [33]. Grain production service was calculated by the gradation on agricultural land quality (GALQ) model, which is derived from the Agricultural

Materials
Land use datasets of the BTH region were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn), including farmland, forestland, grassland, water and wetland, settlement, desert and others. The boundaries of the ecological-functional zone were provided by the Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences (http://www.ecosystem.csdb.cn/). Raster-formatted soil data, that is, soil texture, soil types, organic matter content, reference soil depth, root-restricting layer depth and plant available water content (PAWC), were acquired from the Harmonized World Soil Database. Global digital elevation model data were downloaded from the United States National Aeronautics and Space Administration. Climate data were obtained from China National Meteorological Information Center. The boundaries of the BTH region were obtained from China Data Sharing Infrastructure of Earth System Science (2013). Vector data of watershed were required in water yield assessment; the entire study area was delineated into 122 watersheds using the Integrated Valuation of Ecosystem Services and Tradeoff model (InVEST model) [32]. All the data were uniformly resampled at a resolution of 1000 m. Average annual reference evapotranspiration was provided by the MOD16A3 product from the Numerical Terradynamic Simulation Group (http://www.ntsg.umt.edu/) but 2015 data was unavailable and thus replaced by 2014 data. Annual NPP was provided by the MYD17A3H Version 6 product from the Land Processes Distribution Active Archive Center (https://lpdaac.usgs.gov/).

Methodology
In this study, water yield and soil retention were evaluated by InVEST model, which was a model designed to map and quantify the distribution of ES through a set of standalone but linkable models and a tool for exploring how changes in ecosystems [32]. According to photosynthesis equation, NPP can be convert to the amount of carbon sequestration, therefore, NPP was used to represent the carbon sequestration, the capability of vegetation to absorb carbon dioxide and was calculated Sustainability 2018, 10, 6 5 of 23 by the Carnegie-Ames-Stanford approach [33]. Grain production service was calculated by the gradation on agricultural land quality (GALQ) model, which is derived from the Agricultural Land Classification and Gradation Project by the Ministry of Land and Resources of China. Soil retention, carbon sequestration and grain productive capacity are presented at pixel scale, whereas water yield is shown at sub-watershed scale. To facilitate comparison, we indicated the values of soil retention, carbon sequestration and grain productive capacity at sub-watershed scales. Spearman nonparametric tests were performed to check for the correlations among the ES above to illustrate the regional heterogeneity in the four ecological-functional zones. The main workflow fluxes are presented in Figure 2. Land Classification and Gradation Project by the Ministry of Land and Resources of China. Soil retention, carbon sequestration and grain productive capacity are presented at pixel scale, whereas water yield is shown at sub-watershed scale. To facilitate comparison, we indicated the values of soil retention, carbon sequestration and grain productive capacity at sub-watershed scales. Spearman nonparametric tests were performed to check for the correlations among the ES above to illustrate the regional heterogeneity in the four ecological-functional zones. The main workflow fluxes are presented in Figure 2.

Water Yield Assessment
In the InVEST model, water yield was calculated on the basis of the Budyko curve [34] and annual average precipitation as follows: where Ym denotes the water yield in watershed m; Yxj and AETxj are the water yield and annual actual evapotranspiration on pixel x with land use type j, respectively; and Px denotes the annual precipitation on pixel x. In Equation (2), AETxj/Px approximates the Budyko curve estimated by Zhang et al. [35] as follows: where Rxj denotes the dimensionless Budyko dryness, which is the ratio of potential evapotranspiration to precipitation and ωx is the plant available water coefficient on pixel x. The input data included annual precipitation, annual potential evapotranspiration (ET0), soil depth, PAWC, land use and land cover and root depth. ET0 was calculated by the FAO Penman-Monteith equation [36]. The root depth was referenced from Canadell et al. [37]. This empirical model used for Budyko theory has been tested at larger than pixel scales [38]; thus, the output summed or averaged at the sub-watershed scale has a higher credibility and is more interpretable than that at the pixel scale.

Water Yield Assessment
In the InVEST model, water yield was calculated on the basis of the Budyko curve [34] and annual average precipitation as follows: where Y m denotes the water yield in watershed m; Y xj and AET xj are the water yield and annual actual evapotranspiration on pixel x with land use type j, respectively; and P x denotes the annual precipitation on pixel x. In Equation (2), AET xj /P x approximates the Budyko curve estimated by Zhang et al. [35] as follows: where R xj denotes the dimensionless Budyko dryness, which is the ratio of potential evapotranspiration to precipitation and ω x is the plant available water coefficient on pixel x. The input data included annual precipitation, annual potential evapotranspiration (ET 0 ), soil depth, PAWC, land use and land cover and root depth. ET 0 was calculated by the FAO Penman-Monteith equation [36]. The root depth was referenced from Canadell et al. [37]. This empirical model used for Budyko theory has been tested at larger than pixel scales [38]; thus, the output summed or averaged at the sub-watershed scale has a higher credibility and is more interpretable than that at the pixel scale.

Soil Retention Assessment
Soil retention was calculated by the universal soil loss equation (USLE) [39], which computes the average annual amount of potential soil loss in a particular area on the basis of its land cover Sustainability 2018, 10, 6 6 of 23 characteristics and topography. Soil retention was calculated as soil loss without vegetation cover and soil erosion control minus soil loss for the current land use and soil erosion management. Soil retention can be expressed mathematically as follows: where ∆A is the soil retention; A 0 and A v represent the potential soil erosion without vegetation cover and the soil erosion under the current land cover and management conditions, respectively; R, K, L and S are the rainfall erosivity, soil erodibility, slope length and slope angle factors, respectively; and C and P are the current vegetation cover and erosion control factors, respectively. R factor is calculated by a modified formula based on the data of annual and monthly precipitation [39,40]. Soil texture is a principal factor that affects K; the structure, organic matter and permeability contribute as well. We calculated factor K by the erosion/productivity impact calculator model [41]. The  [42][43][44] according to the actual situation in the BTH region.

Grain Productive Capacity Assessment
On the basis of the relationship between soil and crop production, the GALQ model quantifies the grain productive capacity of agricultural land in China. Different from the real productivities which are not only affected by the quality factors of farmland but also restricted by the regional agricultural input (including fertilizer utilization, pest control and agricultural machinery) and agricultural technology management (including horticultural techniques and labor input), the grain productive capacity focuses on natural endowments of cropping system rather than agricultural mechanization levels. In this model, grain productive capacity is affected by natural endowments, including local climatic conditions (such as light, temperature and water) and soil quality [45]. To emphasize the natural attributes of agricultural land, the grain productive capacity calculated by the GALQ model refers to the highest yield of fine crop variety per unit area that might be obtained assuming all or part of the production factors are in the optimum state and take no account of agricultural input and agricultural technology management; therefore, the value is higher than the real grain yield. The cardinal principle of this model is building a regression equation between the physical quality of the agricultural land and the recorded highest grain yield in local history to estimate maximum grain productivity. Using a database of the experimental samples, this work assessed the grain productive capacity in the BTH region (Table 1). Here, y denotes the grain productive capacity (per unit yield) and x denotes the physical quality of the agricultural land, including cultivated land, yard and facility farmland.

Water Yield
The spatial heterogeneity of the water yield is visualized in Figure 3. High values of water yield were distributed in the eastern region, including Zone II, northeast of Zone III and east of Zone I. The low-value areas were in the northwest, including Zone IV, southwest of Zone III and northwest of Zone I. In general, local climate and land use affected the water yield. In comparison with the region with small forestland and settlement land, the western area of Zone I and the southeast of Zone III, which had forestland and settlement land, evaporated more water. Zone IV, which was part of the Inner Mongolia Plateau region, had a relatively high altitude and the low precipitation led to a low water yield.
The annual average water yield had a remarkable increase from 193.41 mm/km 2 in 2000 to 240.35 mm/km 2 in 2015 ( Figure 4). The water yield in Zone II was higher than that in other regions in all years except 2005; that in Zone IV was lower than other regions in all years except 2010. low-value areas were in the northwest, including Zone IV, southwest of Zone III and northwest of Zone I. In general, local climate and land use affected the water yield. In comparison with the region with small forestland and settlement land, the western area of Zone I and the southeast of Zone III, which had forestland and settlement land, evaporated more water. Zone IV, which was part of the Inner Mongolia Plateau region, had a relatively high altitude and the low precipitation led to a low water yield. The annual average water yield had a remarkable increase from 193.41 mm/km 2 in 2000 to 240.35 mm/km 2 in 2015 ( Figure 4). The water yield in Zone II was higher than that in other regions in all years except 2005; that in Zone IV was lower than other regions in all years except 2010.   low-value areas were in the northwest, including Zone IV, southwest of Zone III and northwest of Zone I. In general, local climate and land use affected the water yield. In comparison with the region with small forestland and settlement land, the western area of Zone I and the southeast of Zone III, which had forestland and settlement land, evaporated more water. Zone IV, which was part of the Inner Mongolia Plateau region, had a relatively high altitude and the low precipitation led to a low water yield. The annual average water yield had a remarkable increase from 193.41 mm/km 2 in 2000 to 240.35 mm/km 2 in 2015 ( Figure 4). The water yield in Zone II was higher than that in other regions in all years except 2005; that in Zone IV was lower than other regions in all years except 2010.

Sediment Retention
Sediment retention had high spatial distribution patterns in Zones I and IV, where the topography was mainly mountain and land coverage was mainly forest and grass. The low value was in Zones II and III, where plain was the main landform ( Figure 5). The sediment erosion of the entire region increased from 990.07 ton/hm 2 in 2000 to 1119.53 ton/hm 2 in 2015 with a slight decline in 2015 ( Figure 6). Zoning in the ecological-functional Zones I and IV had a similar increasing tendency as the entire region and peaked in 2010 (1714.97 ton/hm 2 in Zone I and 1052.17 ton/hm 2 in Zone IV). The soil retention in Zone II decreased and then increased and it reached a historical low in 2005. The only zone whose soil retention steadily increase from beginning to end was Zone III (from 534.19 ton/hm 2 in 2000 to 698.93 ton/hm 2 in 2015).

Sediment Retention
Sediment retention had high spatial distribution patterns in Zones I and IV, where the topography was mainly mountain and land coverage was mainly forest and grass. The low value was in Zones II and III, where plain was the main landform ( Figure 5). The sediment erosion of the entire region increased from 990.07 ton/hm 2 in 2000 to 1119.53 ton/hm 2 in 2015 with a slight decline in 2015 ( Figure 6). Zoning in the ecological-functional Zones I and IV had a similar increasing tendency as the entire region and peaked in 2010 (1714.97 ton/hm 2 in Zone I and 1052.17 ton/hm 2 in Zone IV). The soil retention in Zone II decreased and then increased and it reached a historical low in 2005. The only zone whose soil retention steadily increase from beginning to end was Zone III (from 534.19 ton/hm 2 in 2000 to 698.93 ton/hm 2 in 2015).

Sediment Retention
Sediment retention had high spatial distribution patterns in Zones I and IV, where the topography was mainly mountain and land coverage was mainly forest and grass. The low value was in Zones II and III, where plain was the main landform ( Figure 5). The sediment erosion of the entire region increased from 990.07 ton/hm 2 in 2000 to 1119.53 ton/hm 2 in 2015 with a slight decline in 2015 ( Figure 6). Zoning in the ecological-functional Zones I and IV had a similar increasing tendency as the entire region and peaked in 2010 (1714.97 ton/hm 2 in Zone I and 1052.17 ton/hm 2 in Zone IV). The soil retention in Zone II decreased and then increased and it reached a historical low in 2005. The only zone whose soil retention steadily increase from beginning to end was Zone III (from 534.19 ton/hm 2 in 2000 to 698.93 ton/hm 2 in 2015).

Carbon Sequestration
The carbon sequestration represented by NPP showed a similar spatial pattern as the sediment retention from 2000 to 2015. Zones I and IV, which had more forest and grass than the other zones, had high carbon sequestration values; meanwhile, low values were found in the other regions, which had less forest and grass than Zones I and IV (Figure 7). The highest annual average NPP in the entire BTH region was 324.54 g·C/m 2 in 2000 and then decreased to 278.63 g·C/m 2 in 2015 gradually (Figure 8). At the regional scale, the NPP in Zones I, II and III decreased by 11.64%, 18.67% and 30.88% since 2000, respectively but the tendency in Zone IV decreased initially and then increased.

Carbon Sequestration
The carbon sequestration represented by NPP showed a similar spatial pattern as the sediment retention from 2000 to 2015. Zones I and IV, which had more forest and grass than the other zones, had high carbon sequestration values; meanwhile, low values were found in the other regions, which had less forest and grass than Zones I and IV (Figure 7). The highest annual average NPP in the entire BTH region was 324.54 g·C/m 2 in 2000 and then decreased to 278.63 g·C/m 2 in 2015 gradually (Figure 8). At the regional scale, the NPP in Zones I, II and III decreased by 11.64%, 18.67% and 30.88% since 2000, respectively but the tendency in Zone IV decreased initially and then increased.

Carbon Sequestration
The carbon sequestration represented by NPP showed a similar spatial pattern as the sediment retention from 2000 to 2015. Zones I and IV, which had more forest and grass than the other zones, had high carbon sequestration values; meanwhile, low values were found in the other regions, which had less forest and grass than Zones I and IV (Figure 7). The highest annual average NPP in the entire BTH region was 324.54 g·C/m 2 in 2000 and then decreased to 278.63 g·C/m 2 in 2015 gradually (Figure 8). At the regional scale, the NPP in Zones I, II and III decreased by 11.64%, 18.67% and 30.88% since 2000, respectively but the tendency in Zone IV decreased initially and then increased.

Grain Productive Capacity
Significant differences in the spatial distribution of grain productive capacity are shown in Figure 9. Farmland was mainly distributed in Zones II and III; these farmland regions were of high quality because of the integration of fertile soil, sufficient sunlight and abundant rainfall. By contrast, a small amount of farmland was scattered in Zone I and the grain productive capacity was low. Most of the farmland distributed in Zone IV had medium-low yield.
The average grain productive capacity of the entire region had a slight decrease from 5.09 × 10 5 kg/km 2 in 2000 to 4.98 × 10 5 kg/km 2 in 2015 (Figure 10), which indicated that the quality of farmland in the whole region had degraded. Meanwhile, the grain productive capacity of Zones I, II and III all experienced a decreasing trend, especially that of Zone III, which dropped by 6.58%. By contrast, the production capability of Zone IV had a slight increase of roughly 0.67%.

Grain Productive Capacity
Significant differences in the spatial distribution of grain productive capacity are shown in Figure 9. Farmland was mainly distributed in Zones II and III; these farmland regions were of high quality because of the integration of fertile soil, sufficient sunlight and abundant rainfall. By contrast, a small amount of farmland was scattered in Zone I and the grain productive capacity was low. Most of the farmland distributed in Zone IV had medium-low yield.
The average grain productive capacity of the entire region had a slight decrease from 5.09 × 10 5 kg/km 2 in 2000 to 4.98 × 10 5 kg/km 2 in 2015 (Figure 10), which indicated that the quality of farmland in the whole region had degraded. Meanwhile, the grain productive capacity of Zones I, II and III all experienced a decreasing trend, especially that of Zone III, which dropped by 6.58%. By contrast, the production capability of Zone IV had a slight increase of roughly 0.67%.

Grain Productive Capacity
Significant differences in the spatial distribution of grain productive capacity are shown in Figure 9. Farmland was mainly distributed in Zones II and III; these farmland regions were of high quality because of the integration of fertile soil, sufficient sunlight and abundant rainfall. By contrast, a small amount of farmland was scattered in Zone I and the grain productive capacity was low. Most of the farmland distributed in Zone IV had medium-low yield.
The average grain productive capacity of the entire region had a slight decrease from 5.09 × 10 5 kg/km 2 in 2000 to 4.98 × 10 5 kg/km 2 in 2015 (Figure 10), which indicated that the quality of farmland in the whole region had degraded. Meanwhile, the grain productive capacity of Zones I, II and III all experienced a decreasing trend, especially that of Zone III, which dropped by 6.58%. By contrast, the production capability of Zone IV had a slight increase of roughly 0.67%.

Correlation among ES
Through Spearman correlation analysis, a positive correlation was found between the status quo of water yield and soil retention/grain production capability/carbon sequestration, soil retention and carbon sequestration in each year (Table 2) and the consistent correlation for the four point-in-times separated by five years can guarantee a stable relationship between each two ES; changes in water yield and soil retention in 2000-2015 (r 2 = 0.469 **) were positively correlated as well. The changes in water yield in 2000-2015 were negatively correlated with those of grain production capability and carbon sequestration (r 2 = −0.310 ** and −0.514 **, respectively) but an adverse relationship existed between soil retention and carbon sequestration (r 2 = 0.426 **). It means that there is no uniformity between the correlation of status quo and that of changes for the same pair of ES.
In the entire study area, neither status quo nor the changes in grain production capability and soil retention/carbon sequestration were apparent in 2000-2015. However, this kind of nonobvious and unstable correlation changes at sub regional scale. Grain production capability had a positive correction with carbon sequestration in Zones II and III and a negative correction in Zones I and IV every four years (Table 3). In addition, the changes in grain production capability and carbon sequestration in 2000-2015 indicated a negative correction in Zones II and III and a positive one in Zones I and IV. A negative correlation was found between grain production capability and soil retention in Zones I and IV; no correlation existed in Zones II and III. Meanwhile, the changes in those within 2000-2015 were uncorrelated in Zones I, II and III and had a poor negative correction in Zone IV. The correlation of ES behaved differently according to the various terrains and vegetation. For instance, the grain production capability and soil retention were correlated negatively in the mountainous regions, such as Zone I and IV but not correlation in plain areas such as Zones II and III. As for the production capability and carbon sequestration, it also showed different relationships between them in mountainous and plain areas.

Impact of Land Use and Land Cover Change on ES
Land use and land cover change is the main factor that affects ecosystems. Land use can significantly affect ES [46]. Figures 11-13 reflect the trends of these four ES across different land types during the study period, which were as follows. (1) In farmland, water yield dropped from 2000 to 2005, reached the minimum (195.66 mm/km 2 ) and then increased gradually. Soil retention slightly fluctuated from 422.28 ton/hm 2 in 2000 to 496.02 ton/hm 2 in 2010, whereas carbon sequestration decreased each year. (2) In forestland, water yield and soil retention peaked in 2010 (204.88 mm/km 2 and 1549.17 ton/hm 2 , respectively) and then dropped slightly. Carbon sequestration decreased initially and tended to stabilize. (3) The trends of water yield and carbon sequestration for grassland were the same as those of forestland, whereas soil retention fluctuated up or down at roughly 1200 ton/hm 2 .
Coefficient of variation is used to reflect the trend of water yield in Figure 11. The trend decreased and then tended to stabilize within the range of 0.2-0.6. For the soil retention service (Figure 12), the coefficient of variation of farmland was considerably higher than that of forestland and grassland; therefore, in comparison with the soil retention of forestland and grassland, that of farmland significantly varied by region, being susceptible to natural regional factors. The coefficient of variation for carbon sequestration (Figure 13) remained within 0.2-0.35 and remained stable and this pattern illustrated the carbon sequestration of different land uses was not easily influenced by environmental deviation.    significantly affect ES [46]. Figures 11-13 reflect the trends of these four ES across different land types during the study period, which were as follows. (1) In farmland, water yield dropped from 2000 to 2005, reached the minimum (195.66 mm/km 2 ) and then increased gradually. Soil retention slightly fluctuated from 422.28 ton/hm 2 in 2000 to 496.02 ton/hm 2 in 2010, whereas carbon sequestration decreased each year. (2) In forestland, water yield and soil retention peaked in 2010 (204.88 mm/km 2 and 1549.17 ton/hm 2 , respectively) and then dropped slightly. Carbon sequestration decreased initially and tended to stabilize. (3) The trends of water yield and carbon sequestration for grassland were the same as those of forestland, whereas soil retention fluctuated up or down at roughly 1200 ton/hm 2 . Coefficient of variation is used to reflect the trend of water yield in Figure 11. The trend decreased and then tended to stabilize within the range of 0.2-0.6. For the soil retention service (Figure 12), the coefficient of variation of farmland was considerably higher than that of forestland and grassland; therefore, in comparison with the soil retention of forestland and grassland, that of farmland significantly varied by region, being susceptible to natural regional factors. The coefficient of variation for carbon sequestration (Figure 13) remained within 0.2-0.35 and remained stable and this pattern illustrated the carbon sequestration of different land uses was not easily influenced by environmental deviation.   An index Ci was obtained to analyze the relationships among the different ES in the four ecological-functional zones and determine the intensity of land use change using land use transfer matrices as follows:  An index Ci was obtained to analyze the relationships among the different ES in the four ecological-functional zones and determine the intensity of land use change using land use transfer matrices as follows: An index C i was obtained to analyze the relationships among the different ES in the four ecological-functional zones and determine the intensity of land use change using land use transfer matrices as follows: where C i is the percentage of land use change type i among all change types, A i is the area of change type I and A is the total area of all land use change in 2000-2015. A high value of C i demonstrates that this type of land use change plays a leading role in regional land use changes.
The results (Table 4) indicated the following. (1) The land use changes from farmland, grassland and forestland to settlement (SEL) were the main change types in Zone I and were accompanied with decreases in soil retention, water yield and carbon sequestration; (2) the conversion from farmland to settlement dominated the land use changes in Zone II, thereby leading to a decline in grain production, water yield and carbon sequestration; (3) all four types of ES in Zone III decreased because of the conversion from farmland, water and wetland (WAW) to settlement; and (4) the dominant conversion in Zone IV was slightly different from that of the other zones, which was from grass to farmland, thereby resulting in an increase in grain production and water yield but a decrease in soil retention and carbon sequestration. In the BTH region, urban land expanded by 71% from 1990 to 2000; roughly 74% of the new urban land was converted from arable land [47]. This ratio increased to 81.72% from 2000 to 2015. The development of Zone III, which was the concentration district of the metropolis, mainly depended on secondary and tertiary industries. The urbanization progress directly led to a reduction in the average grain productive capacity in Zone III because of the encroachment of high-quality farmland (Table 3). Meanwhile, the development of Zone II relied mostly on primary industry; thus, this zone had more farmland than Zone III and had inconspicuous decreases in grain productive capacity.

Correlative Mechanisms among ES
ES are closely interlinked in a complex manner and characterized by a nonlinear relationship [48]. Bennett et al. [21] argued that ES are interrelated because of the (1) effects of common drivers or (2) interactions among ES per se. In the present research, water yield and sediment retention were correlated possibly because they shared common hydrological processes and drivers of precipitation and land use/cover change. For instance, the water retention capability of soil was correlated with water yield [32] because the soil loss might be caused by water flows. Specifically, the high water yield was consistently linked with high water retention capability in the soil. The lower the water flow, the lower the soil export. In addition, the growth of vegetation could improve the physical structure and then enhance the water retention capability of the soil. The high vegetation coverage limited soil moisture from achieving saturation and forming water flows [49,50] and consequently reduced soil erosion. This phenomenon explained the positive correlation between carbon sequestration and water yield/soil retention. In the current study, NPP was used to represent carbon sequestration; the higher the NPP, the larger the vegetation coverage [51]. Furthermore, these results also indicated that the areas with abundant vegetation except farmland consistently showed high values of water yield and soil retention (Tables 2 and 3). The positive correlations between the differences in soil retention and water yield/carbon sequestration from 2000 to 2012 indicated the synergies between them on the Loess Plateau [4]. This study also showed negative relationships between water yield and carbon sequestration. In other words, a trade-off existed between water yield and carbon sequestration and this finding agreed with the studies in the Loess Plateau [48,52]. The Loess Plateau and the BTH region are semi-arid and semi-humid regions. With an increase in vegetation coverage, more rainwater retained by leaves evaporate and more rainwater are absorbed by vegetation on the surface and then evaporate. Water yield also decreases.
The study conducted by Li et al. [53] in 2010 on the BTH region showed a weak negative correlation between gain production and carbon sequestration/soil retention. The correlation analysis conducted by the present work in the four ecological-functional zones yielded the following results. (1) Production capability and soil retention were negatively correlated in Zones I and IV but not correlation in Zones II and III. The grain production in the mountain area easily caused soil erosion; thus, the grain production in Zones I and IV, where the terrain was predominantly mountain, was not conducive to soil retention. However, the grain production caused slight soil erosion in Zones II and III because of the flat terrain.
(2) Grain production and soil retention had a weak trade-off in Zone IV and no correlation in Zones I, II and III. (3) In Zones I and IV, the forest and grass had a strong capability for carbon sequestration and a low capability for grain production. In Zones II and III, farmland had a high capability for carbon sequestration and grain production. Grain production and carbon sequestration had a negative correlation in Zones I and IV and a positive correlation in Zones II and III. (4) Grain production and carbon sequestration exhibited synergy in Zones I and IV but a trade-off in Zones II and III. The synergy in Zones I and IV agreed with those found by the study conducted on the Loess Plateau in northern part of middle China [48,52]. The trade-offs in Zones II and III were consistent with those shown in the study conducted on the Sanjiang Plain in northeastern China for the period of 1992-2012 [54].

Implications for Regional Sustainable Development
In 2014, the Chinese President Xi proposed a coordinated development of the BTH region. The generation of an eco-environmental protection strategy should eliminate the restrictions imposed by the administrative boundaries and take the ecological function in its entirety for the further regional sustainable development. Recently, the Chinese national and local governments implemented a series of measures on territory development plans and ecological environmental protection. Although certain evidence indicates that these measures have certain effects, considerable improvement is still needed for strengthening the implementation of environmental policies [55]. These policies have been applied in different ecological-functional zones and the results from other studies are listed in Table 5.
On the basis of the characteristics of land use conversion and the ecological environment problems in the different regions from 2000 to 2015, we combined the spatial correlation, trade-offs and synergies among ES and then proposed suggestions for sustainable regional development. To protect and improve the ecological environment, SLCP plans to stop the farming of sloping farmlands, which is likely to cause soil and water loss and plant trees to restore forest vegetation.
Increase forestland and grassland but reduce farmland The forestland was restored well in the high-terrain-gradient area and the predominant distribution area of unused land and grassland was gradually diminished [56]. The project also had strong influences on the soil quantity improvement and vegetation coverage increase [57].

Beijing-Tianjin Sandstorm Source Control Project (2000-2010) (BTSSCP)
BTSSCP is a sandstorm control policy implemented in Beijing, Tianjin, Hebei, Shanxi and Inner Mongolia to curb the expansion of desertification land and improve the ecological environment.

Increase forestland and grassland
The implementation of BTSSCP has significantly improved the growth conditions of grassland and forestland vegetation [58].
The Three-north Shelterbelt Program (1978-2050) (TNSP) TNSP zone covers the northwest, north and northeast regions of China and consists of three stages (1978-2000, 2001-2020 and 2021-2050); the key objective is to increase regional forest coverage in arid and semi-arid China from 5% to 15%.

Increase forestland
The Chinese government claims that TNSP increased forest cover from 5% in 1978 to more than 13% in 2017. TNSP improved the forest cover of the area effectively and positively affected carbon sequestration [59]. However, Cao et al. and Wang et al. [60,61] strongly argued that the large-scale afforestation failed to address the desertification in some arid and semi-arid regions.
The Policy of Dynamic Equilibrium of the total Cultivated Land (1997-) (PDDCL) To maintain a dynamic balance of cultivated land, provincial governments are required to reclaim new farmland, rehabilitate damaged land, or reuse deserted land to compensate the losses of farmland.

Protect farmland
The farmland expanded to high-terrain-gradient area to compensate the region occupied by the construction of infrastructure facilities in the plain of BTH [56]. (1) According to the demand of population growth for agricultural products, the government set a certain amount of farmland that cannot be occupied.
(2) HSPFCP plans to construct farmlands with concentration, fully supported facilities, stable and high yield, friendly environment and high anti-disaster capability and can adopt to the modern agricultural production and operation modes by the rural land reclamation and readjustment.
Protect farmland and improve the quality of farmland Although PBAL and HSPFCP have protected a certain amount of high-quality farmland and improved their quality to some extent, the fast urbanization in this area has led to significant farmland loss [62,63]. According to Zhang et al. [62], the main reason of the ecosystem service decrease in the BTH region is the increase in artificial land and loss of cropland.

National Wetland Protection (2004-2010) (NWP)
NWP aimed to strengthen the protection of wetland by improving their ecological environment.

Protect wetland
The wetlands remained threatened because of unsustainable usage. The enforcement of the wetland protection law was weak [64]. The municipal government afforested to make the capital an "ecological livable city" by improving the urban ecological environment in the Beijing plain areas.
Increase forestland and grassland but decrease farmland In Beijing, the afforestation efforts covered 1.05 million acres and planted more than 5400 million trees but the farmland decreased at the same time [65,66].
Ecological engineering projects during preparation for 2008 Olympics (2000Olympics ( -2008 To improve the ecological environment during the Beijing Olympic Games period, the capital conducted a series of urban green space construction efforts.
Increase forestland and grassland but decrease farmland The vegetation increased significantly by occupying the farmland [66,67].
(1) The Protection of Basic Agricultural Land (1994-) (PBAL); (2) High-standard Primary Farmland Construction Project (2011-) (HSPFCP) (1) According to the demand of population growth for agricultural products, the government set a certain amount of farmland that cannot be occupied.
(2) HSPFCP plans to construct farmlands with concentration, fully supported facilities, stable and high yield, friendly environment and high anti-disaster capability and can adopt to the modern agricultural production and operation modes by the rural land reclamation and readjustment.
Protect and improve quality of farmland

Farmland protection in Beijing and
Tianjin is difficult to implement thoroughly, hence, the farmlands have been nearly depleted [68,69].
Zone I comprised the deciduous broad-leaved forest ecoregion of the Yanshan-Taihang Mountain, which is the eco-security shield in the northwest of the BTH region. With high mountains, steep slopes and mountainous forestland, the ecosystem in Zone I has degraded seriously over the last few decades, exhibiting ecosystem homogenization, water conservation capacity decrease and soil erosion aggravation. The area with low vegetation coverage shows land desertification and is a high-incidence area for debris flow and geological collapse. Our results demonstrated a negative spatial relationship between grain production and soil retention in Zone I. The grain production in this zone is not effective for soil retention. By contrast, the forestland and grassland can provide high-quality soil retention service. Therefore, the government must continue to implement the Sloping Land Conversion Program (SLCP) and strictly prohibit agricultural reclamation or deforestation on the high-slope areas to compensate for the encroachment in plains via the Policy of Dynamic Equilibrium of the total Cultivated Land (PDDCL). Moreover, certain engineering measures, such as the slope protection project and torrent control works, must be further improved to alleviate water and soil loss. For places with low vegetation and high risk of geologic disasters, livestock grazing and fuel gathering should be stopped to facilitate afforestation.
Zone II covered the agricultural ecoregion of the Northern China Plain and was in the southeastern plain of the BTH region. This area has a flat topography with abundant concentrated and continuous farmland and fertile soil, thereby contributing to its capability to offer grain production and carbon sequestration service. To prevent the encroachment of high-quality farmland, the Protection of Basic Agricultural Land (PBAL) should be implemented strictly and regional development must utilize land use intensively. Newly added settlements should occupy the idle land or the low-quality and scattered farmland. This region should focus on developing modern agriculture for which the High-standard Primary Farmland Construction Project (HSPFCP) should be generalized to proceed with rural land reclamation and readjustment. Our result indicated that farmland with high grain production decreases water yield. The farmland in Zone II should fallow moderately to recharge the groundwater.
Zone III comprised an urban and suburban agricultural ecoregion that included two sub-ecological regions of central urban zones and suburban agricultural ecological zones. This zone covered super urban agglomeration areas, including Beijing and Tianjin. The urban environment is deteriorating and appears as an urban heat island with degraded air quality [70]. Our results demonstrated that the farmland in this region had obviously decreased. The farmland had strong capacities of water yield and carbon sequestration and grain production in Zone III had a positive spatial correlation with carbon sequestration and water yield. As the main land cover type in the BTH region, the farmland not only improves the services of grain production, water yield and carbon sequestration but can also serve as an urban green landscape. Therefore, the laws of farmland protection and regulations should be made consistent, transparent and nondiscriminatory to prevent settlement encroachment. This work recommends that the cultivation of farmland be divided into three circle-layers, namely, circles of urban landscape agriculture, eco-leisure agriculture and high-efficiency facility agriculture.
Zone IV was a typical steppe ecoregion of the Central-Eastern Inner Mongolia Plateau that was situated mainly in the northwest of the BTH region. The grassland had been seriously degraded because of unreasonable human exploitation, exhibiting a decrease in grass cover, land desertification and depletion of soil nutrients. The results of this study showed a negative spatial relationship and weak trade-off between grain production and soil retention in Zone IV. To prevent the further deterioration of grassland, the government should control the extent of agricultural reclamation by prohibiting excessive cultivation and implement conservation tillage to prevent soil wind erosion and desertification. In addition, pasture management needs to ban grazing and implement the strategy of rotation grazing.

Limitations
Uncertainties remain in the ES assessment that is based on models and data in this study. First, various ES are generated in different spatial scales. In our research, the spatial resolution of 1 km was helpful in learning the dynamic change in the research area at the regional scale; however, this system inevitably concealed the land use changes in areas under 1 km 2 . Second, soil retention had not been calibrated by observation data due to limited available data, although the results were consistent with those obtained by Li et al. (2014) and Xu et al. (2011) [53,71]. Third, only the grain productive capacity of farmland was considered in the assessment of food provision service. Grassland provides fibrous feed for livestock and orchards offer fruits for people, thereby eventually providing food for human survival. In further research, the accounting scope should be further expanded to cover additional essential services. Fourth, our study analyzed the interaction of ES in four ecological function zones but spatial heterogeneity remains in each ecological-functional zone. Detailed decisions will be implemented on the basis of the natural environment and socioeconomic status within each region in the future. Finally, the spatial distribution of human habitation was mismatched with that of the ecosystem, thereby resulting in a mismatch between the demand for and supply of ES. The current decisions provided remain based on the ES supply. If possible, a policy of regulating ES for human needs that is based on human settlement distribution will be formulated by the researchers.

Conclusions
In this study, four ES were quantified and mapped explicitly spatially with biophysical models. The spatial correction, trade-offs and synergies among multiple ES were also analyzed in four ecological-functional regions to help determine mechanisms of ES as bases for effective decisions for sustainable regional development. The main conclusions are as follows.

1.
The spatial distribution of water yield was relatively homogeneous but slightly high in Zones I, II and III and all increased generally each year. Sediment retention was high in Zones I and IV and low in Zones II and III. These values all gradually increased every year. The distribution of carbon sequestration was similar to that of soil retention but a downtrend occurred in the entire region rather than in Zone IV. High-quality farmland with high grain production was mainly distributed in Zones II and III and those in Zone III indicated a significant decline.

2.
Spearman correlation analysis indicated a positive spatial relationship between the status quo of water yield and soil retention/grain production capability/carbon sequestration, soil retention and carbon sequestration in the whole region. Furthermore, a positive relationship was observed between grain production capability and carbon sequestration in Zones II and III but a negative one was observed in Zones I and IV. A negative correlation was observed between grain production capability and soil retention in Zones I and IV. Regarding the relationship between trade-off and synergy, our results indicated synergies between soil retention and water yield/carbon and trade-offs between water yield and carbon sequestration/grain production in the whole region. In addition, grain production and carbon sequestration had synergy in Zones I and IV but trade-offs in Zone II and III. A weak trade-off was observed between grain production and soil retention only in Zone IV.

3.
The results of the transfer matrix demonstrated that encroachment on farmland occurred in the BTH region in 2000-2015, except for the grassland occupation by farmland in Zone IV. The government is suggested to implement sequentially the SLCP and prohibit reclamation and deforestation in Zone I. The authors emphasize the protection of high-quality farmland in Zone II to render it the "rice bag" of the BTH region. Furthermore, the farmland should be made multi-functional and its encroachment in Zone III should be prevented. Finally, the farmland should be returned for grass and grazing should be controlled in Zone IV to restore grassland ecology.