Science of the Total Environment 512–513 (2015) 215–226

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Geographic variations of ecosystem service intensity in Fuzhou City, China Xisheng Hu a, Wei Hong b, Rongzu Qiu a, Tao Hong b, Can Chen b, Chengzhen Wu b,c,⁎ a b c

College of Transportation and Civil Engineering, Fujian Agriculture and Forestry University, Fuzhou 350002, China College of Forestry, Fujian Agriculture and Forestry University, Fuzhou 350002, China College of Ecology and Resource Engineering, Wuyi University, Nanping 354300, China

H I G H L I G H T S • • • •

Used ESDA and semivariance to assess geographic variations of ecosystem service Spatial dissimilarities in ecosystem service were observed in Fuzhou City, China. Used the range of semivariance as threshold to develop more precise clusters Land use/cover has great impacts on ecosystem services while varying by locations.

a r t i c l e

i n f o

Article history: Received 8 September 2014 Received in revised form 12 January 2015 Accepted 12 January 2015 Available online xxxx Editor: E. Capri Keywords: Ecosystem services Spatial autocorrelation Semivariance analysis Geographically weighted regression Fuzhou City

a b s t r a c t Ecosystem services are strongly influenced by the landscape configuration of natural and human systems. So they are heterogeneous across landscapes. However lack of the knowledge of spatial variations of ecosystem services constrains the effective management and conservation of ecosystems. We presented a spatially explicit and quantitative assessment of the geographic variations in ecosystem services for the Fuzhou City in 2009 using exploratory spatial data analysis (ESDA) and semivariance analysis. Results confirmed a significant and positive spatial autocorrelation, and revealed several hot-spots and cold-spots for the spatial distribution of ecosystem service intensity (ESI) in the study area. Also the trend surface analysis indicated that the level of ESI tended to be reduced gradually from north to south and from west to east, with a trough in the urban central area, which was quite in accordance with land-use structure. A more precise cluster map was then developed using the range of lag distance, deriving from semivariance analysis, as neighborhood size instead of default value in the software of ESRI ArcGIS 10.0, and geographical clusters where population growth and land-use pressure varied significantly and positively with ESI across the city were also created by geographically weighted regression (GWR). This study has good policy implications applicable to prioritize areas for conservation or construction, and design ecological corridor to improve ecosystem service delivery to benefiting areas. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Ecosystem services, the benefits people obtain from nature resources, are of irreplaceable value for humanity (Costanza et al., 1997; Millennium Ecosystem Assessment (MA), 2005). In the MA, approximately 60% of the ecosystem services are examined to be currently degraded or are being used unsustainably. Such degradation may grow significantly worse during the first half of the 21st century (Vihervaara et al., 2010). Along with the transformation and degradation of ecosystem, an increasing number of studies concerning ecosystem services have been carried out in recent years, because of their ability to reflect and communicate ⁎ Corresponding author at: College of Forestry, Fujian Agriculture and Forestry University, Fuzhou 350002, China. E-mail address: [email protected] (C. Wu).

http://dx.doi.org/10.1016/j.scitotenv.2015.01.035 0048-9697/© 2015 Elsevier B.V. All rights reserved.

human–environmental interactions (Busch et al., 2012). Ecosystem service researches have focused on certain service categories, ecosystem types, and geographical areas, while substantial knowledge gaps remain concerning several aspects, e.g., the spatial variations of ecosystem services (Vihervaara et al., 2010; Syrbe and Walz, 2012). Such spatial correlations of ecosystem services have often been assumed rather than demonstrated (Troy and Wilson, 2006; Turner et al., 2007). This lack of knowledge constrains the effective management and conservation of ecosystem functioning. Fisher et al. (2009) has pointed out that ecosystem services are not homogeneous across landscapes, nor are they static phenomena. Barbier (2012) also gives a good example highlighting the spatial variability of ecosystem services. Understanding both the spatial heterogeneity and homogeneity is indispensible in building efficient policies for governing ecosystem service production and consumption (Barrett

216

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

Fig. 1. Location of study area.

et al., 2011). So, researchers have increasingly sought to provide spatially explicit decision making tools for managing incentive mechanisms in ecosystem service (Nelson et al., 2009). An example provided by Kareiva and Marvier (2007) demonstrates how the mapping of social, biophysical, and ecosystem characteristics can be used to proactively prioritize ecosystem-service-based interventions and management. Despite general agreement on the spatial variability of ecosystem services, it is still not clear which type of spatial relationships (e.g., positive, negative or random) between studied observation units is and what scales do these relationships exist. At present, some methods of spatial statistic including ESDA and semivariance analysis are considered to be extremely effective ways for studying spatial characteristics (Webster and Olivier, 1990; Zawadzki et al., 2009; Zawadzki and Fabijańczyk, 2013). ESDA is a set of Geographical Information System (GIS)-based spatial statistical techniques, which can visualize the spatial agglomeration and anomalies, and reveal the activation mechanism of ecosystem services (McMillen, 2010). Measures of semivariance analysis traditionally have been used for two board purposes: quantification of the scale of variability exhibited by nature pattern of resource distributions and identification of spatial or temporal scale at which a sampled variable exhibits maximum variance (Wallace et al., 2000; He et al., 2007; Zawadzki and Fabijańczyk, 2007). Therefore, the accurate description about the spatial patterns (e.g., positive, negative or random) and gradient changes (scales) of the ecosystem services can be obtained by these two tools. Although ecosystem services not only vary with bio-physical and social properties of the sites, but also with properties of their spatial context (Dalgaard et al., 2007), most studies address their spatial

characteristics within process-related landscape units such as watersheds (Trabucchi et al., 2014), specific habitats (Barbier, 2012), or natural units (Isbell et al., 2011). Nevertheless, cause and effect areas are frequently imprecisely allocated to the processes investigated (Syrbe and Walz, 2012), due to the conflicts resulting from the mismatch between the spatial scale of the process under consideration and the scale at which measurement and observation of census variables are carried out (Anselin, 2001). Identifying the drivers of variability in ecosystem service and quantifying the range at which the variability comes to a steady manner will help scientists and managers to prioritize locations for different goals (Pijanowski et al., 2010; Su et al., 2012). For these reasons, an exploratory research on the spatial relationships between ecosystem services within administrative units and their causes has significant implications for the appropriate policy formulation and implementation. It is also methodologically important because spatial flow planning of services depends on the analysis of spatial correlation at the administrative unit level (Plummer, 2009; Serna-Chavez et al., 2014). Just as the composition and configuration of a landscape's, biophysical features have a direct impact on the type and rate of ecosystem service provisioning, while social heterogeneity is an indirect driver in determining the adoption of conservation incentives for ecosystem services (MA, 2005; Fremier et al., 2013). Spatial complexity of ecosystem services results from the dynamic interaction between the spatial distribution of biophysical cues and variable human actions (Laterra et al., 2012; Fremier et al., 2013). Indirect drivers, such as growing of population concentrations and intensifying of economic activities, can trigger or strengthen direct drivers, such as land-use change and climate

Table 1 Classification system for land-use. Types

Sub-types

Types

Sub-types

1 Cropland

11 Paddy field 12 Irrigable land 13 Dry land 21 Forest 22 Shrub land 23 Orchard, sparse and other woodland 31 Weed 41 River 42 Lake or pond 43 Reservoir

4 Water body

44 Irrigation canals and ditches 45 Wetland or bottomland 51 Urban land 52 Rural settlements 53 Transportation land 54 Other build-up 61 Sand 62 Bare land 63 saline and alkaline land 64 Other

2 Woodland

3 Grassland 4 Water body

5 Construction land

6 Unused land

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

217

Table 2 Standardized equivalent weighting factor of ecosystem services value per hectare for each land types. Land-use type

Gas regulation

Climate regulation

Water supply

Soil formation and retention

Waste treatment

Biodiversity protection

Food

Raw materials

Recreation and culture

Cropland Woodland Grassland Water body Construction land Unused land

14.29 100.00 22.86 0.00 0.00 0.00

32.96 100.00 33.33 17.04 0.00 0.00

2.94 15.70 3.93 100.00 0.00 0.15

37.44 100.00 50.00 0.26 0.00 0.51

9.02 7.21 7.21 100.00 0.00 0.06

21.78 100.00 33.44 76.38 0.00 10.43

100.00 10.00 30.00 10.00 0.00 1.00

3.85 100.00 1.92 0.38 0.00 0.00

0.23 29.49 0.92 100.00 0.00 0.23

change (Weng, 2007). Considering the current trends in population growth and land-use pressures, ecosystem change and vulnerability may be key issues in ecosystem service research (Vihervaara et al., 2010). This study was conducted for the Fuzhou City, one of the fastest growing metropolitan areas in China, however ecosystem services have been seldom taken into account in the policy making process (Hu et al., 2013). Therefore, we employed a benefit transfer approach to assess the provision of ecosystem services (Koschke et al., 2012), then evaluated the ESI of each administrative unit (e.g. town, an administrative and geographical entity which is smaller than county in China) in this area. Then, we employed ESDA and semivariance analysis to detect the geographic variations of such ESI in the study area. The objectives of this research are to (1) examine the characteristic of spatial heterogeneity of ecosystem services; (2) through a regression analysis, we try to identify the drivers of such spatial variation; and (3) provide policy implications applicable to prioritize areas for conservation or construction, and design ecological corridor to improve delivery of ecosystem services to benefiting areas.

2. Materials and methods 2.1. Study area Fuzhou (25°157′N to 26°39′N, 118°08′E to 120°31′E) is the capital and the largest prefecture-level city in the Fujian province, China. It is situated by the East Sea and in the lower reaches of the Minjiang River (Fig. 1.a). It governs 5 districts and 8 counties (Fig. 1.b), consisting of 147 towns and covering a total area of more than 12,153 km2. However, landscapes in this region have changed dramatically and become more fragmented and degraded (Hu et al., 2013). Consequently, a case study of the city can be a representative examination of spatial pattern in ecosystem services, a good example for those of rapidly growing cities in other subtropical parts of the world.

2.2. Data collection The land-use datasets of each town and a vector map of administrative divisions at the town level in 2009 were from the Fuzhou Municipal Bureau of Land and Resources. These datasets were acquired mainly by aerial photograph interpretation aided by ground survey and former land-use maps. They were legal datasets for government planning. The original land-use categories include two levels of national standard types, while they were reclassified into six categories in this study, including cropland, woodland, grassland, water body, construction land, and unused land (Table 1). The demographic and economic data of 2009 at the town level was obtained from statistical yearbooks for Fuzhou City, including variables of population density and per capita industrial output. The land-use and the demographic and economic datasets were then joined into the vector map of administrative divisions using ESRI ArcGIS 10.0 for the subsequent calculation of the values of ecosystem services and spatial analysis. 2.3. Calculation of ecosystem service value Costanza et al. (1997) classified the global biosphere into 16 ecosystem and 17 service function types, and estimated the value of each type of ecosystem services. According to the Costanza et al.'s results, Xie et al. (2003) extracted the equivalent weighting factors of ecosystem services per hectare of terrestrial ecosystems in China. Referring to the equivalent weighting factors, we assigned the values of ecosystem services per unit area of each land-use category in Fuzhou in accordance with the nearest equivalent ecosystems, e.g., woodland equates to forest and unused land equates to wilderness land. We then standardized the values to a relative scale (0–100 value points) (Table 2). These were done by Eq. (1), where Inorm is the standardized equivalent weighting factor of ecosystem service value per hectare for a given land type, whose score is between 0 and 100. I is the assigned value of certain ecosystem service of a given land type. Imin and Imax correspond to the minimum and maximum value of certain ecosystem service of all land types.


 I−Imin  100 I max −I min

Inorm ¼

ð1Þ

In order to obtain an overall performance value for each town, we used a linear additive value function [Eq. (2)] to combine ecosystem services.

Ei ¼

n X

wk sik

ð2Þ

k¼1

wk ¼

m X

wkj =m

ð3Þ

j¼1

Fig. 2. Radar chart displaying the estimated score of each land-use type in terms of ecosystem services according to the standardized equivalent weighting factor for the nine ecosystem services in Table 2.

where Ei represents the overall score (or total utility) of town i; wk is the ecosystem service value of land-use type k (k = 6), which is the sum of the estimated weight of all kinds of ecosystem service wkj divided by

218

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

Fig. 3. Spatial distribution of population density in 2009.

number of ecosystem service categories j (j = 9) (scaled to 0–100 value points) [Eq. (3)], and sik is the area of land-use type k in town i. The weight wk represents the relative importance of certain land-use type for ecosystem service, presented in Fig. 2 according to Table 2. By multiplying the score for each land-use with the corresponding area and adding all services, the overall performance of a town can be calculated in terms of ecosystem services. Considering the spatial variability of ecosystem services (Shi et al., 2012), the ESI was employed to compare ecosystem service function in one town with that of another town and to estimate the spatial heterogeneity at the town level across the study area [Eqs. (4)–(5)].

ESI ¼ Ei =Si

Si ¼

n X

sik

ð4Þ

ð5Þ

k¼1

where Si is the total area of town i.

The global measure of Moran's Ig is as [Eq. (6)]: XX i

  wij ðxi −μ Þ x j −μ

j

1 Ig ¼ N 0 XX X 2 @ wij A ðxi −μ Þ i

j

ð6Þ

i

where wij is the row-standardized contiguity matrix, xi and xj are the ESI at town i and j respectively, and μ is the average level of ESI, and N is the total number of the towns in the study area. Moran's Ig ranges approximately from +1 (for positive spatial autocorrelation) to −1 (negative autocorrelation) and zero expresses the absence of spatial autocorrelation (Anselin, 1995, 2003). The Moran's Ig cannot discover hot-spots and cold-spots between towns. Therefore, the Local Indicator of Spatial Association (LISA) was applied to measure the local spatial association and to indicate the significance of hot-spots and cold-spots (Yang and Wong, 2013). Local Moran's Il statistic [Eq. (7)] was employed to show LISA in this research. Il ¼ X

X

xi −μ 2

ðxi −μ Þ

  wij x j −μ

ð7Þ

j

i

2.4. Geographic variation analysis 2.4.1. Spatial autocorrelation analysis ESDA is employed to detect the existence of spatial autocorrelation for ESI in Fuzhou City. ESDA is a set of techniques used to describe and visualize spatial distributions, to discover patterns of spatial associations, to identify hot-spots and cold-spots and to suggest spatial regimes (Anselin, 1999; Zhang et al., 2011). The analysis focused on two aspects of spatial clustering, namely, the overall “global” spatial clustering and the “local” patterns of ESI distribution.

Both the global Moran's Ig and the local Moran's Il were calculated using ESRI ArcGIS 10.0. When performing the program, the spatial relationships between towns were determined by the method of inverse distance weighting (Łukaszyk, 2004). In the method, neighborhood size determines how many towns are included and a weighted average is taken of the observation values within this neighborhood. In this case, the neighborhood size was set-up automatically using a default value of its radius, which is a Euclidean distance ensuring that each observation has at least one neighborhood. A cluster map, incorporating information about the significance of the local spatial patterns, was then created. In

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

219

Table 3 Descriptive statistics of variables. Variables

Mean

Standard deviation

ESI (score/ha) Indirect drivers

43.780 8.900 5.370 11.3 2.2 56.8

11.717 22.297 9.899 15. 7 2.1 25.9

Direct drivers

Population density (PD) (people/ha) Per capita industrial output (PIO) (1000 yuan/person) Proportion of urban and rural construction land (PCL) (%) Proportion of transportation land (PTL) (%) Proportion of forest land (PFL) (%)

particular, the map resulted in a spatial typology consisting of five categories of towns in terms of ESI: (1) High–High (hot-spots) indicates towns with high level of ESI that are adjacent to towns with high level of ESI (positive spatial autocorrelation); 2) Low–Low (cold-spots) indicates towns with low level of ESI that are adjacent to towns with low level of ESI (positive spatial autocorrelation); (3) Low–High indicates towns with low level of ESI that are adjacent to towns with high level of ESI (negative spatial autocorrelation); (4) High–Low indicates towns with high level of ESI that are adjacent to towns with low level of ESI (negative spatial autocorrelation) and (5) “not significant” indicates towns with no spatial autocorrelation. 2.4.2. Semivariance analysis Semivariance, a central tool of geostatistics was used to measure spatial continuity of neighboring observations (Zawadzki et al., 2005). The value of the experimental semivariance for a vector is derived from calculating one-half the average squared difference between every data pairs separated by a specific lag distance of h (Krige, 1966a; Chen and Feng, 2013). The standard equation for the semivariance is as Eq. (8):

1 X 2 ½zðxi Þ−zðxi þ hÞ 2NðhÞ i¼1 N ðhÞ

γ ðhÞ ¼

ð8Þ

where γ(h) is the experimental semivariance value at distance interval h, describing the degree of autocorrelation present; z(xi) is the measured sample value (i.e. ESI) at town xi; z(xi + h) is the sample value (i.e. ESI) at town xi + h; and N(h) is the total number of sample data pairs separated by distance h. Semivariance functions are usually characterized by the following parameters (Krige, 1966b; Treitz and Howarth, 2000; Zawadzki et al., 2005): Sill (C0 + C) is the maximum level of semivariance, the sum of total variance explained by the spatial structure and Nugget effect; Range (A0) is the distance at where semivariance reaches the maximum, or at which two data pairs are uncorrelated; Nugget effect (C0) is the value of the point where the extrapolated relationship between the two variables intercepts the semivariance axis, represents spatially independent variance; Spatially dependent structural variance (C) is the level of structured variability with the data; and the ratio of Nugget effect (C0) to the Sill (C0 + C) is also an essential predictable measurement of variable (Cambardella et al., 1994), representing the random component proportion in the spatial variance. Linear, spherical, exponential and Gauss models were used to fit the goodness of the experimental semivariance. These parameters of the models were calculated using SAM v 3.1 (Rangel et al., 2010). Semivariance models were cross-validated, coefficient of determination (r2) and the root mean square error of cross-validation (RMSECV) were used to compare the goodness of the models (Davis and Sampson, 2002).

Fig. 4. Spatial distribution of ESI in 2009.

220

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

where, yi is the value of ESI for the ith town, and xi is determinant of the ith town. Similarly, β0 is the constant and βj is the coefficient to be estimated for determinant and ei is the stochastic error term. The GWR model is a modified version of Eq. (9) as shown in Eq. (10): yi ðui ; vi Þ ¼ β0 ðui ; vi Þ þ

Fig. 5. Trend surface analysis of ESI. X and Y are the geographic coordinate axes, Z is the response variable (ESI); the size of the pole is proportional to the value of ESI in each town. The blue curve indicates the variation tendency of ESI in the direction from north to south, while the red curve represents the tendency in the direction from west to east. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.5. Modeling the determinants of ESI 2.5.1. Models Firstly, Ordinary Least Square (OLS) was adopted to examine the factors of ESI. OLS assume that the relationship between dependent and independent variables holds the same everywhere within a given geographical range. Thus it cannot capture the characteristic of spatial non-stationarity for the distribution of ESI. Therefore, the parameter estimations of OLS approach are likely to be biased and inefficient. Then, we employed a varying parameter model (GWR) (Fotheringham et al., 2002) that provides tests for the existence of the spatial interaction between variables across the locations. The OLS regression model was specified as Eq. (9):

yi ¼ β0 þ

X

β j xi þ ei

ð9Þ

X

β j ðui ; vi Þxi þ ei ðui ; vi Þ

ð10Þ

where (ui,vi) represents the spatial location or the geographic coordinates of an observation (i.e., town). In our case, ui and vi are, respectively, the longitude and latitudes of the centroids of the ith town. The GWR fits this model to a group of observations near a given town, thereby yielding a separate parameter estimate for each town. In addition to separately estimating parameters for each observation, GWR puts more weight on observations (i.e., towns) nearer the ith town than those farther away (Fotheringham et al., 2002). The coefficients estimated from the GWR are as Eq. (11):   ^ ðu ; v Þ ¼ X 0 wðu ; v ÞX −1 X 0 wðu ; v ÞY β i i i i i i

ð11Þ

where w(ui,vi) is the diagonal spatial weight matrix that is unique for each observation (i.e., town). The elements of the spatial weight matrix are the essence of the model; these parameter values are based on the spatial location and the bandwidth are used to describe the situation for a spatial location in order to explore the non-stationarity characteristics (Poudyal et al., 2012; Yang and Wong, 2013). In our case, both the OLS and GWR models were estimated using SAM v 3.1 (Rangel et al., 2010). For GWR models, we used Gaussian function to assign spatial weight matrix, and a cross-validation process to set the optimal bandwidth to minimize Akaike Information Criterion (AICc). 2.5.2. Variables Ecosystem services are affected by both direct and indirect drivers (MA, 2005). Population growth and land-use changes are among the key factors (Vihervaara et al., 2010). In this paper, coefficients of density

Fig. 6. The cluster map of spatial association for ESI.

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

221

Table 4 Parameters of spatial heterogeneity of ESI for different models. Model

C0

C

C0 + C

A0

C0 / (C0 + C)

r2

RMSECV

Linear Spherical Exponent Gauss

2.303 14.41 0 56.387

157.217 138.219 154.034 108.207

159.520 152.629 154.034 164.594

48,964.894 45,097.479 56,165.412 53,233.404

0.014 0.094 0.000 0.343

0.376 0.178 0.144 0.209

4.186 4.350 4.259 5.708

of population and per capita industrial output were used to represent population growth; the proportions of three main land types (urban and rural construction land, transportation land and forest land) were used to describe the land-use pressures. Fig. 3 showed the distribution of population density at the town level in 2009. The description of the dependent and independent variables was provided in Table 3. Due to the availability of data, 14 towns in Pingtan County were excluded from the observations, and 133 towns were finally selected as study samples in the regressions.

3.4. Determinants of ESI

Fig. 4 showed that the urban central area had the lowest level of ESI, towns in south-eastern had low or moderate level of ESI, while many towns in extensive western had the highest level of ESI. In general, the trend surface analysis in Fig. 5 indicated that the level of ESI tended to be reduced gradually from north to south and from west to east of the study area.

Table 5 presented the results from the OLS and GWR regressions for ESI. The adjusted R-squared indicated that both of the models fit well in this case, while the values of AICc and Residual showed an improvement of the GWR model over the OLS one. This supported our use of the locally weighted regression model here. Town level regression coefficients estimated by GWR were mapped in Figs. 7 through 9. Figs. 7 and 9 showed that all the towns had significant associations between the proportion of urban and rural construction lands (and the proportion of forest lands) and the ESI, while there were a few portions of the towns in Fig. 8 showing insignificant associations between the proportion of transportation lands and the ESI. Figs. 7 and 8 showed that the constructions of urban, rural and transportation had the reducing tendency to affect the ESI from east to west parts of the study area, while the forest cover exhibited an opposite trend to become more important role in the ESI in the direction (Fig. 9). While, it was unanticipated that the effects of the coefficients of population density and per capita industrial output on ESI were examined to be insignificant for all the towns in the study area.

3.2. Spatial autocorrelation of ESI

4. Discussions

In the analysis of spatial statistics, the default value of neighborhood size was 23,986 m. The value of Moran's Ig statistic was 0.650 and statistically significant at 5% level, indicating a tendency toward the geographical clustering of similar towns with high (or low) level of ESI. Fig. 6 showed the cluster map of ESI. In the map, a large part of towns were noted to have positive spatial autocorrelation (i.e., High–High clusters or Low–Low clusters) significantly at 5% level, while no negative spatial autocorrelation (i.e., Low–High clusters or High–Low clusters) was detected in the study area. The High–High clusters were seated in the west and north sections of Fuzhou, while the Low–Low clusters were mainly situated in the urban central area (see Fig. 1 also).

4.1. Scale problem of geographic variation

3. Results 3.1. Spatial distribution of ESI

3.3. Spatial heterogeneity of ESI Table 4 described the parameters of the fitting models for ESI. Among the four models, the linear model best satisfied the hypothesis according to the value of r2 and RMSECV. Spatial heterogeneity across regions is mainly constituted with the random and autocorrelation parts. According the linear model, C0 / (C + C0) was 0.014, indicating that there were still non-structural (random) factors at some smallscales affecting the level of ESI at the town level, but structural factors were leading sections for spatial heterogeneity in ESI in the study area. Further, the linear model also suggested that the observations were fully independent when the observation distance exceeded 48,965 m.

Pattern is rooted in spatial heterogeneity which in turn stems from variations of spatial dependence (Wu, 2004). Spatial heterogeneity and dependence are ubiquitous across all scales and form the fundamental basis of the ecological phenomena. Scale is the spatial parameterization and temporal parameterization of our perceptive window on reality. The relationship between pattern and scale has been recognized in both geography and ecology for decades (Meentemeyer, 1989; Levin, 1992). The scale problem refers to the results computed from the same data may exhibit various patterns at different scales (e.g., the grain size, lag, or extent of the data sets), or patterns may be best characterized at a distinctive scale (Wu, 1999; Zawadzki et al., 2005). To detect the characteristic scales of landscape patterns, two landscape-scale geographic variation approaches (i.e., autocorrelation and semivariance analyses) are available (Bellehumeur et al., 1997). Anselin (2001) has pointed out that environmental characteristics are often obtained from physical measurements that result from an underlying spatial sampling design. Following this logic, our primary objective was to present an analytical approach enabling the prediction of statistical parameters and features of the spatial variations, which would be observed if a survey for ESI had been designed using different sizes of sampling units (Fig. 10). However, the results showed that all the statistical parameters and features behaved erratically in response to changing sampling sizes. Though, estimation error may result from

Table 5 Parameter estimates from OLS and GWR regressions. Model

Constant

PD (people/ha)

PIO (1000 yuan/person)

PCL (%)

PTL (%)

PFL (%)

Adjusted R-squared

AICc

Residual

OLS GWR

33.299*** 33.625***

−0.007 −0.007

−0.012 −0.013

−0.314*** −0.311***

−0.436*** −0.442***

0.266*** 0.265***

0.975 0.979

552.507 530.799

443.480 355.020

Notes: *** indicates the significance of parameters at the level of 1%.

222

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

Fig. 7. Geographically weighted regression clusters for urban and rural construction lands.

a lack of representativeness of sampling size (Burnett and Blaschke, 2003), these relationships [C0 / (C0 + C)] clearly showed that random fluctuations became less important in our case. Considering one of the purposes of this study was to identify the drivers of the spatial pattern, therefore, we chose the contrast sampling design (an administrative unit at the town level) to match the spatial scale of socio-economic measurements in this research.

4.2. Implication for flows of ecosystem services Fig. 4 suggests that there exist considerable dissimilarities between areas in terms of ESI. This may imply that some ecosystem services are ‘delivered’ from provisioning to benefiting areas through either biophysical or anthropogenic processes (Serna-Chavez et al., 2014). For the sake of the efficient designing for ecosystem services flowing, it is

Fig. 8. Geographically weighted regression clusters for transportation lands.

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

223

Fig. 9. Geographically weighted regression clusters for forest lands.

necessary to investigate the spatial distribution of services (Wunder, 2007; Guariguata and Balvanera, 2009). By identifying spatial heterogeneity and dependence, we are able to describe the distribution of ecosystem services in a more precise and explicit way than traditional inequity indices and choropleth mapping (Yang and Wong, 2013). Hence, this research has some policy implications applicable to prioritize areas for conservation or construction, and design ecological corridor to improve delivery to benefiting areas. The study of spatial variations requires accounting for a threshold distance at which spatial interpolation or processing is valid. According

to the semivariance analysis, the threshold distance (i.e., range) was approximately 50,000 m. Therefore, we processed the spatial autocorrelation analysis again using the threshold distance as neighborhood size instead of the default value (23,986 m). The result in Fig. 11 showed that there were more areas with positive spatial autocorrelation (i.e., High–High clusters or Low–Low clusters) significantly at 5% level than those in Fig. 6, and it was more similar with the distribution of the choropleth mapping (Fig. 4). This finding verifies the previous finding concerning the application of geographic variation analysis (Yang and Wong, 2013), which indicates that more precise information on

Fig. 10. Parameters of spatial heterogeneity respond to sample plot areas.

224

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

Fig. 11. The cluster map of spatial association with a neighborhood size of 50,000 m for ESI.

the distribution of objects can be obtained by such analysis. Fig. 11 revealed existent or latent providing (i.e. High–High clusters) and benefiting (i.e. Low–Low clusters) areas for ecosystem services. These features of the map highlight the possibilities for improving delivery to benefiting areas, e.g., by developing natural resource recreation areas and their connectivity to the benefiting areas. Considering the differences between ecosystem services, studies for specific service are needed to facilitate the analysis of spatial flows of services. 4.3. Proximate forces of ecosystem services Designing appropriate ecosystem-service-based management scheme requires an understanding of the complex interactions between the socio-ecological constraints and ecosystem services (Kremen, 2005). In the study, the intensity of ecosystem services in each town was principally dependent on its land-use structure, i.e., the proportion of urban and rural construction lands, transportation lands and forest lands. This finding is consistent with a previous study, indicating that land-use and land-cover directly affect the ability of biological systems to support human needs by ecosystem services (Vitousek et al., 1997). Also, an unexpected outcome of the study was that population density and per capita industrial output were examined to be insignificant on ESI in the study area. This claim appears to support a conclusion that neither population nor industrialization alone constitutes the sole and major underlying causes of land-cover change worldwide (Lambin et al., 2001). In reality, peoples' responses to economic opportunities, as mediated by institutional factors, drive ecosystem service changes (Li et al., 2010; Helfenstein and Kienast, 2014). Though the role of population density and industrialization on ESI was relatively small when compared to the direct factors of ESI for an political entity (e.g. town), an increasing population means an expansion of populationdependent land-uses (e.g., housing, infrastructure, factories, and shopping centers), and inevitably, this will encroach on ecological lands (e.g., forestlands and wetlands) (Liu et al., 2005; Li et al., 2010).

Additionally, industrialization has also intensified ecological land losses (Deng et al., 2008). Therefore, population density and industrialization could be main drivers of degradation of ecosystem services locally. As the analysis moved beyond the simple regression models and adopted a spatially varying parameter model, it successfully detected local level variation between variables and ESI. The GWR model confirmed that the relationship varies by locations, which demonstrates that the magnitudes of negative/positive relationships do not hold everywhere; however, in the study area where there are overall stable changes in gradient from east to west of the area, this is closely tied to the economy, urbanization and developing processes. We also located some areas in where a variable (e.g. proportion of transportation land) has insignificant effects on ESI. Identified areas where have different effects of variables on ESI, could be targeted by government to implement town level based mitigation initiatives; as these entities (i.e., towns) must contend with both biophysical and social factors at the administrative unit level (Poudyal et al., 2012).

5. Conclusions The analyses performed to establish the geographic variations for the ESI in Fuzhou in 2009 can be generalized as: there exist considerable dissimilarities between towns, with a reducing tendency in the directions from north to south and from west to east as well, and the urban central area was the trough in terms of ESI. An idealized linear semivariance model also showed that for the range of lag distance of less than 50,000 m between observations, these would be dependent. Using the range of lag distance as neighborhood size instead of default value, we can get more precise information on patterns of spatial associations for ESI in the study area. Compared with the OLS model, the results of the GWR method are more reasonable and acceptable for considering the characteristic of spatial non-stationarity for the distribution of ESI in the latter. The

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

analysis results by GWR model showed that the major factor of ESI for each town was its land-use and land-cover. Hence, this research has some policy implications applicable to prioritize areas for conservation or construction, and design ecological corridor to improve delivery to benefiting areas. Considering the differences between ecosystem services, studies for specific service are needed to facilitate the analysis of spatial flows of services. Acknowledgments This research was funded by the National Natural Science Foundation of China (no. 41201100), to which we are very grateful. We are also very grateful for the support provided by the Natural Science Foundation of Fujian Province (no. 2012J01071) and the Science and Technology Major Project of the Hall of Science and Technology of Fujian (no. 2012NZ0001). References Anselin, L., 1995. Local indicators of spatial association: LISA. Geogr. Anal. 27, 93–115. Anselin, L., 1999. Interactive techniques and exploratory spatial data analysis. In: Longley, P.A., Godchild, M.F., Maguire, D.J. (Eds.), Geographical Information Systems. John Wiley & Sons, New York, pp. 253–266. Anselin, L., 2001. Spatial effects in econometric practice in environmental and resource economics. Am. J. Agric. Econ. 705–710. Anselin, L., 2003. GeoDa 0.9 User's Guide. University of Illinois, Spatial, Analysis Laboratory, Urbana-Champaign, IL. Barbier, E.B., 2012. A spatial model of coastal ecosystem services. Ecol. Econ. 78, 70–79. Barrett, C.B., Travis, A.J., Dasgupta, P., 2011. On biodiversity conservation and poverty traps. Proc. Natl. Acad. Sci. 108 (34), 13907–13912. Bellehumeur, C., Legendre, P., Marcotte, D., 1997. Variance and spatial scales in a tropical rain forest: changing the size of sampling units. Plant Ecol. 130 (1), 89–98. Burnett, C., Blaschke, T., 2003. A multi-scale segmentation/object relationship modelling methodology for landscape analysis. Ecol. Model. 168 (3), 233–249. Busch, M., La Notte, A., Laporte, V., Erhard, M., 2012. Potentials of quantitative and qualitative approaches to assessing ecosystem services. Ecol. Indic. 21, 89–103. Cambardella, C.A., Moorman, T.B., Parkin, T.B., Karlen, D.L., Novak, J.M., Turco, R.F., Konopka, A.E., 1994. Field-scale variability of soil properties in central Iowa soils. Soil Sci. Soc. Am. J. 58 (5), 1501–1511. Chen, L.J., Feng, Q., 2013. Geostatistical analysis of temporal and spatial variations in groundwater levels and quality in the Minqin oasis, Northwest China. Environ. Earth Sci. 70 (3), 1367–1378. Costanza, R., d'Arge, R., de Groot, R., Farber, S., Grasso, M., Hannon, B., Limburg, K., Naeem, S., O'Neill, R.V., Paruelo, J., Raskin, R.G., Sutton, P., van den Belt, M., 1997. The value of the world's ecosystem services and natural capital. Nature 387, 253–260. Dalgaard, T., Kjeldsen, C., Hutchings, N., Happe, K., Osuch, A., Damgaard, M., Zander, P., Piorr, A., 2007. Multifunctional farming, multifunctional landscapes and rural development. In: Mander, Ü., Wiggering, H., Helming, K. (Eds.), Multifunctional Land Use. Meeting Future Demands for Landscape Goods and Services. Springer Verlag, Berlin (183–193 pp.). Davis, J.C., Sampson, R.J., 2002. Statistics and Data Analysis in Geology. 3rd edition. John Wiley & Sons, New York (57–61 pp.). Deng, X.Z., Huang, J.K., Rozelle, S., Uchida, E., 2008. Growth, population and industrialization, and urban land expansion of China. J. Urban Econ. 63, 96–115. Fisher, B., Turner, R.K., Morling, P., 2009. Defining and classifying ecosystem services for decision making. Ecol. Econ. 68 (3), 643–653. Fotheringham, A.S., Brunsdon, C., Charlton, M., 2002. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. John Wiley & Sons, New York (27–64 pp.). Fremier, A.K., DeClerck, F.A., Bosque-Pérez, N.A., Carmona, N.E., Hill, R., Joyal, T., Keesecker, L., Klos, P.Z., Martinez-salinas, A., Niemeyer, R., Sanfiorenzo, A., Welsh, K., Wulfhorst, J.D., 2013. Understanding spatiotemporal lags in ecosystem services to improve incentives. Bioscience 63 (6), 472–482. Guariguata, M.R., Balvanera, P., 2009. Tropical forest service flows: improving our understanding of the biophysical dimension of ecosystem services. For. Ecol. Manag. 258 (9), 1825–1829. He, Z.B., Zhao, W.Z., Chang, X.L., 2007. The modifiable areal unit problem of spatial heterogeneity of plant community in the transitional zone between oasis and desert using semivariance analysis. Landsc. Ecol. 22 (1), 95–104. Helfenstein, J., Kienast, F., 2014. Ecosystem service state and trends at the regional to national level: a rapid assessment. Ecol. Indic. 36, 11–18. Hu, X.S., Wu, C.Z., Hong, W., Qiu, R.Z., Qi, X.H., 2013. Impact of land-use change on ecosystem service values and their effects under different intervention scenarios in Fuzhou City, China. Geosci. J. 17 (4), 497–504. Isbell, F., Calcagno, V., Hector, A., Connolly, J., Harpole, W.S., Reich, P.B., Scherer-Lorenzen, M., Schmid, B., Tilman, D., van Ruijven, J., Weigelt, A., Wilsey, B.J., Zavaleta, E.S., Loreau, M., 2011. High plant diversity is needed to maintain ecosystem services. Nature 477 (7363), 199–202. Kareiva, P., Marvier, M., 2007. Conservation for the people. Sci. Am. 297 (4), 50–57.

225

Koschke, L., Fürst, C., Frank, S., Makeschin, F., 2012. A multi-criteria approach for an integrated land-cover-based assessment of ecosystem services provision to support landscape planning. Ecol. Indic. 21, 54–66. Kremen, C., 2005. Managing ecosystem services: what do we need to know about their ecology? Ecol. Lett. 8 (5), 468–479. Krige, D.G., 1966a. A study of gold and uranium distribution patterns in the Klerksdorp gold field. Geoexploration 4 (1), 43–53. Krige, D.G., 1966b. Two-dimensional weighted moving average trend surface for ore evaluation. J. S. Afr. Inst. Min. Metall. 67, 13–79. Lambin, E.F., Turner, B.L., Geist, H.J., Agbola, S.B., Angelsen, A., Bruce, J.W., Coomes, O.T., Dirzo, R., Fischer, G., Folke, C., George, P., Homewood, K., Imbernon, J., Leemans, R., Li, X.B., Moran, E.F., Mortimore, M., Ramakrishnan, P.S., Richards, J.F., Skanes, H., Steffen, W., Stone, G.D., Svedin, U., Veldkamp, T.A., Vogel, C., Xu, J.C., 2001. The causes of land-use and land-cover change: moving beyond the myths. Glob. Environ. Change 11 (4), 261–269. Laterra, P., Orúe, M.E., Booman, G.C., 2012. Spatial complexity and ecosystem services in rural landscapes. Agric. Ecosyst. Environ. 154, 56–67. Levin, S.A., 1992. The problem of pattern and scale in ecology. Ecology 73, 1943–1967. Li, M.S., Mao, L.J., Zhou, C.G., Vogelmann, J.E., Zhu, Z.L., 2010. Comparing forest fragmentation and its drivers in China and the USA with Globcover v2.2. J. Environ. Manag. 91, 2572–2580. Liu, X.H., Wang, J.F., Liu, M.L., Meng, B., 2005. Spatial heterogeneity of the driving forces of cropland change in China. Sci. China D 48 (12), 2231–2240. Łukaszyk, S., 2004. A new concept of probability metric and its applications in approximation of scattered data sets. Comput. Mech. 33 (4), 299–304. McMillen, D.P., 2010. Issues in spatial data analysis. J. Reg. Sci. 50 (1), 119–141. Meentemeyer, V., 1989. Geographical perspectives of space, time, and scale. Landsc. Ecol. 3 (3–4), 163–173. Millennium Ecosystem Assessment (MA), 2005. Ecosystems and Human Well-being: Synthesis. Island Press, Washington, DC (137 pp.). Nelson, E., Mendoza, G., Regetz, J., Polasky, S., Tallis, H., Cameron, R., Chan, K.M., 2009. Modeling multiple ecosystem services, biodiversity conservation, commodity production, and tradeoffs at landscape scales. Front. Ecol. Environ. 7 (1), 4–11. Pijanowski, B.C., Iverson, L.R., Drew, C.A., Bulley, H.N., Rhemtulla, J.M., Wimberly, M.C., Bartsch, A., Peng, J., 2010. Addressing the interplay of poverty and the ecology of landscapes: a Grand Challenge Topic for landscape ecologists? Landsc. Ecol. 25 (1), 5–16. Plummer, M.L., 2009. Assessing benefit transfer for the valuation of ecosystem services. Front. Ecol. Environ. 7 (1), 38–45. Poudyal, N.C., Johnson-Gaither, C., Goodrick, S., Bowker, J.M., Gan, J.B., 2012. Locating spatial variation in the association between wildland fire risk and social vulnerability across six southern states. Environ. Manag. 49 (3), 623–635. Rangel, T.F.L.V.B., Diniz-Filho, J.A.F., Bini, L.M., 2010. SAM: a comprehensive application for Spatial Analysis in Macroecology. Ecography 33, 46–50. Serna-Chavez, H.M., Schulp, C.J.E., van Bodegom, P.M., Bouten, W., Verburg, P.H., Davidson, M.D., 2014. A quantitative framework for assessing spatial flows of ecosystem services. Ecol. Indic. 39, 24–33. Shi, Y., Wang, R.S., Huang, J.L., Yang, W.R., 2012. An analysis of the spatial and temporal changes in Chinese terrestrial ecosystem service functions. Chin. Sci. Bull. 57 (17), 2120–2131. Su, S.L., Xiao, R., Jiang, Z.L., Zhang, Y., 2012. Characterizing landscape pattern and ecosystem service value changes for urbanization impacts at an eco-regional scale. Appl. Geogr. 34, 295–305. Syrbe, R.U., Walz, U., 2012. Spatial indicators for the assessment of ecosystem services: providing, benefiting and connecting areas and landscape metrics. Ecol. Indic. 21, 80–88. Trabucchi, M., O'Farrell, P.J., Notivol, E., Comín, F.A., 2014. Mapping ecological processes and ecosystem services for prioritizing restoration efforts in a semi-arid Mediterranean river basin. Environ. Manag. 53 (6), 1132–1145. Treitz, P., Howarth, P., 2000. High spatial resolution remote sensing data for forest ecosystem classification: an examination of spatial scale. Remote Sens. Environ. 72 (3), 268–289. Troy, A., Wilson, M.A., 2006. Mapping ecosystem services: practical challenges and opportunities in linking GIS and value transfer. Ecol. Econ. 60 (2), 435–449. Turner, W.R., Brandon, K., Brooks, T.M., Costanza, R., Da Fonseca, G.A., Portela, R., 2007. Global conservation of biodiversity and ecosystem services. Bioscience 57 (10), 868–873. Vihervaara, P., Rönkä, M., Walls, M., 2010. Trends in ecosystem service research: early steps and current drivers. Ambio 39 (4), 314–324. Vitousek, P.M., Mooney, H.A., Lubchenco, J., Melillo, J.M., 1997. Human domination of earth's ecosystems. Science 277, 494–499. Wallace, C.S.A., Watts, J.M., Yool, S.R., 2000. Characterizing the spatial structure of vegetation communities in the Mojave Desert using geostatistical techniques. Comput. Geosci. 26, 397–410. Webster, R., Olivier, M.A., 1990. Statistical Methods in Soil and Land Resource Survey. Oxford University Press, New York (316 pp.). Weng, Y., 2007. Spatiotemporal changes of landscape pattern in response to urbanization. Landsc. Urban Plan. 81, 341–353. Wu, J.G., 1999. Hierarchy and scaling: extrapolating information along a scaling ladder. Can. J. Remote Sens. 25 (4), 367–380. Wu, J.G., 2004. Effects of changing scale on landscape pattern analysis: scaling relations. Landsc. Ecol. 19 (2), 125–138. Wunder, S., 2007. The efficiency of payments for environmental services in tropical conservation. Conserv. Biol. 21 (1), 48–58. Xie, G.D., Lu, C.X., Leng, Y.F., Zheng, D., Li, S.C., 2003. Ecological assets valuation of the Tibetan Plateau. J. Nat. Resour. 18 (2), 189–196 (in Chinese with English abstract).

226

X. Hu et al. / Science of the Total Environment 512–513 (2015) 215–226

Yang, Y., Wong, K.K., 2013. Spatial distribution of tourist flows to China's cities. Tour. Geogr. 15 (2), 338–363. Zawadzki, J., Fabijańczyk, P., 2007. Modeling the spatial variability of soils magnetic susceptibility in the area of Upper Silesian Industrial Area. Ecol. Chem. Eng. Chem. Inzynier. Ekol. 14 (S3), 379–391. Zawadzki, J., Fabijańczyk, P., 2013. Geostatistical evaluation of lead and zinc concentration in soils of an old mining area with complex land management. Int. J. Environ. Sci. Technol. 10 (4), 729–742.

Zawadzki, J., Cieszewski, C.J., Zasada, M., Lowe, R.C., 2005. Applying geostatistics for investigations of forest ecosystems using remote sensing imagery. Silva Fenn. 39 (4), 599–618. Zawadzki, J., Magiera, T., Fabijańczyk, P., 2009. Geostatistical evaluation of magnetic indicators of forest soil contamination with heavy metals. Stud. Geophys. Geod. 53 (1), 133–149. Zhang, Y., Xu, J.H., Zhuang, P.J., 2011. The spatial relationship of tourist distribution in Chinese cities. Tour. Geogr. 13 (1), 75–90.