Coupled analysis on land use, landscape pattern and nonpoint source pollution loads in Shitoukoumen Reservoir watershed, China

Sustainable Cities and Society 51 (2019) 101788

Contents lists available at ScienceDirect

Sustainable Cities and Society journal homepage: www.elsevier.com/locate/scs

Coupled analysis on land use, landscape pattern and nonpoint source pollution loads in Shitoukoumen Reservoir watershed, China ⁎

T

⁎

Lei Zhanga,b, Wenxi Lub, , Guanglei Houc, , Huiyan Gaoa, Hongquan Liua, Yining Zhenga a

College of Urban and Rural Construction, Hebei Agricultural University, Baoding, 071001, China College of New Energy and Environment, Jilin University, Changchun, 130021, China c Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun, 130102, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Land use Landscape pattern Nonpoint source pollution Pollution control

Landscape pattern is a major factor affecting nonpoint source pollution (NPS) of watercourse, and understanding the impact of landscape pattern on NPS is critical for sustainable land management. This study undertook a comprehensive analysis, at sub-basin scales, of the interactions among land use, landscape pattern and NPS load in Shitoukoumen Reservoir watershed in the time period between 1980 and 2005. FRAGSTAT was used to quantitatively analyze the temporal and spatial change of landscape metrics, whilst the responses of streamflow, NPS TN and TP load to land uses were simulated by SWAT model. Redundancy analysis (RDA) was used to analyze the interactions between NPS load and landscape pattern dynamics at sub-basin scales. It was found that cropland and forest were the main forms of land use, increased by 38.0 km2 and 8.5 km2 between 1980 and 2005. Thirteen landscape metrics were chosen by factor analysis (FA) method. The magnitudes of changes in streamflow and NPS TN, TP load were found to vary significantly in the sub-basins, decelerated by forest patterns that are unfragmented, have a high value for the largest patch proportion. These analyses provided scientific support for preparing NPS control strategies from the perspective of landscape ecology.

1. Introduction Non-point source (NPS) pollution caused by various agricultural activities is difficult to quantify (Ongley, Zhang, & Yu, 2010). Worldwide, combatting NPS pollution is becoming an increasingly urgent task, with challenges such as unpredictability, complex process mechanisms, and extensive pollution sources, to be tackled by researchers and management authorities (Yang, Huang, Sun, & Zhang, 2017). In addition to regional climate and topography, pesticide and fertilizer application, irrigation and farming methods, and waste management practices, non-point source pollution is significantly affected by landscape types of land-use and macroscopic landscape spatial pattern (Chen, Fu, Xu, & Gong, 2003). Most of the early research obtained data through historical data analysis, typical sample field monitoring and experimental simulation methods to explore the simple relationship between land use types and water quality or non-point source pollution processes. Small-scale field observations tend to have better controllability, and through first-hand monitoring data, more reliable and clear conclusions can be obtained. However, on-the-spot monitoring often requires a lot of manpower and material resources, and its results have

obvious regional characteristics, limiting its application. With the development of GIS and remote sensing technology, relevant scholars tend to use different scale basins as the basic unit, using RS, GIS and mathematical models to explore the relationship between land use and river sediments, nutrient elements (Lee, Hwang, Lee, Hwang, & Sung, 2009; Pauleit, Ennos, & Golding, 2005). The use of mathematical models to simulate the occurrence and transport of pollutants can effectively compensate for the lack of workload and limited objective conditions in field monitoring (Adu & Kumarasamy, 2018; Hao, Zhang, Wang, & Ouyang, 2012). However, in this kind of research, whether it is an empirical model or a mechanism model, more consideration is given to the relationship between land use types, topography and other factors and non-point source pollution, ignoring the impact of spatial pattern of land use, especially the combination and configuration of land use on non-point source pollution. Many studies had reported about the relationship between water quality and the composition of land use types, such as cropland and urban, were related with stream pollutants positively, while forest and grasslands that were less influenced by anthropogenic activities had negative correlations (Teixeira, Teixeira, & Marques, 2014). Chen, Fu,

⁎

Corresponding authors. E-mail addresses: [email protected] (L. Zhang), [email protected] (W. Lu), [email protected] (G. Hou), [email protected] (H. Gao), [email protected] (H. Liu), [email protected] (Y. Zheng). https://doi.org/10.1016/j.scs.2019.101788 Received 24 January 2019; Received in revised form 23 July 2019; Accepted 15 August 2019 Available online 16 August 2019 2210-6707/ © 2019 Elsevier Ltd. All rights reserved.

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

temporal variations of land use types, landscape pattern metrics and non-point source pollution load based on the Fragstats software, SWAT model, (2) reveal the quantitative relationship between landscape pattern metrics and NPS pollution load at sub-basin scale based on RDA (Redundancy analysis) method.

and Zhao (2006) proposed the theoretical framework of "source" and "sink", which indicated that some landscapes that can inhibit and absorb non-point source pollutants are categorized as sink landscapes, whereas source landscapes refer to other landscapes that contribute substantially to pollutants, and discovered a location-weighted landscape index (LWLI), which takes the landscape types, area, spatial location, and terrain feature into account, to quantitatively assess the impact on source and sink landscape patterns on the NPS pollution process. This method provided a basis for landscape configuration optimization to control the migration and transformation process of NPS in heterogeneous landscape. Recently, many studies have been conducted to quantified relationships of “source” and “sink” landscape patterns and non-point source pollution formation process (Xu & Zhou, 2008). Jiang, Chen, and Chen (2013) established the Grid Landscape Contrast Index based on Location-weighted Landscape Contrast Index according to the “source-sink” theory, which reflected the contribution of landscape in a certain space (e.g. 100 m☓100 m grid) to non-point source pollution of water. Zhou and Li (2015) proposed Slope-HRUs landscape index (SHLI), which could reflect the relationship between landscape pattern and soil erosion processes to a certain extent. de Oliveira, Maillard, and de Andrade Pinto (2017) introduced a Land Cover Pollution Index (LCPI) that combined different LULC categories by inferring on their known ability to serve either as a pollution “source” or “sink”. This method focused on the effect of land use type composition on water quality or NPS pollution, and resulted in a suitable landscape pattern (Chen et al., 2006). But it lacks of recognizing the spatial change of landscape configuration, which is also one of the main factors affecting non-point pollution. And the use of landscape metrics is the main quantitative approach to evaluate the characteristics of landscape composition information, spatial structure, spatial configuration, etc. (Leemans, 2013). Shen, Hou, Li, and Aini (2014) integrated multiple stepwise regression analysis and redundancy analysis to explore the quantitative association between landscape metrics, at both the landscape and class levels, and water quality in the highly urbanized Beiyun River Watershed. Zhang, Liu, and Zhou (2018) established the multiple linear regression models to explore the relationship between landscape metrics and water quality indices at different scales. Cheng et al. (2018) used partial least square regression to assess the impact of land use patterns on water quality at different ecofunctional regions in the Songhua River basin. Wu & Lu, 2019 adopted redundancy analysis, nonparametric deviance reduction approach, bootstrap sampling and other statistical methods to reveal the quantitative relationship between landscape pattern metrics and water quality variables. But these metrics were not of ecological significance and most research was in watershed scale, thus an integrated research method was needed to fully explain the contributions of landscape changes to ecological processes (NPS pollution) at diff ;erent spatial scales. Shitoukoumen Reservoir is a large reservoir at the middle reaches of the Yinma River in a tributary of the Songhua River. It is a major resource of drinking water for Changchun City that has over seven million residents. In recent years, the water quality has degraded, as a result of agricultural non-point source pollution and the main pollutants are TP and TN according to the monitoring data from Changchun Environment Protection Bureau (CCEPB) (Zhang, Lu, An, Li, & Yi, 2012). The serious non-point source pollution not only is a serious threat to drinking water and agricultural irrigation water security, but also restricts the sustainable economic and social development of Shitoukoumen Reservoir watershed. Many scholars have found that the relationship between landscape characteristics and NPS pollution processes provides important information for effectively addressing NPS management planning and the optimizing design of regional land use type and landscape patterns to improve the NPS pollution loads in an agricultural watershed. And better understanding of source contributions and control selection is critical for preparing water quality improvement strategies. So the objectives of this study are to (1) analyze the spatial and

2. Materials and methods 2.1. Description of the study area Shitoukoumen Reservoir is located in Jiutai City, Jilin Province (125°45'E, 43°58'N). It is a large reservoir in the middle of the Yinma River, which is the major first-level branch of the Songhuajiang River covering an area of 4944 km2. Shitoukoumen reservoir catchment lies in the North Temperate Zone with a climate of continental seasonal temperate semi-humid monsoon. The Shitoukoumen reservoir catchment lies in the North Temperate Zone, with a climate of continental, seasonal, temperate and semi-humid monsoon. The location is dry and windy in spring, warm and rainy in summer, dry with early frost and rapid cooling in autumn and with long, cold, multi-northwest winds in winter. The annual average temperature is 5.3 °C. The coldest average temperature is −17.2 °C in January, and the hottest average temperature is 23.0 °C in July. The annual rainfall is 369.9–667.9 mm, and it is mainly concentrated from May to September, which accounts for 80% of annual precipitation. The main soils types are dark brown forest soil, meadow soil, lessivage dark brown soil and black soil. The forest dominated by artificial woods accounts for 36.46% of total area in the catchment, including Juglans mandshurica, Fraxinus mandshurica and Phellodendron amurense (Tang, Zhang, Yang, & Gao, 2009; Zhang, Lu, An, Li, Yi et al., 2012). Land use data between 1980 and 2005 were obtained by the Research Center for Remote Sensing and Geographic Information, Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences (Fig. 1). The land use was classified into six types, namely cropland, forest, grassland, water, residential land, and unused land (Fig. 1). In our previous study, we simulated the runoff and sediment load in the condition of different sub-watershed by the SWAT model and analyzed the relationship between the numbers of sub-watershed and the response of the simulated results in Shitoukoumen Reservoir catchment. Results revealed that the reasonable number of sub-watershed was 25. So taking into account the spatial change of land use and the relationship of land use and nonpoint source pollution load, twenty-five sub-basins were divided by watershed delineation method according to our previous study (Lu, Yi, Zhang, & Li, 2010; Zhang, Lu, An, Li, & Gong, 2012) (Fig. 1). 2.2. Methodology The general framework of this study is illustrated in Fig. 2. First, the land use dynamic change and spatio-temporal change of landscape metrics in 1980 and 2005 were calculated in sub-basin scale with ArcGIS tool and FRAGSTATS software. Second, the streamflow, NPS nutrient yield were simulated and estimated in sub-basin by using the SWAT model under 1980 and 2005 land use scenarios. Third, the correlations of the NPS pollution loads with different land use type and landscape metrics were calculated by redundancy analysis (RDA). Finally, based on the RDA results, some suggestions for NPS pollution control and land use pattern planning were made. 2.2.1. land use dynamic change analysis The quantitative change of land use can be represented by land use dynamicity. The single land use dynamicity can be calculated as:

K=

Ub − Ua 1 × × 100% Ua T

(1)

where K denotes the dynamicity of one land use type over the given 2

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Fig. 1. Location and land use maps (1980, 2005) of the Shitoukoumen Reservoir Watershed (SRW).

Fig. 2. The general framework of this study. 3

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

period, Ub is the area (km2) of a specific land use type at the end of study period, and Ua is the area (km2) at the beginning of the study period, T is the interval length, the number of years. The transition matrix reflects the structure and states of transition change of land use types in two stages. ArcGIS 10.2 software was used to analyze the land use area change and the transition matrix analysis between land use types. Firstly, spatial analysis was carried out to describe the patterns of land cover changes over time and to measure the change rates. Based on the information, a land cover type transformation layer (1980–2005) was then generated and analyzed by using the crosstabulation function of the ArcGIS. Then a transition matrix of land use types was produced. Consequently, the changes of the overall land uses in relation with the gains and losses in each category could be compiled (Liu et al., 2009).

metrics were listed in Table 1 (McGarigal et al., 2012). All the metrics were next subjected to factor analyses using the varimax rotation. To identify the most useful subset of metrics, metrics that constituted factors that produced Eigenvalues > 1 were selected. Then combining the individual landscape structure components extracted by FA into groups based on factor pattern similarity. We selected one to three metrics from the same group according to the loading and the ecological meaning of metric.

2.2.3. Non-point source pollution simulation The SWAT model (Arnold, Srinivasan, Muttiah, & Williams, 1998) was used to simulate the stramflow and nonpoint source pollution (NPS) load (TN, TP) in each sub-basin. The detailed establishment of SWAT model, calibration and validation procedure, and simulation results were provided in our previous paper (Zhang, Lu, An, Li, Gong et al., 2012). In our previous study, the study catchment was first divided into multiple sub-watersheds based on a DEM, which were then further subdivided into hydrologic response units (HRUs) based on topography, land use and soil types. After HRUs had been defined, weather data including rainfall, maximum and minimum temperature, relative humidity, solar radiation and wind speed were imported to the SWAT database. The calibration and validation results showed that SWAT model was successfully applied in the inlet of Shitoukoumen Reservoir to simulate streamflow, TN and TP through satisfactory NashSutcliffe coefficient of efficiency (Ens), coefficient of determination (R2)

2.2.2. Landscape metrics selection and calculation More than 100 different landscape metrics have been recently developed and used to characterize landscape patterns. FRAGSTATS (McGarigal, Cushman, & Ene, 2012) and correlation analysis and factor analysis (FA) were used to clarify interrelationships and reduce redundancy of landscape metrics. This procedure involved first calculating pair-wise Spearman’s correlation coefficients among the landscape metrics in order to identify the strongly correlated metrics. We preliminarily chose 40 metrics based on previous landscape studies (Cushman, McGarigal, & Neel, 2008). The name and description of Table 1 Descriptions of selected Landscape metrics. Metrics (acronym)

Descriptions

Number of Patches (NP) Patch Density (PD) Largest Patch Index (LPI) Landscape Shape Index (LSI) Total Edge (TE) Patch Area Distribution (AREA_MN, _AM) Radius of Gyration Distribution (GYRATE_MN, _AM) Shape Index Distribution (SHAPE_MN, _AM) Fractal Index Distribution (FRAC_MN, _AM) Perimeter-Area Ratio Distribution (PARA_MN, _AM) Related Circumscribing Circle Distribution (CIRCLE_MN, _AM) Contiguity Index Distribution (CONTIG_MN, _AM) Euclidean Nearest Neighbor Distance Distribution (ENN_MN, _AM) Perimeter-Area Fractal Dimension (PAFRAC) Contagion (CONTAG) Interspersion & Juxtaposition Index (IJI) Percentage of Like Adjacencies (PLADJ) Aggregation Index (AI)

Numbers of patches of the corresponding patch type (class) Number of patches relative to total landscape area (m2) Percentage of the landscape comprised by the largest patch Reflect the complex shape of the patches that make up the landscape Sum of the lengths (m) of all edge segments Area (m2) of the patch The mean distance (m) between each cell in the patch and the patch centroid Patch perimeter (m) divided by the square root of patch area (m2), 2 times the logarithm of patch perimeter (m) divided by the logarithm of patch area (m2) The ratio of the patch perimeter (m) to area (m2) 1 minus patch area (m2) divided by the area (m2) of the smallest circumscribing circle

Patch Cohesion Index (COHESION) Connectance (CONNECT) Effective Mesh Size (MESH) Landscape Division Index (DIVISION) Splitting Index (SPLIT) Patch Richness(PR) Relative Patch Richness (RPR) Patch Richness Density (PRD) Shannon's Diversity Index (SHDI) Simpson's Diversity Index (SIDI) Modified Simpson's Diversity Index (MSIDI) Shannon's Evenness Index (SHEI) Simpson's Evenness Index (SIEI) Modified Simpson's Evenness Index (MSIEI)

The average contiguity value for the cells in a patch minus 1 Distance (m) to the nearest neighboring patch of the same type, based on shortest edge-to-edge distance The all patches in the landscape are included in the regression of patch area against patch perimeter. Tendency of land use types to be aggregated Measures the distribution of patch adjacencies and isolates the interspersion or intermixing of patch types Percentage of cell adjacencies involving the corresponding patch type that are like adjacencies The number of like adjacencies involving the corresponding class, divided by the maximum possible number of like adjacencies involving the corresponding class The physical connectedness of the corresponding patch type Number of functional joinings between all patches of the corresponding patch type, divided by the total number of possible joinings between all patches of the corresponding patch type Sum of patch area squared, summed across all patches of the corresponding patch type, divided by the total landscape area (m2) 1 minus the sum of patch area (m2) divided by total landscape area (m2), quantity squared, summed across all patches of the corresponding patch type Total landscape area (m2) squared divided by the sum of patch area (m2) squared, summed across all patches of the corresponding patch type Number of different patch types present within the landscape boundary PR divided by the maximum potential number of patch types specified by the user Number of different patch types present within the landscape boundary divided by total landscape area (m2) Minus the sum, across all patch types, of the proportional abundance of each patch type multiplied by that proportion 1 minus the sum, across all patch types, of the proportional abundance of each patch type squared Minus the logarithm of the sum, across all patch types, of the proportional abundance of each patch type squared Measure of patch distribution and abundance Measure of the distribution of area among patch types An even distribution of area among patch types results in maximum evenness

MN: arithmetic mean; AM: weighted average. 4

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

between 1980 and 2005, the land use transformation networks were drawn by Netdraw software as shown in Fig. 4. The direction of the arrow is the transformation direction of the land use types and the relative thickness of arrow is the transformation percentage of different land use types to the land use pointed by arrow. According to the land use transformation networks in Fig. 4, from 1980 to 2005, 48.3% of grassland, 26.8% of residential land and 48.5% of unused land were transitioned to cropland. In sub-basin 1 and 2, the most obvious change of land use type were residential land, which was transitioned to cropland with change rate of 57.0% and 55.6% for sub-basin 1 and 2, respectively. The land use types of sub-basin 3 were unchanged. In subbasin 4, nearly 100% of grassland, 43.2% of forest and 28.8% of residential land were transitioned to cropland from 1980 to 2005. In subbasin 5, the area of water, grassland, unused land and residential land all largely decreased and were transitioned to cropland with change rate of 96.9% (water), 80.5% (grassland), 51.3%(unused land) and 44.7(residential land). 100% of grassland, 60.8% of water and 97.3% of unused land were transitioned to cropland in sub-basin 6, 7 and 8, respectively. In sub-basin 9, all grassland changed into residential land, and 38.5% of unused land changed into cropland. 44.3% of residential land changed into cropland in sub-basin 10 and 49.9% of water changed into forest in sub-basin 11. All cropland transitioned to residential land in sub-basin 12 and all grassland transitioned to forest in sub-basin 14. There was no obvious change in sub-basin 13, 15, 21, 23 and 24. In sub-basin 16, 30.1% of grassland and 42.0% residential land changed into cropland. Large area of grassland were transitioned to cropland in sub-basin 17(100%), 19(100%), 22(76.9%), 25(59.2%). And all unused land changed into cropland in sub-basin 18 and 20.

values. So the spatial variations of nonpoint source pollution load under two land use scenarios could be simulated and analyzed by SWAT model. 2.2.4. Correlation analysis of the landscape pattern with non-point source pollution loads A DCA (Detrended Correspondence Analysis) was first performed to select the most suitable ordination method, as recommended by Lepš and Šmilauer (2003). As the length of the gradient was less than 3 (0.274), it was selected redundancy analysis (RDA), an ordination method based on the linear model. The analysis was undertaken using CANOCO 4.5 (ter Braak & Šmilauer, 2002). Correlation triplot were drawed to evaluate the correlation of landscape metrics at landscape and class level with non-point source pollution loads. 3. Results and discussions 3.1. Analysis of land use dynamic changes 3.1.1. Change of land use dynamicity Land use types in different sub-basions from 1980 to 2005 were showed in Fig. 3. In 1980, the area of cropland and residential land in the whole watershed accounted for 53.04% and 5.70%, respectively, and the area of forest accounted for 37.55% of the total basin area. So cropland and forest were the main land use types in Sitoukoumen Reservior watershed. The area percentages of cropland in sub-basins located that in the middle and lower reaches are all greater than 50%. Forest is mainly distributed in the upstream and accounts for up to 70% in sub-basin 10, 15 and 17. In 2005, the spatial distribution of land use types was similar with that in 1980. Between 1980 and 2005, the area of residential land decreased by 58.5 km2, while cropland, grassland and forest increased by 38.0 km2, and 12.1 km2 and 8.5 km2, respectively. The increase of cropland in sub-basin 4, 5, 6, 7, 8, 11, 12, 13 and 20 contributed to the increase of the whole region (WR). The land use dynamicity results showed that the grassland showed dramastic change in sub-basins and WR. The dynamicity of grassland was 5.2%, 306.1%, 17.4%, 8.1% and 9.1% in WR, sub-basin 14, 9, 16 and 18, respectively.

3.2. Landscape pattern analysis 3.2.1. Correlation between landscape metrics The change rate of landscape metrics values for two land use scenarios (1980 and 200) in each sub-basin was used for correlation analysis and factor analysis (FA), as shown in Table 2. According to factor analysis results, 7 components had Eigenvalues > 1 after varimax rotation and these components cumulatively explained 91.08% of the variation in the landscape pattern metrics. The first factor generated strong loading for NP (0.902), PD (0.902), TE (0.897), LSI (0.909), AREA_MN (−0.910), PARA_AM (0.905), CONTIG_AM (−0.841),

3.1.2. Land use transfer analysis Based on transition matrix of land use types calculated by ArcGIS

Fig. 3. Land use type change in different sub-basions. 5

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Fig. 4. Land use type transformation network in different sub-basins. F: Forest; W: Water; R: Residential land; G: Grassland; C: Cropland; U: Unused land.

which describes the isolation/proximity of patches.

PLADJ (−0.839), AI (−0.839), and moderate loading for SHAPE_AM (0.698), FRAC_AM (0.776), ENN_MN (−0.620), which describes the area/density/edge (NP, PD, TE, LSI, AREA_MN), shape (PARA_AM, FRAC_AM, CONTIG_AM, SHAPE_AM), isolation and contagion (ENN_MN, PLADJ, AI) of patches. The second factor was strongly associated with LPI (0.954), AREA_AM (0.949), GYRATE_AM (0.950), COHESION (0.926), DIVISION (−0.821), MESH (0.949), SPLIT (−0.924), which describes the area/perimeter (LPI, AREA_AM, GYRATE_AM), connectivity (COHESION) and interspersion (DIVISION, MESH, SPLIT) of patches. The third component was significantly related to PR (0.906), PRD (0.905), RPR (−0.906) and SHEI (−0.839) and moderately related to CONTAG (0.788), IJI (−0.550) and MSIEI (−0.735), which describes the diversity (PR, PRD, RPR, SHEI, MSIEI) and contagion (CONTAG, IJI) of patches. The fourth component had highest positive loading on SIDI (0.810), and moderate loading on CIRCLE_AM (0.669), SHDI (0.736), MSIDI (0.793), SIEI (0.746), which describes the diversity (SIDI, SHDI, MSIDI, SIEI) and shape of patches. The fifth component had strong positive loading on SHAPE_MN (0.910), and moderate positive loading on GYRATE_MN (0.728) and FRAC_MN (0.754), which describes the area/density/edge of patches. The sixth component was characterized by strong loading on PARA_MN (-0.872), CIRCLE_MN (0.807), and moderate loading on CONTIG_MN (0.773), PAFRAC (0.736), which describes the shape of patches. The seventh component only had the highest positive loading on ENN_AM (0.781),

3.2.2. Changes of landscape pattern metrics in sub-basins Based on FA analysis results, AREA_MN, PARA_AM, PLADJ, LPI, COHESION, MESH, PR, CONTAG, SIDI, SHAPE_MN, PARA_MN, ENN_AM, CONNECT were chosen as representative metrics for the subsequent analysis. The spatial and temporal change rates of these metrics in landscape and class level were shown in Fig. 5. The most obvious change of landscape pattern metrics was AREA_MN, averagely reducing 49.9% in landscape-level, but averagely increasing 21.4% for cropland and 88.7% for forest in class-level between two land use scenarios. The next were PARA_AM, PARA_MN, SHAPE_MN, of which the change rate were respectively 24.4%, 17.1%, and -14.0% in landscape-level. And the change rates of ENN_AM, PARA_AM, CONNECT were respectively 14.7%, −18.0%, 19.0% for cropland. The SHAPE_MN, PARA_MN, ENN_AM and MESH changed obvious with 17.8%, −16.2%, 13.2% and 14.3% for forest in classlevel. The results illustrated that the area, shape and contagion of patches changed significantly in sub-basins from 1980 to 2005. 3.3. Nonpoint source pollution simulation The watershed streamflow, NPS total N, P yields were estimated by SWAT for two land use scenarios. And the annual change of streamflow, 6

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Table 2 Rotated component matrix. Metrics

Component 1

NP PD LPI TE LSI AREA_MN AREA_AM GYRATE_MN GYRATE_AM SHAPE_MN SHAPE_AM FRAC_MN FRAC_AM PARA_MN PARA_AM CIRCLE_MN CIRCLE_AM CONTIG_MN CONTIG_AM PAFRAC ENN_MN ENN_AM CONTAG PLADJ IJI CONNECT COHESION DIVISION MESH SPLIT PR PRD RPR SHDI SIDI MSIDI SHEI SIEI MSIEI AI Eigenvalues % of total variance Cumulative %

**

0.902 0.902** 0.091 0.897** 0.909** −0.910** 0.026 −0.547 −0.09 −0.023 0.698** −0.105 0.776** 0.141 0.905** 0.086 0.182 −0.230 −0.841** 0.422 −0.620 −0.268 −0.432 −0.839** −0.216 0.297 −0.068 0.246 0.026 0.142 0.304 0.305 0.304 0.522 0.509 0.534 0.155 0.375 0.247 −0.839** 10.914 27.285 27.285

2

3

4

5

6

7

−0.068 −0.068 0.954** −0.162 −0.138 0.035 0.949** −0.219 0.950** −0.223 0.628 −0.309 0.500 0.082 −0.149 −0.156 0.159 −0.327 −0.051 −0.152 0.232 0.022 0.128 −0.068 −0.100 −0.001 0.926** −0.821** 0.949** −0.924** −0.132 −0.13 −0.132 −0.299 −0.227 −0.24 −0.161 −0.187 −0.135 −0.076 7.521 18.803 46.088

0.015 0.015 −0.021 0.138 0.108 −0.006 −0.003 0.039 0.014 0.099 0.113 0.048 0.122 0.100 0.121 −0.068 0.149 0.098 −0.022 −0.116 0.036 0.110 0.788** −0.018 −0.550** 0.153 0.045 0.073 −0.003 0.058 0.906** 0.905** 0.906** 0.124 −0.011 −0.01 −0.839** −0.47 −0.735** −0.010 5.063 12.658 58.746

0.060 0.060 −0.154 0.281 0.302 −0.041 −0.176 0.156 0.049 0.004 0.238 −0.166 0.217 −0.034 0.294 0.019 0.669** 0.107 −0.346 0.083 0.205 −0.028 −0.347 −0.333 0.053 0.130 0.004 0.299 −0.176 −0.028 0.205 0.207 0.205 0.736** 0.810** 0.793** 0.469 0.746** 0.576 −0.325 4.649 11.622 70.368

−0.359 −0.359 0.004 0.169 0.153 0.362 −0.035 0.728** −0.14 0.910** 0.045 0.754** 0.045 −0.252 0.160 0.030 −0.150 0.117 −0.195 0.159 0.217 0.201 0.004 −0.204 −0.596 0.638 −0.144 0.210 −0.035 0.117 0.107 0.108 0.107 0.025 0.083 0.070 −0.093 0.04 −0.023 −0.204 3.653 9.132 79.500

−0.115 −0.115 −0.060 0.092 0.104 0.119 −0.155 0.190 −0.091 0.153 −0.022 0.220 0.087 −0.872** 0.100 0.807** 0.26 0.773** −0.054 0.736** −0.019 −0.207 −0.184 −0.043 0.133 0.089 −0.152 0.004 −0.155 0.145 0.001 0.001 0.001 0.056 0.004 0.036 0.064 0.011 0.036 −0.04 3.027 7.567 87.067

0.107 0.107 0.020 0.003 −0.052 −0.044 −0.01 −0.096 0.091 0.048 −0.150 0.188 −0.04 0.252 −0.029 0.007 0.359 0.101 0.199 −0.069 0.362 0.781** −0.111 0.201 0.402 0.116 −0.221 −0.200 −0.01 −0.073 0.073 0.075 0.073 0.187 −0.108 −0.099 0.063 −0.175 −0.185 0.211 1.604 4.010 91.077

**p < 0.01.

3.4. Correlation of land use and landscape metrics with non-point source pollution loads

TN and TP loads were calculated based on sub-basin level (Fig. 6). The spatial and temporal change rate of streamflow showed a remarkable decline in most of sub-basins, especially in sub-basin 1, 17, 19, 24, averagely decreasing by 9.3%, 2.0%, 4.3% and 3.1%, respectively from 1980 to 2005 for two land use scenarios. While it also increased by 1.1%, 3.0%, 2.3%, 1.8%, 1.2%, respectively in sub-basin 11, 12, 13, 20, 22. From Fig. 3, in sub-basin 1, 17, 19, 24, the area of cropland reduced and forest increased, while the area of cropland and forest showed opposite change trends in sub-basin 11, 12, 13, 20 under the two land use scenarios, which indicated that the increase of forest and decline of cropland contributed to the decrease of streamflow. The spatial distributions of TN and TP loadings were very similar in sub-basins for the two land use scenarios, as shown in Fig. 6. TN and TP losses were obvious in sub-basin 1, 17, 19, 24 of the watershed, respectively with an average decrease of 24.7%, 6.4%, 10.0%, 7.3% for TN and 19.7%, 4.4%, 7.5%, 3.9% for TP. Similar with streamflow, TN and TP loadings also showed significant increase in sub-basin 11, 12, 13, 20, 22, and the average change rate were respectively 6.8%, 14.7%, 12.1%, 8.3%, 12.2% for TN and 6.1%, 13.3%, 10.2%, 8.2%, 11.1% for TP. The change of TN and TP loads was mostly attributed to intensive agriculture and degradation of forest.

3.4.1. Relationship of land use change, NPS loads and sites The change rate of NPS TN and TP load, area change of land use types in each sub-basin between two land use scenarios were used to RDA analysis. Table 3 showed the results obtained for the first two axes, and the canonical coefficients of the environmental variables (area change of different land use types) that define these two axes. The first two axes retained 66.9% of the total variance. The F-ratio of the first axis was 38.38 with a p-value of 0.002, and the first axis had significant correlation (0.82) with NPS load. The most significant environmental variables in the first axis are the area change of forest (F) and cropland (C), with a correlation of 0.770, −0.751, respectively. Fig. 7 presented RDA ordination diagram (triplot) showing area change of different land use type (red solid arrows), change of NPS loads (blue solid arrow) and sub-basin sites (black hollow circles). The angles among arrows denote the degree of correlation between individual variables, and the smaller the angle, the larger the correlation. In addition, arrows pointing in the same direction indicates positive correlation, on the contrary negative correlation. The length of the arrow represented the contribution of each environmental indexes 7

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Fig. 5. Change rate of landscape pattern metrics in landscape- and class-level.

Fig. 6. Spatial variations in streamflow, TN load and TP load from 1980 to 2005 under two land use scenarios for each sub-basin.

8

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Table 3 The eigenvalues and correlations between landscape, NPS load and ordination axes of RDA. Environmental variables

Axis 1

Axis 2

Environmental variables

Landscape level

Class level Crop

Grassland cropland Residential land Forest Water Unused Grassland – – – – – – Eigenvalues Cumulative percentage variance of NPS data Of NPS-landscape data

−0.015 −0.751** −0.093 0.770** 0.094 0.109 −0.015 – – – – – – 0.669 66.9 99.9

0.054 0.064 −0.140 −0.033 −0.012 −0.054 0.054 – – – – – – 0 66.9 100

LPI AREA_MN SHAPE_MN PARA_MN PARA_AM ENN_AM PLADJ CONNECT COHESION MESH CONTAG PR SIDI Eigenvalues Cumulative percentage variance of NPS data Of NPS-landscape data

Forest

Axis 1

Axis 2

Axis 1

Axis 2

Axis 1

Axis 2

0.038 0.050 0.180 −0.193 −0.083 0.197 0.175 −0.100 −0.230 −0.026 0.129 0.049 −0.362 0.654 65.4 99.5

0.054 −0.143 −0.339 0.168 −0.053 0.182 0.093 −0.187 0.017 0.022 −0.126 −0.055 0.066 0.003 65.7 99.9

−0.115 −0.329 −0.067 0.281 0.351 −0.207 −0.322 −0.097 −0.076 −0.233 – – – 0.688 68.8 99.2

−0.322 0.228 0.540* −0.414 0.312 0.643* −0.252 0.444 −0.286 −0.345 – – – 0.005 69.3 99.9

0.719** −0.063 −0.390 0.237 −0.472 −0.274 0.167 0.063 0.117 0.797** – – – 0.751 75.1 99.5

0.129 −0.183 −0.231 −0.132 −0.026 0.012 −0.079 0.001 −0.111 0.121 – – – 0.003 75.4 99.9

*p < 0.05. **p < 0.01.

has a positive correlation with TN and TP load, indicating that cropland is the main source of NPS pollution. Cropland can be considered as the source of NPS pollution. The reason may be due to increasing inputs of inorganic fertilizer, organic manure, and pesticides that contribute to increased regional crop production. Residual nitrogen and phosphorus nutrient elements are easily lost with rainfall runoff. In addition, conventional agricultural activity (i.e., tillage or excessive fertilization) increases the risk of soil erosion and excessive nutrient applications can result in NPS pollution to stream water following surface runoff (Shi et al., 2017).

(land use or landscape metrics) to the variance among the NPS variables. Obviously, the triplot showed that the increase of streamflow (Q), NPS TN and TP load (N, P) had the closest positive association with area increase of cropland, while negatively related to area increase of forest, implying that the change of NPS pollution load depend on the area of forest and cropland. For the spatial distributions of NPS load, sub-basin 12, 13, 20, 18, 5, 7, 11, 8 were positively correlated with cropland but negatively with forest. These sub-basins all showed a remarkable area increase of cropland and decrease of forest (Fig. 3). The sub-basin 1, 19, 24, 9 were positively correlated with forest but negatively related with cropland. These sub-basins all showed a remarkable increase of forest but decrease of cropland (Fig. 3). Land use change affects the ecological processes such as hydrological cycle, soil erosion, nutrient transformation in the basin, and changes the amount of pollutants entering rivers, lakes and other water bodies, and finally causes changes of the ecological environment quality (Lee et al., 2009). So some studies have dealt with the influences of land use change on streamflow, sediment, and nutrient budgets, as well as water quality (Bin, Xu, Xu, Lian, & Ma, 2018; Fu, Zhao, Chen, Liu, & Lv, 2005). Non-point source pollution takes the water circulation process of the basin as the carrier and driving force. As the relationship between rainfall and runoff changes, the generation and migration process of pollutants will also change (Wu, Long, Liu, & Guo, 2012). So different types of land use have a significant impact on the process of rainfall-runoff, the generation of non-point source pollution, and the interception of pollutants. In this study, the change percentage of forest area is negatively correlated with NPS TN and TP load, which is consistent with most similar researches (Lee et al., 2009; Shi, Zhang, Li, Li, & Xu, 2017). Forest has a strong interception effect on the runoff and sediment, which are the carrier for pollutants transport, so forest land reflects a positive role in controlling water quality by filtration, absorption and translation of pollutants (Lee et al., 2009). Forest land is identified as a sink for NPS pollution. The change rate of cropland area

3.4.2. Relationship of landscape metrics, NPS loads and sites Table 3 showed the results obtained for the first two axes, and the canonical coefficients of the environmental variables (change of landscape metrics) that define these two axes. The first two axes retained 65.7% of the total variance. The F-ratio of the first axis was 20.8 with a p-value of 0.21. All the landscape metrics weakly correlated with the axes. The RDA triplot of landscape metrics, NPS load, sub-basin location was shown in Fig. 8. TN, TP load were closely positive related with SIDI, COHESION, PARA_MN, which indicated that the NPS load were influenced by diversity, connectivity and shape of patches. The streamflow was positively correlated with CONNECT, which indicated that the connectivity of patches has important effect on streamflow. According to the relationship of sites and landscape metrics, sites could be grouped into 3 groups. Sub-basin 11, 12, 13, 23, 20, 5, 7, 15, 20, 2, 17,14 were positively related to SIDI, COHESION, PARA_MN. Sub-basin 18, 22, 4 were positively related to Connect. Sub-basin 1, 19, 21, 25, 6, 9, 24, 16, 8 were positively related Shape_MN, ENN_AM, PLADJ, AREA_MN, CONTAG. 3.4.3. Relationship of class-level metrics, NPS loads and sites According to the above analysis, the change of NPS load is closely Fig. 7. Triplot of land use type, NPS load and sites by RDA. The arrows represented with red lines are the environmental variables: Cropland (C), Forest (F), Grassland (G), Residential land (R), Water (W), Unused land (U).

9

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Fig. 8. Triplot of landscape pattern metric, NPS load and sites by the redundancy analysis. P_M: PARA_MN; P_A: PARA_AM; COH: COHESION; CON: CONNECT; A_M: AREA_MN; S_M: SHAPE_MN; E_A: ENN_AM; PL: PLADJ.

of cropland. The streamflow was positively and closely correlated with COHESION, MESH and LPI, and sub-basin 20, 7, 21 were also related to these metrics, indicating that the increase of streamflow in these subbasins is related to the connectivity and aggregation of cropland patches. The sub-basin 1, 4, 23, 17, 24, 19, 14, 16, 25 were positively correlated with PARA_AM, PARA_MN, while negative with AREA_MN, PLADJ, indicating that the cropland in these sub-basins were changed into smaller and scattered patches from 1980 to 2005, resulting in the decrease of NPS load. For forest landscape metrics, the first two axes retained 75.4% of the total variance. The F-ratio of the first axis was 42.2 with a p-value of 0.008. LPI and MESH was strongly correlated with first axis, with correlation of 0.719 and 0.797. As presented in Fig. 9, NPS load were positively related to PARA_AM, ENN_AM, while negatively and significantly related to LPI, MESH, PARA_MN, PLADJ, COHSION. Subbasin 12, 13, 20, 16, 22, 18 were positively correlated with PARA_AM, ENN_AM, and sub-basin 1, 24, 9, 19 were positively correlated with LPI, MESH, PARA_MN, PLADJ, COHSION. The results indicated that the

related to the area change of cropland and forest, but it is not clearly about the relationship of landscape pattern and sub-basins. So the classlevel metrics in each sub-basin were calculated by Fragstats 4.2, and RDA was used to analyzed the relationship of metrics, NPS load and sites. The RDA ordination results which depicted the correlation between ordination axes and cropland landscape metrics were listed in Table 3. The first two axes retained 69.3% of the total variance. The F-ratio of the first axis was 30.9 with a p-value of 0.018. The second axis Among the landscape metrics, the ENN_AM and SHAPE_MN had strong correlation of 0.643 and 0.540 with the second axis. As shown in Fig. 9, the TN and TP load was positively correlated with the AREA_MN and PLADJ, MESH, while negatively related to PARA_MN and PARA_AM, indicating that the increase of area, decrease of disperse, lower of shape complexity of cropland patches contributed to the increase of NPS load. The sub-basin 15, 12, 5, 2, 10 were correlated with above-mentioned metrics, so the increase of NPS load in these sub-basins is attributed to the increase of concentrated large area

Fig. 9. Triplot of class-level metric, NPS load and sites by the redundancy analysis. P_M: PARA_MN; P_A: PARA_AM; COH: COHESION; CON: CONNECT; A_M: AREA_MN; S_M: SHAPE_MN; E_A: ENN_AM; PL: PLADJ. 10

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

significantly and negatively correlated to the NPS TN and TP, reflecting that agglomeration distribution of large-scale forests play a significant role in NPS pollutants control.

increase of connectivity and aggregation of forest patches contributed to the decrease of NPS load. The streamflow had positive correlation with AREA_MN, and sub-basin 23, 17 were also positively related to AREA_MN, but the correlation was not significant. The landscape pattern strongly influences the ecological process, especially the combination of obstacles, channels and highly heterogeneous regions largely determines the exchange and flow of energy and matter in the landscape (Urban, O’Neill, & Shugart, 1987). The influence of landscape pattern changes on the occurrence, migration and transformation of non-point source pollutants are mainly reflected in the changes of ecosystem material and energy flow processes caused by changes in the landscape composition and configuration of the watershed (Turner, Gardner, & O’Neill, 2001). Optimizing the landscape pattern is a very economic and effective measure to reduce the generation of non-point source pollutants at the source, and continuously intercept, absorb and promote the conversion of NPS pollutants to nonhazardous forms during the process of pollutants transport. Therefore, exploring the relationship between landscape pattern and NPS pollutants is of great significance for watershed landscape planning NPS pollution control. According to the RDA results of landscape metrics and NPS load, the largest patch index (LPI) and effective mesh size (MESH) of forest land were significantly and negatively correlated to the NPS TN and TP. The LPI provides a measure of the size of the largest patch of a given type as a percentage of the total landscape area, which is the characterizations of dominance. MESH represents an intensive and area-proportionately additive measure, which is the characterizations of fragmentation. When many small patches of land with different uses occupy the watershed, LPI and MESH will decrease. Therefore, our results showed that a watershed with a more fragmented forest landscape pattern would result in more NPS pollution. Hao et al. (2012) found that LPI had a negative influence on NPS pollution, and a watershed with a more fragmented and complex landscape pattern would result in more nutrient loss. Rajaei et al. (2017) also indicated that preferable forest patterns for better water quality appear to be those that are unfragmented, have a high value for the largest patch proportion. The forest has an important regulatory role in the hydro-ecological functions of the basin. Forest can promote the redistribution of rainfall, affect soil water movement, change the conditions of production and flow, and thus play a role in improving water quality (Shi et al., 2017). The more broken the landscape, the more scattered the plaque, the more unfavorable to inhibit the process of soil erosion, runoff, etc., thus causing the loss of pollution elements (Uuemaa, Roosaare, & Mander, 2005). Conversely, the agglomeration distribution of large-scale forests can better intercept nutrient elements (Winston, Hunt, Osmond, Lord, & Woodward, 2011). So reasonable land planning, reducing the complexity of the landscape pattern, increasing the connectivity between the landscapes, can reduce the output risk of NPS TN and TP.

Acknowledgements This study was supported by Hebei Province Natural Science Fund (E2019204334). Hebei Agricultural University Science & Engineering Fund Project (No. LG201818). Hebei Agricultural University Introduces Talent Research Project (No. YJ201813). The National Key Research and Development China (No. 2016YFC0500204). References Adu, J. T., & Kumarasamy, M. V. (2018). Assessing non-point source pollution models: A review. Polish Journal of Environmental Studies, 27(5), 1913–1922. Arnold, J. G., Srinivasan, R., Muttiah, R. S., & Williams, J. R. (1998). Large area hydrologic modeling and assessment part I: Model development. Journal of the American Water Resources Association, 34, 73–89. Bin, L. L., Xu, K., Xu, X. Y., Lian, J. J., & Ma, C. (2018). Development of a landscape indicator to evaluate the effect of landscape pattern on surface runoff in the Haihe River Basin. Journal of Hydrology, 566, 546–557. Chen, L. D., Fu, B. J., Xu, J. Y., & Gong, J. (2003). Location-weighted landscape contrast index: A scale independent approach for landscape pattern evaluation based on “source-sink” ecological processes. Acta Ecologica Sinica, 23(11), 2406–2413 in Chinese. Chen, L. D., Fu, B. J., & Zhao, W. W. (2006). Source-sink landscape theory and its ecological significance. Acta Ecologica Sinica, 26(5), 1444–1449 in Chinese. Cheng, P. X., Meng, F. S., Wang, Y. Y., Zhang, L. S., Yang, Q., & Jiang, M. C. (2018). The impacts of land use patterns on water quality in a Trans-Boundary River Basin in Northeast China based on eco-functional regionalization. International Journal of Environmental Research and Public Health, 15(9), 1872. Cushman, S. A., McGarigal, K., & Neel, M. C. (2008). Parsimony in landscape metrics: Strength, universality, and consistency. Ecological Indicators, 8, 691–703. de Oliveira, L. M., Maillard, P., & de Andrade Pinto, E. J. (2017). Application of a land cover pollution index to model non-point pollution sources in a Brazilian watershed. Catena, 150, 124–132. Fu, B. J., Zhao, W. W., Chen, L. D., Liu, Z. F., & Lv, Y. H. (2005). Eco-hydrological effects of landscape pattern change. Landscape and Ecological Engineering, 1, 25–32. Hao, F. H., Zhang, X., Wang, X., & Ouyang, W. (2012). Assessing the relationship between landscape patterns and nonpoint-source pollution in the Danjiangkou reservoir basin in china. Journal of the American Water Resources Association (JAWRA), 48(6), 1162–1177. Jiang, M., Chen, H., & Chen, Q. (2013). A method to analyze “sourceesink” structure of non-point source pollution based on remote sensing technology. Environmental Pollution, 182, 135–140. Lee, S. W., Hwang, S. J., Lee, S. B., Hwang, H. S., & Sung, H. C. (2009). Landscape ecological approach to the relationships of land use patterns in watersheds to water quality characteristics. Landscape and Urban Planning, 92, 80–89. Leemans, R. (2013). Ecological systems. New York: Springer. Lepš, J., & Šmilauer, P. (2003). Multivariate analysis of ecological data using CANOCO. Cambridge: Cambridge University Press. Liu, D., Wang, Z., Song, K., Zhang, B., Hu, L., Huang, N., ... Jiang, G. (2009). Land use/ cover changes and environmental consequences in Songnen Plain, Northeast China. Chinese Geographical Science, 19(4), 299–305. Lu, W. X., Yi, Y. P., Zhang, L., & Li, D. (2010). Influence of different watershed subdivision numbers on simulation results of SWAT model. Water Resources and Power, 28(10), 23–25 in Chinese. McGarigal, K., Cushman, S., & Ene, E. (2012). FRAGSTATS v4: Spatial pattern analysis program for categorical and continuous maps. Amherst, MA: University of Massachusettes. http://www.umass.edu/landeco/research/fragstats/fragstats.html. Ongley, E. D., Zhang, X. L., & Yu, T. (2010). Current status of agricultural and rural nonpoint source pollution assessment in China. Environmental Pollution, 158, 1159–1168. Pauleit, S., Ennos, R., & Golding, Y. (2005). Modeling the environmental impacts of urban land use and land cover change: a study in Merseyside, UK. Landscape and Urban Planning, 71, 295–310. Rajaei, F., Sari, A. E., Salmanmahiny, A., Delavar, M., Bavani, A. R. M., & Srinivasan, R. (2017). Surface drainage nitrate loading estimate from agriculture fields and its relationship with landscape metrics in Tajan watershed. Paddy and Water Environment, 15(3), 541–552. Shen, Z. Y., Hou, X. S., Li, W., & Aini, G. (2014). Relating landscape characteristics to nonpoint source pollution in a typical urbanized watershed in the municipality of Beijing. Landscape and Urban Planning, 123, 96–107. Shi, P., Zhang, Y., Li, Z. B., Li, P., & Xu, G. C. (2017). Influence of land use and land cover patterns on seasonal water quality at multi-spatial scales. Catena, 151, 182–190. Tang, Y. L., Zhang, G. X., Yang, Y. S., & Gao, Y. Z. (2009). Identifying key environmental factors influencing spatial variation of water quality in Upper Shitoukoumen Reservoir Basin in Jilin Province, China. Chinese Geographical Science, 19(4), 365–374. Teixeira, Z., Teixeira, H., & Marques, J. C. (2014). Systematic processes of land use/land cover change to identify relevant driving forces: Implications on water quality.

4. Conclusions This study analyzed the interactions of land use dynamic change, landscape pattern and nonpoint source pollution load using the FRAGSTATS, SWAT and RDA methods in Shitoukoumen Reservoir watershed (SRW). All the analysis was based on sub-basin scale, to identify the spatial change and relationship of land use, landscape and NPS load. The forest and cropland were the main land use type in SRW. The biggest area change was conversion to cropland, increased by 38.0 km2 from 1980 to 2005. The FRAGSTATS and correlation analysis and factor analysis (FA) were used to clarify interrelationships and reduce redundancy of landscape, and 13 metrics were identified. The change of selected metrics illustrated that the area, shape and contagion of patches changed significantly in sub-basins. The increase of NPS TN, TP load is positively correlated with increase of cropland area, but negatively correlated with increase of forest land area. And the largest patch index (LPI) and effective mesh size (MESH) of forest land were 11

Sustainable Cities and Society 51 (2019) 101788

L. Zhang, et al.

Jialing River Watershed, China. Journal of Hydrology, 475, 26–41. Xu, S., & Zhou, H. (2008). The landscape dynamics of “source” and “sink” and its quantification method. Research of Soil and Water Conservation, 15(6), 64–71 in Chinese. Zhang, L., Lu, W. X., An, Y. L., Li, D., & Yi, Y. P. (2012). Multivariate statistical techniques for the evaluation of water quality of Shitoukoumen reservoir catchment. Environmental Monitoring in China, 28(2), 119–123 in Chinese. Yang, B., Huang, K., Sun, D., & Zhang, Y. (2017). Mapping the scientific research on nonpoint source pollution: A bibliometric analysis. Environmental Science and Pollution Research, 24(5), 4352–4366. Zhang, L., Lu, W. X., An, Y. L., Li, D., & Gong, L. (2012). Response of non-point source pollutant loads to climate change in the Shitoukoumen reservoir catchment. Environmental Monitoring and Assessment, 184, 581–594. Zhang, X., Liu, Y. Q., & Zhou, L. (2018). Correlation analysis between landscape metrics and water quality under multiple scales. International Journal of Environmental Research and Public Health, 15(8), 1606. Zhou, Z. X., & Li, J. (2015). The correlation analysis on the landscape pattern index and hydrological processes in the Yanhe watershed, China. Journal of Hydrology, 524, 417–426.

Science of the Total Environment, 470–471, 1320–1335. ter Braak, C. J. F., & Šmilauer, P. (2002). CANOCO reference manual and CanoDraw for windows. User’s guide: Software for canonical community ordination (v. 4.5). Microcomputer power. Ithaca, New York. Turner, M. G., Gardner, R. H., & O’Neill, R. V. (2001). Landscape ecology in theory and practice. New York: Springer. Urban, D. L., O’Neill, R. V., & Shugart, H. H. J. (1987). Landscape ecology: A hierarchical perspective can help scientists understand spatial patterns. BioScience, 37(2), 119–127. Uuemaa, E., Roosaare, J., & Mander, Ü. (2005). Scale dependence of landscape metrics and their indicatory value for nutrient and organic matter losses from catchments. Ecological Indicator, 5, 350–369. Winston, R. J., Hunt, W. F., III, Osmond, D. L., Lord, W. G., & Woodward, M. D. (2011). Field evaluation of four level spreader–vegetative filter strips to improve urban storm-water quality. Journal of Irrigation and Drainage Engineering, 137(3), 170–182. Wu, J. H., & Lu, J. (2019). Landscape patterns regulate non-point source nutrient pollution in an agricultural watershed. Science of the Total Environment, 669(15), 377–388. Wu, L., Long, T., Liu, X., & Guo, J. (2012). Impacts of climate and land-use changes on the migration of non-point source nitrogen and phosphorus during rainfall-runoff in the

12