Analyzing Three-Decadal Patterns of Land Use/Land Cover Change and Regional Ecosystem Services at the Landscape Level: Case Study of Two Coastal Metropolitan Regions, Eastern China

Abstract

Rapid urbanization, land scarcity, and accompanying ecological deterioration in China have received growing attention. In this paper, two fast-growing metropolitan regions, Greater Shanghai and Greater Hangzhou, were selected as case studies to quantify the impact of land use/land cover (LULC) change on regional ecosystem services value (ESV) at the landscape scale since the late 1970s. The results show that in both regions, dramatic LULC change, especially recent land development at the urban fringes, led to a steady decline in the available area of productive agricultural land, natural land and semi-natural land. This consequently caused remarkable landscape fragmentation along the urban-rural gradient as measured by five class-level landscape metrics. It was estimated that in Greater Shanghai, regulating, supporting, provisioning, and cultural ESVs decreased by 32.05%, 17.89%, 53.72%, and 17.06%, respectively. In Greater Hangzhou, these values decreased by 27.82%, 23.86%, 28.62%, and 22.85%, respectively. In addition, the relationship is quantified between zonal buffer-based ESV and class-level landscape metrics. Further analysis shows that spatiotemporal patterns of zonal ESVs along the urban-rural gradient in these two regions exhibited unbalanced patterns of ecological services delivery.

1. Introduction

Human induced land use/land cover (LULC) change, particularly accelerated urbanization, has played a key role in the transformation of landscapes and ecosystems worldwide. One of the issues of concern underlying LULC change is consequent landscape fragmentation, which changes structure and pattern of ecosystems and decreases ecosystems’ functions through a wide range of ecological functions and processes. As a result, landscape fragmentation not only affects ecosystems provisioning to meet the basic demands of consumption by society, but also impairs ecosystems’ buffering capacity for human and natural communities, including flood protection, climate regulation, and the control of diseases and pests [1,2,3,4,5,6,7]. Therefore, when revisiting the human-nature relationship and finding the solution for sustaining ecosystems’ functions, trends in landscape fragmentation and declining ecosystems’ functions have increasingly attracted attention. The dynamics of landscape patterns can be measured to indirectly depict LULC change and ecological consequences [8,9,10,11]. The identification and measurement of varying ecosystem services linked to changing landscape patterns may help quantify the environmental cost-benefit of different land planning decisions, and thus allow decision-makers to better understand different trade-offs for efficient ecosystem management [12,13,14]. Therefore, in this context, incorporating dynamics of landscape and ecosystem services into land use planning provides a practicable way for decision-makers to efficiently manage ecosystems and land use, especially in setting program priorities, choosing among environmental options, and communicating the importance of their actions to the public [15]. Unfortunately, published articles that quantified the relationship between landscape patterns and ecosystem services were relatively scarce at the time of this research, though some case studies focusing on landscape pattern and ecosystem services were carried out [16,17,18,19,20,21,22].

China has experienced unprecedented urbanization since the late 1970s, as evidenced by the booming megacities along the coastal economic regions. Uncontrolled urban expansion, explosive population growth, loss of arable land, and pronounced environmental problems have bottlenecked this country’s sustainable development [23,24,25,26,27,28,29,30]. Therefore, practicable policies for land use management and urban planning, which address these emerging environmental challenges and further aim at developing rational solutions for sustainability in human-dominated ecosystems, are urgently needed for a rapidly urbanizing China.

This study focuses on Greater Shanghai and Greater Hangzhou, which are currently the largest and second largest metropolitan regions in the Yangtze River basin, respectively [31]. Since the 1990s they have both undergone population booms and rapid urbanization, as witnessed by agricultural land in the urban periphery being remarkably converted to urban settlements, infrastructure, and industrial parks. Local governments can maintain a land bank through land reclamation and creation of new land by enclosing tidal areas. However, it seems that there is no long term solution embodied in official policies regarding land use choice and urban design, nor is baseline guidance issued to balance the conflict between insufficient buildable land and land development. Therefore, to address these challenges, the objectives of this study are: (1) understanding rapid urbanization and the consequences of LULC change; (2) quantitatively assessing the status of regional ecosystem services at the landscape level. Our goal is to help promote ecosystem resilience and the maintenance of sustainable ecosystem provisioning in fast-growing metropolitan regions in China and the other developing countries.

2. Materials and Methods

2.1. Study Areas

Figure 1 shows the location of Greater Shanghai and Greater Hangzhou.

Greater Shanghai has been the economic center of China since the 1900s. It is located between latitudes 31°32′ N–31°27′ N and longitudes 120°52′ E–121°45′ E. This region has a northern subtropical monsoon climate, with an average annual temperature of approximately 15 °C. Temperatures average 28 °C in the summer and 4 °C in the winter. The average annual precipitation is approximately 1000–1200 mm, with approximately 60% of the rainfall being received during spring and autumn. Topographically, the region is mainly located on an alluvial terrace, with an average elevation of 4 meters. The Huangpu River is the major river encompassing the region. This region covers an area of 10,000 km², with 18.18 million residents [32].

Greater Hangzhou is the capital city of Zhejiang province, with a recorded history of approximately 2200 years. This city is well known for traditional culture, agricultural products, natural landscape, historical resorts, and relics. It is situated between latitudes 29°50′ N–30°32′ N and longitudes 119°41′ E–120°43′ E. This region has a northern subtropical monsoon climate with an annual temperature ranging 15.7–17.2 °C. Annual precipitation varies from 1352 to 1601.7 mm, of which approximately 80% is received during spring and autumn. The Qiantang River is the major river encompassing eastern Greater Hangzhou. Generally, most of this region is located on a flood plain, with a surface elevation ranging from 2 to 10 m, while hilly and mountainous parts account for 28.8% of the region. For a long time, its metropolitan statistical area was limited within the city proper. However, its administrative jurisdiction was redefined in 2000 by merging neighboring Yuhang County and Xiaoshan County. At present, it covers an area of approximately 3320 km², with 4.99 million residents [33].

2.2. Data Sources

In this study, both regional socioeconomic data and multi-temporal remotely sensed data were used. Regional socioeconomic data were extracted from the statistical yearbooks [32,33,34]. The multi-temporal remotely sensed data, including Landsat Multispectral Scanner (MSS), Thematic Mapper (TM), and Enhanced Thematic Mapper Plus (ETM+) imagery spanning the study period were also used (see

Table A1. Satellite images used in this study.

Study Area	On-Board Sensor	Path/Row	Acquisition Date (DD-MM-YY)	Resolution (m)
Greater Shanghai	MSS	127/038, 127/039	4-8-1979	60
	TM	118/038, 118/039	18-5-1987	30
	TM	118/038, 118/039	11-4-1997	30
Greater Hangzhou	TM	118/038, 118/039	24-3-2008	30
	MSS	128/039	5-7-1978	60
	TM	119/039	23-7-1991	30
	ETM+	119/039	11-10-2000	30
	ETM+	119/039	24-4-2008	30

Note: This table only gives the spatial resolution of reflective bands.

).

2.3. Satellite Imagery Preprocessing, Classification, Accuracy Assessment, and Post-Classification

Eight Landsat MSS/TM/ETM+ satellite images, which were clear or nearly free of cloud contamination, were used for this study (see Table A1 in Appendix A). For Landsat ETM+ images, one SLC-off scene dated 24 April 2008 was repaired with a gap-filled method known as self-adaptive local regression model for multi-temporal imagery [35]. The images were geometrically rectified and georeferenced to the WGS84 UTM map projection system, using the 1:250,000 digitalized administrative maps of greater Shanghai and greater Hangzhou. The georeferenced images were combined to produce false-color images for visual interpretation. Based on our prior knowledge from previous studies and field surveys, a local land cover classification scheme including developed land, cropland, forest, shrub, water, tidal land, and bare land (see

Table A2. Land cover classification scheme in this study.

Land Cover Type	Description
Developed land	Visually detectable urban and rural settlements, commercial areas, transportation lines, and industrial parks.
Cropland	Paddy fields, fallow lands after harvest, and dry lands.
Forest	Natural and artificial woodlands.
Shrub	Wild scrubland and forest nurseries.
Water	Rivers, creeks, reservoirs, lakes, fishponds, and dikes.
Tidal land	Sandy flat periodically inundated by tides.
Bare land	Bare rocks, gravel pits, quarries, mines, permanently enclosed tidal land, and vacant land after clearing vegetation for urban development.

) was defined, according to guidance released by China National Committee of Agricultural Division [36] for classification of land cover. Subsequently, the ‘supervised signature extraction with maximum likelihood’ algorithm was used to perform classification of land covers. For each image, we chose at least 100 training sites to ensure that all spectral classes representing each LULC category were adequately captured in the training statistics. Imagery preprocessing and classification was performed with GEOSTAR 3.0^® image processing software.

We performed an accuracy assessment by undertaking the following steps. First, the key ancillary data, including the national land use dataset (1990, 1995, 2000, 2005, 2008), local historical aerial photos (2000, 2003, 2007), and 1:250,000 digitized LULC maps (1988, 1991, 1996, 2000) [37,38,39] were used to collect the reference data. Second, on each classified map, for each land cover class, 50 sample sites were selected, using the random stratified method. Thus, for each classified map the accuracy was assessed with a total of 250 samples. Third, reference data and false-color maps were combined with the classification maps to improve the overall accuracy of the classified images [40]. The User’s accuracy, Producer’s accuracy, overall accuracy, and Kappa statistic of LULC maps of the study areas were shown in

Table 1. Accuracy assessment of LULC maps of Greater Shanghai.

LULC Type	Developed Land		Cropland		Forest		Water		Tidal Land		Bare Land
	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)
Developed land	75.00–87.10	62.50–71.15	0	0	0	0	0–7.32	0–7.67	3.57–5.56	4.55–6.67	14.63–25.00	15.79–22.92
Cropland	0	0	62.5–80.49	76.92–100.00	10.20–35.42	16.67–32.69	0	0	0–2.44	0–2.70	0–10.20	0–10.42
Forest	0	0	0–20.45	0–23.08	76.09–100.00	67.31–83.33	0	0	0–4.35	0–4.55	0	0
Water	0–2.33	0–2.86	0	0	0	0	77.78–86.05	84.85–100.00	11.63–23.91	11.36–29.73	0	0
Tidal land	2.77–8.00	2.44–4.65	0	0	0–4.65	0–6.67	0–11.63	0–15.55	67.44–80.00	54.05–75.00	8.11–12.00	5.77–13.16
Bare land	2.44–13.89	2.13–12.20	0	0	0	0	0	0	3.23–12.5	2.27–11.11	75.00–90.24	62.50–71.15
UA (%)	76.88–82.81
PA (%)	76.93–83.48
OA (%)	76.40–84.20
Kappa statistic	0.72–0.79

Note: R (%) and C (%) denote percent of Row and percent of Column, respectively. UA, PA, and OA denote User’s accuracy, Producer’s accuracy, and Overall accuracy, respectively. All of the data in this table mean the intervals of statistics across the study period.

and

Table 2. Accuracy assessment of LULC maps of Greater Hangzhou.

LULC Type	Developed Land		Cropland		Forest		Shrub		Water		Tidal Land		Bare Land
	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)	R (%)	C (%)
Developed land	68.90–84.80	83.09–92.01	0–1.71	0.76–4.29	0–1.78	0–2.33	0–0.67	0.65–4.24	0–4.11	0–2.89	0–5.15	0–8.21	0–4.90	0–5.32
Cropland	1.75–11.59	0–4.66	69.90–95.81	83.75–92.31	2.79–13.10	3.61–11.44	0–11.26	0–9.71	0.56–9.58	0.32–2.35	2.11–13.16	0–3.17	0–6.80	0–23.61
Forest	0–4.68	0–2.40	1.25–5.97	1.22–4.29	75.86–96.17	80.47–92.88	0–3.36	0–7.79	0–1.58	0–2.31	0–1.58	0	0–0.91	0–0.69
Shrub	0.58–9.76	0–0.65	0–5.24	0–3.84	0–10.34	0–4.03	81.46–98.34	79.87–96.74	0–3.42	0–1.54	0–3.68	0–6.72	0–5.88	0–0.70
Water	0–5.49	0–9.56	0.19–1.73	0.23–5.95	0–1.33	0–3.25	0–1.32	0–3.25	81.57–96.05	73.08–96.38	0–30.00	1.14–9.51	0–2.94	0–4.90
Tidal land	0–6.43	0–4.41	0–0.76	0.77–5.73	0	0–1.49	0–5.96	0–4.00	0.49–10.57	0–21.79	49.12–96.32	74.63–95.43	0–2.90	0–31.47
Bare land	0–10.55	0–3.27	0–9.43	0–3.44	0–1.63	0–1.49	0–1.34	0–1.63	0–1.72	0–0.64	0–13.24	0–4.48	83.33–100.00	59.44–100.00
UA (%)	83.19–88.31
PA (%)	82.54–86.99
OA (%)	80.70–88.56
Kappa statistic	0.78–0.85

Note: R (%) and C (%) denote percent of Row and percent of Column, respectively. UA, PA, and OA denote User’s accuracy, Producer’s accuracy, and Overall accuracy, respectively. All of the data in this table mean the intervals of statistics across the study period.

. As shown, the overall accuracy for two regions is acceptable, according to the values recommended by Jassen et al. [41]. Finally, LULC change detection was performed through post-classification comparisons [42]. Land use change matrices quantitatively representing the overall LULC change in the study areas were established. Thereafter, based on the gain-loss within and between LULC categories of land use change matrix, LULC change maps were produced.

2.4. Computation of Class-Level Metrics for Measuring Landscape Fragmentation

Our interest was mainly focused on measuring landscape fragmentation under the pressure of human activities, especially during the recent rapid transition of land use attributable to urban growth and socioeconomic development. Five class-level metrics, including patch density (PD), percentage of landscape (PLAND), mean patch size (MPS), largest patch index (LPI), and landscape shape index (LSI), were employed to measure average fragmentation for the whole landscape (See Table A2). To depict spatiotemporal patterns in the landscape, a series of concentric zones with different buffering distances from the city centers of the two metropolises were adopted. Note that in this study, the choice of zonal width or buffer distances from the city centers is largely based on our deep understanding of the urban growth patterns of the study areas [43,44]. Initially, we set a variety of concentric buffers with even distance intervals ranging from 1 and 10 km, aiming to better capture variation in landscape structures along the urban-rural gradients. Narrower zonal buffers may provide more detailed information, but they may also present redundant results [10]. Alternatively, we tried to merge the zonal buffers with uneven distance intervals and found that they showed robust trends corresponding to an increase in distance intervals. Take Greater Hangzhou as an example, the distance intervals ranging from 6 to 15 km are suitable in the urban fringe, distance intervals ranging from 15 to 25 km are acceptable in the exurban and rural areas given the class-level landscape metrics changed little within these intervals. These findings provided the baselines for drawing the zonal buffers. Therefore, for Greater Shanghai, zonal buffers were drawn at 0–6 km, 6–12 km,12–21 km, 21–35 km, and >35 km (35–79 km) from the city center. For Greater Hangzhou, zonal buffers were drawn at 0–3 km, 3–6 km, 6–15 km, 15–25 km, and >25 km (25–40 km) from the city center (see Figure A1 and

Table A4. Description of zonal buffers from the city cores of Greater Shanghai and Greater Hangzhou.

Region	Buffers Distance (km)	Synopotical Description
Greater Shanghai	0–6	The city core of downtown Shanghai.
	6–12	The newly in-filling urban area between the inner and outer rings.
	12–21	The urban finge with intensive settlements and industrial parks.
	21–35	The rapidly urbaning areas with intensive settlements, industrial parks, harbors, and airport
	>35	The low-density developed rural areas with sparsely distributed towns and villages. Aside from some settlements and industrial parks, this zonal buffer is characterized with cropland and tidal land.
Greater Hangzhou	0–3	The city core of downtown Hangzhou.
	3–6	The newly in-filling urban area between the city core and neiboring towns.
	6–15	The rapidly urbanizing areas expanding eastward and southward between the city core and well-developed urban area of Xiaoshan district.
	15–25	The mixing middle-density and low-density developed rural areas with sparsely distributed towns and villages. Aside from the well-developed urban area of Yuhang district, this zonal buffer is characterzied with hilly terrain, cropland, and river network.
	>25	The low-density developed rural areas with sparsely distributed towns and villages. This zonal buffer is characterzied with hilly and mountaineous terrain and tidal land.

). Subsequently, for each zonal buffer, all of the previously described metrics were computed to depict fragmentation distribution along the urban-rural gradient to capture landscape configuration by using Fragstats 4.2 [45].

2.5. Computation of Regional ESV

Quantitatively measuring ecosystem services helps to better understand the importance of ecosystem function. However, how to measure ecosystem services has been a hotly-debated issue. Since Costanza et al. [46] used restoration cost to derive economic values for the ecosystem services of global biomes, this method has received special attention. However, criticism arose because simply assigning the ESV coefficients for a specific land cover with its global average may cause biased results. To fix such a problem, Xie et al. [47] developed an enhanced method for valuing China’s terrestrial ecosystem services by surveying 200 Chinese ecologists, who scored the per ha ESV coefficients for typical terrestrial ecosystems. By comparing the surveyed results and Costanza et al. ‘s assumption, both of the biased ESVs for wetland (overestimated) and cropland (underestimated) were adjusted when downscaling the global level to Chinese localized context. Therefore, Xie et al.’s enhanced method is considered more practicable and has thus been widely adopted [11,23,24,25,26,27]. Therefore, we adopted this enhanced method in this study by modifying the ESV coefficients of soil, water and tidal land. According to four categories of ecosystem services delineated by MEA [5], aforementioned ESV was regrouped into four categories as with regulating (gas regulation and climate regulation), supporting (soil formation and retention, waste purification, and biodiversity protection), provisioning (water supply, food production, and raw material), and cultural (recreation and culture). The annual ESV of each land cover category is shown in

Table A3. Formula and description of class-level landscape metrics.

Formula	Description
$PD=\frac{n_i}{A}\times 10000$	where ni is counts of land cover patch (class) i, A is total area of all patches (m2).
$PLAND=P_i=\frac{{{\displaystyle \sum }}_{i=1}^j{a}_{ij}}{A}\times 100\%$	where Pi is percentage of patch type (class) i within the landscape, aij is area (m2) of patch ij, and A is area of all patches (m2).
$LPI=\frac{\begin{array}{c}j\\ \underset{i=1}{Max}\left(a_{ij}\right)\end{array}}{A}$	where aij is area (m2) of patch ij and A is total area of all patches (m2).
$LSI=\frac{e_i}{\min e_i}$	where ei is class i’s total length of edge for given grids; min ei is class i’s minimum total length of edge.
$MPS=\frac{A_i}{N_i}$	where Ai is class i’s area (m2), and Ni is number of class i.

.

The estimated regional ESV was calculated as follows,

$$ESV={\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{j=1}^n{A}_i\times V{C}_{ij}$$

(1)

where A_i is area (ha) of land cover i, VC_ij is the value coefficient of ecosystem service function for type j (RMB Yuan/ha) combined with land cover i, m is the number of land covers, and n is the number of ecosystem service functions, respectively (see

Table A5. Annual ESV of each land cover category (RMB Yuan/ha).

Ecosystem Services Category	Ecosystem Service Functions	Land Use Category
		Forest	Cropland	Water	Shrub	Bare Land	Tidal Land
Regulating	Gas regulation	3097.00	442.40	0.00	1769.70	0.00	0.00
	Climate regulation	2389.10	787.50	407.00	1588.30	0.00	203.50
Supporting	Soil formation and retention	3450.90	1291.90	8.80	2371.40	17.70	13.30
	Waste purification	1159.20	1451.20	16086.60	1287.20	8.80	8047.70
	Biodiversity protection	2884.60	628.20	2203.30	1756.40	300.80	1252.10
Provisioning	Water supply	2831.50	530.90	18,033.20	1681.20	26.50	9029.90
	Food production	88.50	884.90	88.50	177.00	8.80	48.70
	Raw material	2300.60	88.50	8.80	1194.60	0.00	4.40
Cultural	Recreation and culture	1132.60	8.80	3840.20	570.70	8.80	1924.50

in Appendix A).

2.6. Statistical Analysis

To quantify the relationship between dependent variable (zonal buffer-based ESV) and independent variables (class-level landscape metrics) along urban-rural gradient, the multiple linear regression model is a useful tool. However, it is very important to diagnose co-linearity effects prior to establishing reasonable regressions, given that significant correlation may exist between variables. The stepwise regression model was employed to exclude co-linearity effects between variables because of its popularity in statistical packages [48]. In this study, according to result of normal distribution test, zonal annual ESV was found to be skewed. Then, Box-Cox transformation for zonal annual ESV was performed and a natural log-transformation was adopted. Thus, the stepwise regression model can be written as follows,

$$\text{Ln}Y=\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2+\dots {\beta }_n{X}_n+\epsilon$$

(2)

where Y is zonal annual ESV, α is the constant, β_i (i = 1, 2, 3, n) are the partial regression coefficients for significant independent variables (class-level landscape indices), and ε is random error.

This model starts as an empty model and then adds or removes a variable for each step, according to the user-defined criterion for introducing a new variable to the model or removing a variable from the model. Herein, threshold values for alpha-to-enter and alpha-to-remove were both set at 0.10 by comparing criteria ranging between 0.10 and 0.20 for introducing or removing a variable at each step to ensure that the model contained the most significant potential variables. All of these statistical processes were performed with Data Processing System (DPS) version 12.05 [49].

3. Results

3.1. LULC Change

Figure 2 and Figure 3 generally show LULC dynamics of the study areas. As shown, in Greater Shanghai, cropland and water are dominant land covers, followed by developed land. At present they account for 91.88% of regional land covers. Forest, tidal land, and bare land only account for a very small proportion of regional land cover. It was found that developed land grew by 1230 km² in 1979–2008, showing an exponential growth pattern. The continuous growth of developed land mainly occurred within the 6–12 km, 12–21 km, and 21–35 km zonal buffers, where the newly developed land should be responsible for approximately 67.35%, 9.12%, 4.80%, and 2.19% loss of cropland, forest, water, and tidal land, respectively. In Greater Hangzhou, cropland, developed land, and forest are the dominant land covers, which account for 78.12% of present land cover. It is noteworthy that developed land grew by 542.50 km² in 1978–2008. Further analysis shows that most of the presently developed land mainly occurred within the 3–6 km, 6–12 km, and 12–25 km zonal buffers, where presently developed land should be responsible for approximately 38.95%, 24.40%, 9.89%, and 7.14% loss of cropland, shrub, forest, and water, respectively.

3.2. Variation of Class-Level Landscape Patterns

Figure 4 shows the varying class-level landscape patterns in zonal buffers of Greater Shanghai. The overall declining trends in PLAND and LPI for developed land from the city core to the exurban and rural areas were observed. As shown, in 1979–1987 the spatial extent of developed land, in particular urban build-up land, was limited to within a 0–6 km distance from the city core. The 6–12 km distance from the city core was typical urban periphery, which was characterized by a mixture of urban and rural landscapes. During this period the overall changes in PLAND and LPI for the other land covers such as cropland and water show stability, with only slight fluctuations. In 1987–1997, the spatial extent of developed land within the 0–6 km and 6–12 km distances from the city core increased remarkably. It is noteworthy that within the 6–12 km distance from the city core the relative growth rate of developed land was much higher. In contrast, cropland falling within this zonal buffer showed an increase in PD and a remarkable decline in PLAND and LPI. The national strategy of ‘opening Pudong new area and developing Shanghai’, which initiated in 1990s and triggered rapid urbanization, should be responsible for more fragmented cropland of this zonal buffer due to intensive urban encroachment during this period. In 1997–2008, the spatial extent of developed land within the 0–6 km and 6–12 km distances from the city core remained relatively stable in 1997–2008. Significant growth of developed land occurred in the 12–21 km and 21–35 km distances from the city core. Generally, judged by distance-based varying landscape metrics of all land covers, it can be concluded that outward expansion of developed land occurring in different stages resulted in more compact urban areas within the city proper and less fragmented urbanizing areas (the 6–12 km and 12–21 km distances from the city core), as evidenced with more fragmented semi-natural and natural landscapes along the urban-rural gradient.

Figure 5 shows the varying class-level landscape patterns in Greater Hangzhou. The overall decreases in PLAND (percentage of landscape) and LPI (largest patch index) for developed land with increasing buffer distances from the city core were detected in 1978–2008. In contrast, changes in PLAND and LPI for the other land covers were relatively stable across the study period. There was remarkable fluctuation of developed land within different zonal buffers. As shown, developed land within the 0–3 km buffer distance initially increased in 1978–1991, but a consecutive slowdown occurred across the period 1991–2008. Generally, developed land within the 3–6 km buffer distance increased remarkably across all stages. This indicated that the in-filling development of intensively developed land was predominant. The much lower LSI for all land covers within the 0–3 km and 3–6 km buffer distances indicate substantial growth of developed land throughout the study period, with replacement of the other land covers. It is also noteworthy that significant growth in developed land within the 6–15 km and 15–25 km buffer distances occurred from 2000 onward. In 2000, Greater Hangzhou’s administrative jurisdiction was reshaped by merging the former Yuhang County and Xiaoshan County. This consequently accelerated the rapid expansion of urbanized and urbanizing areas. The increasing PLAND, LPI, and LSI as well as decreasing PD of developed land were detected within these two buffer distances in 2000–2008. For the 15–25 km and >25 km buffer distances, there were upward trends in LSI for developed land, cropland, shrub, water, and bare land across the study period. In contrast, the LSI of forest increased in 1979–1991 but decreased in 2000–2008. This may be explained by the recent recovery of wild and semi-natural vegetation due to the local government’s policy for restoration and conservation of forests.

3.3. Variation in Spatiotemporal Pattern of ESV

Table 3. Annual ESV in Greater Shanghai and Greater Hangzhou (unit: million RMB Yuan).

Ecosystem Service Category	Greater Shanghai				Greater Hangzhou
	1979	1987	1997	2008	1978	1991	2000	2008
Regulating	1101.19	916.58	844.39	747.97	780.93	658.26	616.03	563.7
Supporting	9409.11	8514.57	8548.44	7725.88	1894.61	1620.96	1527.04	1442.6
Provisioning	8264.66	7437.53	7489.05	6844.91	1343.87	1111.41	1055.24	1043.8
Cultural	1595.93	1419.73	1425.49	1323.73	246.69	189.13	180.86	190.31
Sum	21,966.82	19,708.14	19,732.86	17,966.22	4512.79	3768.89	3560.03	3430.7

shows an overall decline in ESVs for each specific ecosystem function of the study areas. In Greater Shanghai, there is an uneven pattern in the change rates of ESVs for each specific ecosystem function throughout the study period. Obviously, the overall change rates in the ESVs of raw material and gas regulation declined remarkably compared to those of the other specific ecosystem functions. In contrast, Greater Hangzhou exhibits a relatively even pattern for change rates of ESVs for each specific ecosystem function.

Table 4. Variation of zonal ESV (million RMB Yuan) in Greater Shanghai.

Zonal Buffer	Year				Stage Change Rate
	1979	1987	1997	2008	1979–1987	1987–1997	1997–2008	1978–2008
0–6 km	59.84	39.87	35.25	41.05	−33.36%	−11.61%	16.48%	−31.39%
6–12 km	236.82	164.71	98.65	109.82	−30.45%	−40.11%	11.32%	−53.63%
12–21 km	780.60	650.36	660.73	484.66	−16.69%	1.60%	−26.65%	−37.91%
21–35 km	3748.93	3522.41	3682.43	3311.89	−6.04%	4.54%	−10.06%	−11.66%
>35 km	15,544.69	13,911.03	13,830.31	12,695.06	−10.51%	−0.58%	−8.21%	−18.33%

indicates the variation in zonal ESV with distance from the city core for Greater Shanghai. Similar to Table 4, the overall declining trends in zonal ESVs were detected in several stages. However, the maximum decline in zonal ESVs occurred within the 6–12 km buffer distance, followed by those within the 12–21 km and 0–6 km buffer distances. The minimum decline in zonal ESV occurred within the 21–35 km buffer distance, followed by another zonal ESV within the >35 km buffer distance.

Table 5. Variation of the zonal ESV (million RMB Yuan) in Greater Hangzhou.

Zonal Buffer	Year				Stage Change Rate
	1978	1991	2000	2008	1979–1991	1991–2000	2000–2008	1978–2008
0–3 km	16.49	6.79	8.93	13.90	−58.82%	31.52%	55.66%	−15.71%
3–6 km	92.68	62.35	55.42	67.41	−32.73%	−11.11%	21.63%	−27.27%
6–15 km	746.42	671.23	585.21	524.87	−10.07%	−12.82%	−10.31%	−29.68%
15–25 km	1287.31	1191.25	1115.90	989.62	−7.46%	−6.33%	−11.32%	−23.12%
>25 km	2123.18	1931.19	1965.73	1912.66	−9.04%	1.79%	−2.70%	−9.92%

shows the variation in the zonal buffer-based ESV with differing distances from the city core of Greater Hangzhou. Overall, declining trends in zonal buffer-based ESVs were detected, though some zonal buffer-based ESVs within the 0–3 km and 3–6 km buffer distances decreased during 1978 and 2000 and then rebounded in several stages. As a whole, the maximum decline in zonal buffer-based ESVs occurred within the 6–15 km distance from the city core, followed by those occurred within the 3–6 km, 15–25 km, and 0–3 km zonal buffers. In contrast, the minimum decline in zonal buffer-based ESV occurred within the >25 km buffer distance.

Figure 6a,b generally show a power growth trend in the zonal mean ESV and its percentage of the study areas. As shown, in Greater Shanghai the zonal mean ESV within the 0–6 km, 6–12 km, 12–21 km, 21–35 km, and >35 km buffer distances accounted for 0.24%, 0.83%, 3.50%, 19.38%, and 76.05% of regional mean ESV, respectively. Consequently, zonal ESVs within the 0–21 km and 0–35 km buffer distances only accounted for 4.53% and 23.95% of regional mean ESV, respectively. In contrast, the 21–35 km and >35 km zonal buffers are typically exurban and rural areas with much higher ESV. In Greater Hangzhou the zonal mean ESV within the 0–3 km, 3–6 km, 6–15 km, 15–25 km, and >25 km buffer distances accounted for 0.39%, 2.17%, 17.50%, 30.18%, and 49.77% of regional mean ESV, respectively.

Based on these results, further analysis combining regional population growth and spatial agglomeration should produce more insight into the equity and rationality of spatial allocations of regional ESV. In Greater Shanghai, the spatial extent within the 0–35 km buffer distance is characterized by highly urbanized and urbanizing areas with dense population. Approximately 80% of regional total population resides in this extent. Accordingly, zonal buffer-based ESVs within this extent only accounted for 33.95% of the regional total ESV. In Greater Hangzhou, approximately 70.37% of regional total population resides within the 0–15 km buffer distance from the city proper, of which zonal buffer-based ESVs only accounted for 20.06% of regional total ESV. Additionally, Figure 7 indicates that over the study period the total population in Greater Shanghai has increased by 66.09% (approximately 7.51 million) and in Greater Hangzhou the total population has increased by 46.06% (approximately 1.35 million). Accordingly, ESV per capita in Greater Shanghai and Greater Hangzhou decreased by 50.81% and 46.60%, respectively. Thus, when measured with ESV per capita, the status of ecosystem services for Greater Shanghai and Greater Hangzhou degraded remarkably.

3.4. Relationship between the Landscape Pattern and Allocation of Esv

Table 6. Multi-linear regression models quantifying the relationship between zonal buffer-based LnESV and class-level landscape metrics.

Multi-Linear Regression Model	R2
Ln ESVSH = 6.142 − 0.042 LPID + 0.001MPSC + 0.001 MPSF + 0.007 PLANDW+ 0.002 LSIT	0.993 **	(3)
Ln ESVHZ = 2.355 − 0.011 LPID + 0.013 LSIC + 0.036 PLANDF − 0.002 MPSB + 0.178 MPST	0.990 **	(4)

Note: For the dependent variables, subscripts SH and HZ denote Greater Shanghai and Greater Hangzhou, respectively. For the independent variables, subscripts D, C, F, W, B, and T denote developed land, cropland, forest, water, bare land, and tidal land, respectively. ** denotes significant at p = 0.01.

shows the multi-linear regression models depicting the relationship between zonal buffer-based class-level landscape metrics and the allocation of ESV in Greater Shanghai and in Greater Hangzhou.

Equation (3) shows that in Greater Shanghai there are significant negative associations between log-transformed zonal ESV (Ln ESVSH) and LPI of developed land, whereas the MPS of cropland and forest, PLAND of water, and LSI of tidal land exhibit significant positive associations with log-transformed zonal ESV. Equation (4) shows that in Greater Hangzhou, there are significant negative associations between log-transformed zonal ESV (Ln ESVHZ) and the LPI of developed land and MPS of bare land. In contrast, the LSI of cropland, PLAND of forest, and MPS of tidal land show significant positive effects in determining variation of log-transformed zonal ESV. Overall, the partial coefficients of independent variables in abovementioned equations are somewhat different, indicating their role in determining spatiotemporal pattern of log-transformed zonal ESV. In Greater Shanghai, spatial extent and land development intensity is much higher than that in Greater Hangzhou, this can explain how developed land played the most pronounced role in decreasing log-transformed zonal ESV. However, the hilly and mountainous terrain as well as numerous cultural resorts in Greater Hangzhou shelter vast cropland, forests, shrubs, and tidal lands, which limit large-scale land development in the southwestern, northwestern, and eastern sub-regions. This can better explain how these semi-natural and natural landscape elements played the relative higher role in trading-off negative effect of urban expansion in Greater Hangzhou, though overall ESV in these two metropolitan regions decreased remarkably due to rapid urbanization. Moreover, observed together with Figure 3, Figure 4, Figure 5 and Figure 6 and Table A1, Table A2, and Table A3 apparently the trends in PLAND, LPI, and MPS of developed land are reverse to that in zonal ESV along the urban-rural gradient, while trends in LSI of all landscape elements are similar to that in zonal ESV. This result indicates that intensive human activities, especially urban encroachment caused more fragmented landscape and adversely affected the allocation of ecosystem services delivery along urban-rural gradients. In contrast, the rural areas away from the city proper were only slightly affected, as evidenced by relatively higher ESV and less fragmented landscape patterns.

4. Discussion

4.1. Revisiting the Cause-Effect Relationship between LULC Change, Landscape Fragmentation, and Ecosystems’ Functioning

Rapid expansion of Greater Shanghai and Greater Hangzhou exemplified the dilemma in China’s eastern developed regions. Both of these two fast-growing metropolitan regions changed from a typical single-core compact city to a mixture of compact and sprawling morphologies with multiple nuclei. The implementation of the government-oriented strategies for industrial restructuring and urban land development accelerated the spatial configuration between the downtowns and urban fringes [43,44,49,50]. As addressed, dramatic LULC change in Greater Shanghai and Greater Hangzhou, particularly recent expansion of developed land at the cost of semi-natural and natural lands, has resulted in a fragmented landscape and unbalanced patterns of ecological services delivery at a regional level, considering the relationship between spatial allocation of ecosystem services and human needs. Some specific ecosystem functions affected by LULC change and landscape fragmentation, such as pollination disturbance, habitat loss, landscape connectivity loss, which can be elaborated with one or several assumptions and empirical interpretations, such as pressures from land-use change and intensification [51,52], species’ functional traits influenced sensitivity to human-dominated land use [53], the isolating effects of different patterns [54], and human-induced shifts in the functional structure of biological communities with possible repercussion on important ecosystem functions and services [55]. However, when focusing on interpreting the mechanism underlying the cause-effect relationship between anthropogenic large-scale LULC change, landscape fragmentation, and influenced ecosystems’ functioning, the aforementioned assumptions and empirical interpretations may be problematic. Given the complexity of ecosystem processes and ecosystems’ response to human activities, it is still a challenge to well document the cause-effect relationship between LULC change, landscape fragmentation, and influenced ecosystems’ functioning. Discussion towards better understanding of the patterns and drivers of ecosystem processes and ecosystems’ response to human activities, should be put into a broad context of natural and socioeconomic factors.

4.2. Implications for Policies towards Sustainable Land Use and Ecosystem Management

Analysis of land use patterns and landscape dynamics affecting ecosystem services could provide a basis for informed decision making towards sustainable ecosystems management [11,19,21,56]. As addressed in Section 3, LULC change and land fragmentation measured with landscape metrics helped quantitatively illustrate the spatiotemporal pattern of deteriorated ecosystem services delivery during recent rural-urban transition. Such findings, which dynamically embodied the linkage between ESV cost/benefit supplying and LULC change [11,57,58,59], can not only be used for retrospective review of previous LULC change and its consequent ecological loss, but be borrowed as baseline for varying land development scenarios. These results can further base scientific soundness for local authorities’ decision-making process throughout making regulations, decisions, and policies towards effectively guiding land use and regulating the stakeholders’ land development intensity. However, it seemed that ongoing urban expansion, loss of non-urban land, and deterioration of regional ecosystem services in the study areas have not influenced local governments’ land use policy. There is official encouragement for cropland to be developed for modern manufacturing, real estate, and tertiary industries, given the increasing pressures from land scarcity and enterprises competing for land use [26,43,44,60,61,62,63]. The absence of systematic and practical approaches combining ESV cost/benefit supplying and LULC change in local governments’ decision-making should be responsible for the failure of land use policy. Therefore, to ensure sustainable land use and ecosystem services, landscape analysis, ecosystem services approach, and land use planning should be integrated to support governmental decision-making, which should focused on: (1) cost-benefit analyses of land development and ecological footprint; (2) trade-offs for balancing conflicts between land development and conservation of ecosystems; (3) shortcomings and limitations of current policy and regulation; and (4) to what extent specific policy alternatives address these issues as well as to envision space for further improvement. Only successfully achieving these goals, will these policy developments affect ecosystem conditions and determine the quality of provision of ecosystem services at multi-scales.

4.3. Limitation of This Study

In this study, Xie et al.’s [47] method for measuring ESV and class-level landscape metrics were combined to quantitatively examine the impact of LULC change on ESV at landscape scale. Such an approach can partly answer the questions of demand-supply pattern and landscape heterogeneity. However, two limiting factors, which may cause biased ESV and further misunderstanding of relationship between landscape pattern and ecosystem services, should be addressed here. The first one is the method for valuing ecosystem services. With the increasing application of land covers extracted from remotely sensed data, pixel-based land covers have been widely used as the proxy for valuing ecosystem services. However, it is noteworthy that such an approach neglected landscape heterogeneity and inevitably mixed the differential pixel-based ESV for the same land cover. For instance, both clear water in headwater and polluted rivers in cities were classified as water, but their roles in ecosystem functioning and individual ESVs were quite different. Therefore, future research focusing ecosystem services assessment should avoid the use of global estimates by simply assigning the coefficients to each pixel of land covers. The second one is uncertainty in understanding of landscape fragmentation. Landscape fragmentation, which partly represents landscape heterogeneity and means spatial and geometric linkage and split of landscape components, plays the key role in determining pattern-process based ecosystem service delivery via changing the transfer of matter, energy, and organism [20]. Herein, landscape heterogeneity and scale-dependent edge effects between varying landscape components should be a noteworthy issue. In our study, the concentric buffers were used as sampling frames for measuring landscape fragmentation and quantifying the relationship between class-level landscape fragmentation and zonal ESVs. If the dense concentric buffers with short interval (e.g., 10–100 m) were drawn, then it may be easily to detect the sharp edge effects of each buffer zone but another problem such as redundant results may arise. Alternatively, a series of spatially uneven buffer zones were drawn with relatively larger distance scale (e.g., 3–6 km). To some extent, the sharp edge effects may be smoothed due to relatively larger distance scale. In Section 3.4 we quantified the relationship between class-level landscape fragmentation and zonal ESVs along the urban-rural gradient. Our regression models can make sense in interpreting the cause-effect relationship between them. However, some sharp edge effects of the boundary of the buffers cannot be completely removed yet. Consequently, the averaged zonal ESVs neglecting anisotropic difference of landscape components may produce puzzling results when other researchers use the same data set if they have different understanding of landscape fragmentation. Therefore, we appeal for future research to focus on development of a standardized procedure so that different case studies using the same methods can be comparable.

5. Conclusions

In this study, Greater Shanghai and Greater Hangzhou were selected as examples to address the troubling issues of rapid urbanization, land scarcity, and accompanying ecological deterioration in eastern China since the late 1970s. It was observed that outward expansion of developed land in these two fast-growing metropolitan regions resulted in substantially declining productive agricultural land, natural land and semi-natural land. This led to remarkable landscape fragmentation and deteriorated regional ecosystem functioning. We highlighted that in both regions, the status of regional ecosystem services degraded remarkably, largely due to unplanned and poorly managed urban sprawl. Embedded in the complex ecological-economic-geographical processes, surging urban expansion and population growth in both regions will inevitably demand more land for development, and thus will exacerbate regional ecosystem services, which feed the development boom of the human-dominated ecosystem. Therefore, on regional and national levels, future policies on land use and urban development must reject any land development motives purely towards achieving economic goals and impairing ecosystem function and services.