Measuring the zonal responses of nitrogen output to landscape pattern in a flatland with river network: a case study in Taihu Lake Basin, China.

Landscape pattern changes induced by rapid urbanization and intensified agricultural activities have exerted great pressure on regional water purification services. Relationship between landscape metrics and nitrogen-related ecosystem services has been a major concern of many scholars and has been widely used for guidance for land use and cover (LULC) management. However, clear zonal differences may exist, especially in highly developed reticular river network area, thus limiting our understanding of nitrogen output (NOP) to landscape pattern in the details. The spatial distribution of regional NOP was obtained based on the InVEST model. The zonal responses of NOP to landscape patter were examined under hydraulic subregions and subbasin scale. The results show that the unit value of average NOP in the Taihu Lake Basin (TLB) was 146.14 (kg/km²), and the total output reached 23677.92 t in 2020. The simulation NOP showed reasonable agreement with verified water quality observations in the lake inlet stations, with an R² of 0.76. In terms of space composition, merely cropland have significant effects on NOP in the whole basin scale, while the explanatory variables include cropland and developed land in Pudong (PD), Puxi (PX), Wuchengxiyu (WC), and Hangjiahu (HJ) regions. In Huxi (HX) and Yangchengdianmao (YC) regions, cropland and forest are the significant impact types, while in (Zhexi) ZX region, cropland, developed land, and forest are significant impact types. In the space configuration, the percentage of landscape (PLAND) or largest patch index (LPI) of cropland showed positive effects about NOP, whether in the whole basin or the hydraulic subregions. Edge density (ED) (−3.48), number of patches (NP) (−3.91), and percentage of like adjacencies (PLAND) (−2.80) of the forest exhibit negative correlations with NOP, in the HX, ZX, and YC region, respectively. It displays diversiform in the response of NOP to the landscape metric of developed land, which speculate that the heterogeneity of developed land can also have a constraint on NOP, in the highly urbanized areas with less forest area. In addition, the total nitrogen output of the TLB needs to be controlled, especially in HJ region which was identified as the sensitive area of pollution sources with the largest NOP and should be paid more attention to. Compared with the administrative management unit, it is more reasonable to control and manage the pollution sources by referring to the hydraulic subregions and subbasin units. Senior managers are required to strengthen communication and cooperation with hydraulic subregions across administrative regions. However, when managing NOP through the landscape modifications, measures should be taken to reduce the aggregation of nitrogen sources and increase the fragmentation of nitrogen sinks. As for high aggregation developed and agricultural land regions, the types of land used should be enriched to help the sustainable development.

Keywords: Landscape pattern metrics; Nitrogen output; InVEST model; Stepwise multiple regression analysis; Taihu Lake Basin

As an important link between natural ecosystems and human well-being, ecosystem services have become an important issue being faced in human sustainable development (Olander et al. [29]). As people's awareness of the importance of ecosystem services has increased, relevant managers have also used ecosystem service methods for the regional sustainable management (Giri et al. [14]). Rapid urbanization has changed the types of ground cover and resulted in various problems, such as environmental pollution and ecological degradation. The emergence of these problems was inseparable from the evolution of landscape patterns (Duarte et al. [10]; Duan et al. [9]).

The landscape pattern metrics are highly concentrated on landscape pattern information and have been widely used to reflect the composition and spatial configuration of the landscape structure (Hu and Zhang [21]; Wu [42]). The influence of regional landscape patterns on its ecological environment has always been a research hotspot (Bozorgi et al. [2]). Previous studies emphasized that the composition and configuration of land use and cover (LULC) types could cause heterogeneity in the landscape and further affect the regional ecosystem services (Gao et al. [12]; Ning et al. [28]).

The research on the relationship between the landscape pattern metrics and ecosystem services was relatively comprehensive. In order to further quantify the spatial discrepancy, researchers often choose various spatial units. Previous studies mainly consider the whole basin (Hao et al. [18]) or focused on the buffer zone of the river (Hou et al. [20]), or the buffer zone of the water quality monitoring station (Xu et al. [46]; Zhang et al. [49]), or the administrative boundary unit (Bai et al. [1]). However, less attention has been paid based on hydraulic subregions, with mesoscale horizon. Hydraulic subregions can reflect the similarities between topography and landforms, hydrometeorology, and natural disasters and reflect the discrepancy of hydrodynamics. At the same time, subbasins with a smaller unit area than hydraulic subregions can further quantitatively depict. At present, the division of subbasin is generally based on the surface runoff overflow model (Chang et al. [3]). However, the subbasin division of plain river network area has certain difficulties, such as flat terrain, dense river networks, and lack of hydrodynamic under natural conditions. Nevertheless, dividing subbasin units is a prerequisite for accurately portraying the hydrological process for watershed research (Giri [13]). The importance of basin zoning is mainly from the aspect of water regime management, pollutant transport, and even socio-economic gathering development.

For the simulation of nitrogen output, multiple watershed process models have also been proposed, such as the soil and water assessment tool (SWAT), the hydrological simulation program-Fortran (HSPF), and the agricultural non-point source (AGNPS). Although the above model focuses on the detailed process mechanism, there are still many parameters to be verified in the process of pollutant migration and transformation (Wang et al. [38]). At the same time, high-precision and large amounts of input data also bring challenges to the application of these models (Wang and Shui [36]). The last and most important thing is that as a plain river network area, the TLB often lacks natural hydrodynamic conditions to form runoff, and there are often backflow phenomena. Therefore, it is not applicable to the above model. As an important tool for studying ecosystem services, the InVEST model provides technical support for quantitatively describing the spatial heterogeneity of ecosystem services. At the same time, it also has the advantages of high science, intuitiveness, and easy operation. It is widely used in the research of ecosystem service assessment and ecological management (Xu et al. [45]).

As a concentrated area of economic development in the Taihu Lake Basin (TLB), it had a pollution load that exceeded the carrying capacity of the environment. Water quality degradation was the main environmental problem, especially due to excessive nitrogen caused by agricultural non-point source pollution (Chen et al. [5]). Paying attention to the spatial output layout of nitrogen is a prerequisite for agricultural non-point source pollution control. In terms of spatial heterogeneity expression of ecosystem services, the InVEST model showed certain advantages (Yu et al. [47]).

In order to explore the relationship between landscape metrics and nitrogen output (NOP) from the perspective of mesoscale, which include hydraulic district and subbasin scale, our manuscript proposed the following research objectives based on previous research content: (1) the spatial distribution of NOP in the TLB, (2) the differences in landscape metrics at different hydraulic subregions, and (3) measuring the zonal responses of NOP to landscape pattern, in the whole watershed and hydraulic regions.

The TLB was located in the lower reaches of the Yangtze River, which covered an area about 36,900 km2(Fig. 1). It consisted of Jiangsu provinces, Zhejiang provinces, and Shanghai municipality, thereinto incorporated large- and medium-sized cities (e.g., Hangzhou, Wuxi, Suzhou, and Changzhou cities). According to the hydraulic subregions, the TLB was divided into Huxi (HX), Wuchengxiyu (WC), Yangchengdianmou (YC), Taihu (TH), Zhexi (ZX), Hangjiahu (HJ), Pudong (PD), and Puxi (PX) regions (Fig. 1). It is an economically developed region with a high urbanization rate in Shanghai municipality and Suzhou City. The TLB has a subtropical monsoon climate. The average annual temperature increased from north to south at 15–17°C. The average annual rainfall was 1200 mm, and the precipitation was mostly concentrated from April to June. The TLB presents a topographical distribution of high elevation in the southwest and low in the east. The western part is hilly and covered with mixed forest of deciduous and evergreen deciduous. The Yili Mountains are mainly distributed with subtropical evergreen broad-leaved forest, which is located in the southwest of TLB, while the southern part of the watershed was a plain area with a dense river network. The rice and economic crops, such as bayberry and citrus, grow along the periphery of the lake. Red and yellow-brown soils are abundant, and the main type of cropland is paddy soil.

Graph: Fig.1 (a) Location of Yangtze River Delta, China. (b) Location of Taihu Lake Basin in Yangtze River Delta. (c) Hydraulic subregion boundary of Taihu Lake Basin and locations of water quality station.

The data used in our research include LULC data, meteorological data, water quality data, DEM data, and vegetation biomass data in 2020. The LULC data came from the 30-m global land cover data called GlobeLand30 (http://www.globallandcover.com), which has been proven with good applicability in part of the TLB (Hu et al. [22]; Pan et al. [30]). The main land use types were divided into cropland, forest, grassland, water, developed land, bare land, wetland, and shrubland (Fig. 2). The meteorological data mainly came from the China Meteorological Data Network (http://data.cma.cn/). DEM data were obtained from the geospatial data cloud (http://www.gscloud.cn/), where the ASTER GDEM 30M resolution digital elevation data was available. The water quality observation data were obtained from the Environmental Monitoring Center of Jiangsu Province. The monitoring sites are located in the inlet rivers of the lake, of which main indicators include ammonia nitrogen (NH4_N) and total nitrogen (TN) in 2020 (Fig. 1).

Graph: Fig. 2 Land use and cover in Taihu Lake Basin

The landscape metrics can reflect the spatial configuration status of LULC types (Milovanovi 2020). The present study selected 9 metrics, which showed correlation with water quality at the catchment scale, referring to previous research (Bai et al. [1]; Duarte et al. [10]). Meanwhile, from the perspective of the algorithm, considered culling indicators have collinearity. The chosen landscape metric can reflect the allocation status of LULC, reflecting their aggregation and spatial heterogeneity. The indicators, such as largest patch metric (LPI), connectance (CONNE) and aggregation (AI), can reflect the spatial aggregation. From the perspective of reflecting the spatial heterogeneity, the main indicators include number of patches (NP), edge density (ED), fractal dimension metric (FRACMN), landscape division index (DIVISI), splitting index (SPLIT), and percentage of like adjacencies (PLADJ). Similarly, the selection of the above indicators mainly needs to consider the focus of the description of the landscape metrics from three perspectives: area, shape, and degree of aggregation. The description of specific metrics is shown in Table 1. The data preprocessing was based on ArcGIS 9.2, and the landscape metrics were calculated through FRAGSTATS 4.2 (https://www.umass.edu/landeco/research/fragstats/fragstats.html).

Table 1 Description of landscape metrics

Note: In Equations (1)–(10), eik is the total length (m) of edge in landscape involving patch type i. aij is an area of patch ij. A is the total landscape area (m2). pij is perimeter (m) of patch ij. ni is several patches in the landscape of the corresponding patch type. cijk is joining between patches j and k of the corresponding patch type i, which was based on a user-specified threshold distance. gii is the number of like adjacencies between pixels of patch type i, which was based on the double-count method. max → gii is the maximum number of like adjacencies between pixels of patch type i. pi is the proportion of landscape comprised of patch type i. gik is the number of adjacencies between pixels of patch types i and k based on the double-count method

The InVEST (Version 3.3.3) (https://naturalcapitalproject.stanford.edu/software/invest) was used for map nitrogen source from watersheds and their transport to the stream in nutrient delivery ratio module. Consistent with the export coefficient literature, load values for each LULC class are derived from empirical measures of nutrient export (Reckhow et al. [32]). The amount of nutrient (e.g., fertilizer and livestock waste) and atmospheric deposition and their amendatory factor could be used to obtain the load parameters. The module calculated the purification ability by measuring pollutant interception through storage effects of regional vegetation and soil on raster pixel scale (Han and Li [17]). The higher the NOT of a unit pixel means that the purification capacity is lower for total nitrogen. The NOP load of pixel i (Xi) could be expressed as below:

11 $$X_i={\mathit{\text{load}}}_i\times {\mathit{RPI}}_i$$

Graph

where Xi is the NOP load of pixel i (kg/year) and loadi consists of surface runoff load and underground runoff load, which were calculated by the parameter of nitrogen load (Load-n), which quantifies the ratio of dissolved nutrients over the total amount of nutrients. RPIi is expressed as the runoff potential of pixel i, which could be obtained using precipitation.

This module is based on basic data such as DEM, water production, land classification data, watershed boundary, NOP load coefficient, and plant filtration efficiency. The plant biophysical table is mainly used to describe the nutrient load of each LULC type, which is defined according to the LULC and InVEST user's guide as shown in Table 2.

Table 2 Biophysical information of land use and cover type in InVEST model.

Note: Load-n is nitrogen load. Eff-n is the nitrogen transport efficiency. Crit-len is the maximum transport distance of nitrogen. Root depth is the root length of vegetation

For observational data, use water quality monitoring data and streamflow data to obtain NOP values. The specific algorithm is to use the annual cumulative NH4_N and TN to multiply the annual streamflow data to obtain the total nitrogen output of the observation site. Thereinto, the average streamflow data comes from the "Annual Report on Water Regime of the Taihu Lake Basin and Southeastern Rivers in 2019." For the simulated output, use the water quality monitoring station to obtain the NOP value of the grid which was output by the model, and use the 1km2 area within the site as the minimum verification unit.

Determination correlation coefficients (R2), root mean square error (RMSE), and average absolute errors (MAE) are used to evaluate the simulation effect. The determination correlation coefficient can describe the degree of linear correlation between the simulated and observed value (Equation (12)), which can reflect the mutual relationship and correlation direction between the two variables. However, it cannot accurately indicate the degree of correlation between the two variables. The RMSE (Equation (13)) and the MAE (Equation (14)) are measures of the difference between the simulated and the observed value. Among them, RMSE demonstrated sensitive to maximum and minimum errors of the data, while MAE could reflect the actual situation of the simulated value error well:

12 $$R^2=1-\frac{{\sum }_i{\left(\widehat{y_i}-y_i\right)}^2}{{\sum }_i{\left(\overline{y_i}-y_i\right)}^2}$$

Graph

13 $$\mathit{\text{RMSE}}=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}$$

Graph

14 $$\mathit{MAE}=\frac{1}{n}{\sum }_{i=1}^n\left|\widehat{y_i}-y_i\right|$$

Graph

where $\widehat{y_i}$ represents the simulated value and yi represents the observed data. $\overline{y_i}$ is mean of observed data. The larger the R2 and smaller the RMSE and MAE values means that the model accuracy is desirable.

The division of subbasin based on hydraulic subregions can further portray the particular information under similar topography, geomorphology, and hydraulic conditions. The subbasin division process of the TLB includes revise original DEM image according to the hydrographic net vector diagram in each hydraulic subregions, making it to get more closer to the real terrain. With reference to the river network density method (Gao et al. [11]), in order to ensure the retention and display of river network information, a standard for the catchment area threshold has been set. The threshold of the accumulation of confluence in HJ, HX, and ZX region are 25000 m3, PD and PX are 15000 m3, and WC and YC are 20000 m3. The average area of subbasin unit is about 269 km2. The result of the division is shown in Fig. 3.

Graph: Fig.3 (a) DEM. (b) Descript hydraulic subregions boundary and hydrographic net. (c) Descript subbasin boundary.

Stepwise regression analysis selects independent variables to establish the optimal regression equation. With the confidence level of 0.05 as the independent variable entry condition, the stepwise regression analysis was used in R software to construct a multiple regression equation of the NOP and landscape pattern at the different watershed scales. Both the dependent (NOP) and independent variables (landscape pattern) are analyzed based on LULC data. Therefore, the response relationship analysis emphasizes the relationship between the transportation process of pollutants and the shape, size, and arrangement of landscape patches. It can quantitatively describe the relationship between the migration of pollutants and the degree of heterogeneity and aggregation of the landscape.

The water quality monitoring data of five inlet rivers in 2020 were used for verification about NOP simulation (Fig. 4). The mean value of the simulated and the observed value is 254.93 kg/km2 and 267.52 kg/km2, respectively, which showed insignificant difference. The simulated NOP have reached a generally accepted performance, with R2, RMSE, and MAE of 0.76, 34.81 kg/km2, and 68.39 kg/km2, respectively. However, there are significant deviations in the extreme value, in particular the low value.

Graph: Fig. 4 Boxplot of simulated and measured NOP

The maximum nitrogen output of the TLB was 410.26 kg/km2, and the average value is 146.14 kg/km2 in 2020 (Fig.5). The NOP varied greatly among different subregions (Table 3). The highest NOP was found in HJ region, with an average value of 212.18 kg/km2. The lowest NOP in the TLB was located in TH region, with an average value of 39.55 kg/km2. The highest output value of the whole region is 455.84 kg/km2 in HJ region. In terms of the total NOP of the subregions, the higher areas are HJ and HX regions, and their total NOP is 6544.55 t and 5213.09 t, respectively (Table 3). In the plain river network area of the TLB, crops and citrus orchards were mainly planted, and the application of chemical fertilizers was the main source of nitrogen (Wang et al. [37]).

Graph: Fig. 5 Nitrogen output in the Taihu Lake Basin

Table 3 Regional statistics of nitrogen output in the Taihu Lake Basin

Main LULC types in the TLB include crops, developed land, water, and forest. In 2020, the total area of crops is 16431.27 km2 (44.01%), the area of developed land is 9890.07 km2 (26.48%), the water area is 5724.95 km2 (15.33%), and the area of forest is 4542.61km2 (12.16%) (Fig. 6).

Graph: Fig. 6 The main land use and cover type composition of Taihu Lake Basin

In terms of LULC composition in subregions, HX and HJ regions occupy a larger area than other subregions, which were worth noting that the areas of cropland have reached more than 3000km2. For developed land, the high areas are in the HJ and WC regions, which are about 2006.59 km2 and 1699.68 km2, respectively. Regarding for water, the TH region far exceeds the other regions, with an area about 2326.82km2. Combined with the LULC composition of TH region, water is the main LULC type (97.21%). The water areas of the HJ, YC, and HX regions are 980.83 km2, 896.48 km2, and 887.47km2, respectively. The largest forest area is in ZX region, with about 3274.66km2, followed by the HX region, with about 766.70km2. The diverse types of regional LULC are affected by the various human production activities and the regional economic development.

Landscape metric reflects the regional numerical aggregation degree and spatial heterogeneity. It exhibits the highest aggregation characteristics, showing the largest LPI (95.55), CONNE (1.32), AI (99.31), and PLADJ (99.22) in the TH region (Table 4). The second level of region aggregation is PX region, which reflected LPI (60.01), CONNE (0.56), and AI (96.12). Meanwhile, the PX region has the highest FRACMN (1.09), which also reflects that the regular shape of the regional landscape may be due to artificial developed land.

Table 4 Taihu Lake Basin landscape metric

In the HX and HJ regions have the regional agglomeration levels in the third and fourth, respectively. Among them, LPI (55.14), CONNE (0.13) and AI (95.47) in HX region, and LPI (48.27), CONNE (0.16), and AI (94.62) in HJ region.

In terms of spatial heterogeneity, it shows the largest NP (49660) and ED (53.70), where it covers a wide area of forest (57.39%) in ZX region. It shows strong landscape heterogeneity, e.g., DIVISI (0.95) and SPLIT (19.46) in the YC region, which have the area of developed land, cropland, and water analogously accounting for 36.98%, 42.48%, and 20.09%, respectively.

NOP has the highest sensitivity to the landscape pattern at the watershed scale, with the best fit (R2 = 0.75) using PLAND and LPI of cropland as explanatory variable. PLAND and LPI reveal the area of landscape type information, reflecting the degree of aggregation, which mean the higher its value, the lower diversity of the landscape. The regression equation is as follows: NOP = 3.11 + 2.99 cropland_PLAND +3.02 cropland_LPI, p < 0.05.

The landscape metrics also have the greatest impact on NOP in the subregion (Table 5). In whole of the subregions, both PLAND and LPI of cropland showed a positive correlation with NOP, which is consistent with the whole watershed scale.

Table 5 Stepwise multiple regression analysis of landscape metric and nitrogen output

Note: CroP means cropland, Fore means forest land, and Dev means developed land. ** represents a confidence level of p at 0.01, and * represents a confidence level of p at 0.05

At the subregional scale, the landscape metrics of forest with NOP show a negative relationship. The ED (−3.48), NP (−3.91), and PLAND (−2.80) of the forest all exhibit negative correlations with NOP, in the HX, ZX, and YC regions, respectively. The ED and NP can reflect the heterogeneity of the landscape patch. Measures such as increasing tree planting around developed and cropland could be taken to improve the water purification service.

The response performance of the landscape metric of developed land to NOP is diverse. The PLAND metric showed a positive effect in WC (+2.37), PD (+ 4.46), and HJ (+2.54) regions, and the DIVISI (+3.23) metric showed a positive effect in the ZX region. Zhang et al. ([48]) pointed out that the DIVISI of developed land is the main driving factor of ammonia nitrogen, which has a positive effect. The constraint effect of landscape metric on NOP is PLADJ which showed negative effects in WC (−3.78) and HJ (−6.13) regions. Meanwhile, LPI (−2.64) in WC region and SPLIT (−3.21) in PD region and CONNE (−3.54) in PX region showed negative effects. Both DIVIS and SPLIT reflect the degree of dispersion of patch but show diversiform influence. DIVIS metric shows a positive effect in the ZX region, while SPLIT metric shows a negative effect in the PD region.

The constituent elements behave discrepantly in different regions. The response factors to NOP include cropland, developed land, and forest in ZX region, while it includes cropland and developed land in PD, HX, WC, PX and HJ regions, and the response land types are cropland and forest in HX and YC regions.

Our manuscript shows that the simulation and measured values fit better, while the numerical range of the observed is greater than the simulated value (Fig.4). The possible reason is that the river entrance observation site is located by the lake, and the terrain is flat, especially the geographically adjacent between them, which reflects the lack of sensitivity of the model simulation as well.

The InVEST model has advantages in the spatial expression of NOP. Combined with the 2018 TLB Health Report (http://www.tba.gov.cn/slbthlyglj/), it was pointed out that the high-value areas of the pollution load of total nitrogen into the lake are HX and HJ regions, which is consistent with the spatial distribution shown in our manuscript. Shallow lakes in the above areas, such as Caohu and Chenghu Lake, need to pay special attention to their ecological environment.

The digital results of NOP are consistent with Bai et al. ([1]) and Chen et al. ([4]), which pointed out that the total nitrogen output of the TLB in 2000, 2005, and 2010 was 29683 t, 28,432 t, and 27,119t, respectively. However, the total nitrogen output in the TLB simulated by the SWAT model is quite different. For example, the manuscript pointed out that the NOP in the TLB was 33–66t at the river outlet estimated from 1995 to 2012 (Xu et al. [43]). The reason is speculated that the mechanism of the model is different. The InVEST model is mainly based on the LULC type to reflect the level of NOP and lacks detailed descriptions of nitrogen in water and terrestrial human production activities. At the same time, InVEST emphasizes the spatial distribution of NOP in the entire watershed, while SWAT focuses on changes in the content of substances in the river.

The amount of chemical fertilizer per unit area of cropland shows PD and PX regions up to 83087 kg/km2, followed by HJ, WC, HX, and YC region which is 38,318 kg/km2, 33,790 kg/km2, 31,448 kg/km2, and 29,910 kg/km2, respectively. The spatial distribution of this value is similar to the spatial distribution of NOP. It can be found that agricultural nitrogen fertilizer application is still the main source of nitrogen (Li et al. [26]). In fact, the NOP value has been significantly reduced. The main reasons may considered the absorption of crop and vegetation, soil retention, and water purification (Jabłońska et al., 2020).

The InVEST model has certain uncertainty when calculating NOP. First, the model is mainly based on the LULC data and empirical parameters. Although factors such as topography and precipitation are considered, the determination of empirical parameters is subjective, and they mainly come from the reference values provided by the model manual (Chen et al. [4]). Second, the quantitative description of the process of runoff generation and convergence is straightforward, which lack of consideration of various forms of pollutants, compared with model, e.g., SWAT. Third, the model lacks consideration of factors such as plant growth and microbial activities, as well as the process effect by temperature or precipitation (Chen et al. [4]).

Although the InVEST model has uncertain deficiencies in the description of the process mechanism, it can also reflect the spatial distribution characteristics of NOP in the TLB and provide a data reference for the treatment of pollutant sensitive areas.

The research revealed that the cropland had positive effects on NOP and forest had negative effects. These results were consistent with those of most previous study (Zhang et al. [49]). Cropland are still the main factor in the deterioration of water quality in major agricultural production areas, and agricultural non-point source pollution has always been one of the hot spots of concern (Peng and Li [31]).

However, the response of landscape metric to NOP shows more details in regional heterogeneity; for example, in the unit of hydraulic subregions, each region has at least two land use types that show the response relationship to NOP. Diverse LULC type combination expresses disparate in response to NOP.

Included large area of the cropland and developed (i.e., PD, PX, WC and HJ) had significant impact on NOP. The above-mentioned areas show some commonalities in terms of social, economic and terrain. They mainly include Shanghai municipality (including PD and PX regions), Wuxi City (located in WC regions), and Jiaxing City (located in HJ regions). In terms of social economy, the population density of Shanghai is 3,830 (people/km2), the highest GDP is 385.530 billion yuan (about 3.8% of the China national GDP), and the per capita GDP is 157,279 yuan. The population density of Wuxi City is 1,376 (people/km2), and the GDP is 179,848 billion yuan. The population density of Jiaxing City is 861 (people/km2), the GDP is 537.032 billion yuan, and the per capita output value is 147,649 yuan. NOP is closely linked to the rapid development of the regional population and social economy (Li et al. [25]). Historical data showed that TN concentration was significantly associated with populations of humans and GDP, which could explain 72% and 70% as explanatory variance (p < 0.05), respectively (Li et al. [25]). In terms of topography, Shanghai, the eastern plain of Wuxi City, and Jiaxing City have an average elevation of 3–5m, with dendritic drainage patterns.

Although there are certain commonalities in before-mentioned areas, there are still differences between the specific landscape indexes. In WC, PD, and HJ regions, the PLAND of cropland and developed land performs a positive impact on NOP, while only LPI of cropland shows a positive impact in PX region. Compared to the PLADJ, using LPI is more sensitive to NOP changes, which shows that the distribution of cropland of PX region is more concentrated. In the negative impact indicators of NOP, PLADJ, SPLIT, and CONNE of developed land reflect the heterogeneity between patches, which has a restrictive effect (Su et al. [33]). Compared with cropland as a high-intensity nitrogen source, the heterogeneity of developed land can reduce the total amount of nitrogen source and block the transport path (Wang et al. [40]).

The responded LULC types to NOP include cropland, developed land, and forest in ZX region, while HX and YC regions are cropland and forest. Compared with the response of a single crop landuse in the watershed, while forests landuse show sensitive responses to NOP after zoning. Forest as the main nitrogen sinks, which has an inhibitory effect on NOP (Ding et al. [6]; Groffman et al. [15]). There are different requirements for forest configuration (Hilary et al. [19]). In YC region, the PLAND of forest reflects the constraint effect to NOP, which should focus on the area of forest, while ZX region should focus on increased fragmentation of forest. HX region should pay attention to the shape of forest.

The distribution of land types in space and area seemingly determines the factors which can respond to the dependent variable. Regression analysis can only analyze based on the existing landscape types, which have enough area, whereas it cannot consider the impact of landscape types that are not enough or only have a small area.

Results reveal that cropland is still the main LULC type whether on a watershed scale or at a subregional scale. Nevertheless, compared with HX, ZX regions can reflect the forest landscape metric for NOP which has a negative effect. We need to pay more attention to ecological lands, such as the area of forest landuse in WC, PD, PX and HJ regions, even it is not shown in the regression equation. When the types of land are gradually enriched, especially increasing the proportion of non-agricultural land and natural land will improve regional water purification services and achieve regional sustainable development (Duarte et al. [10]).

The results of the manuscript may bring enlightenment as follows: (1) focus on the diversity of land use types. Diverse land use types have certain benefits in reducing the transmission path of nutrients, increasing the decomposition of microorganisms, and weakening the accumulation of pollution sources (Duan et al. [8]). In the process of controlling the pollution of nitrogen output, while taking measures to increase the area percentage of nitrogen sinks, the intensity of pollution sources per unit area should also be regulated (Wang et al. [40], [41]). (2) Concerned about the management of hydraulic subregions, which involving comprehensive cross-border regulation and control. Hydraulic subregion is a natural unit formed under watershed topography and natural and hydrological conditions, which can advantageously explain the migration path of pollutants. However, some hydraulic subregion units often contain at least two municipal administrative regions. This requires managers to coordinate and cooperate with different administrative districts within the same water conservancy unit. Definitely, the relationship between upstream and downstream also requires coordination and cooperation.

There is some uncertainty in the results of the manuscript, mainly summarized into two aspects: (1) the uncertainty of the data. This includes DEM data for subbasin extraction. The higher the accuracy of the data, the more detailed the subbasin units can be divided. In terms of nitrogen output simulation, the InVEST model input data, such as meteorological data using point interpolation, cannot accurately reflect the regional precipitation and temperature conditions. The verification data of nitrogen output is the automatic monitoring data deployed by the government, which may be affected by many factors such as weather. (2) Uncertainty of models and algorithms. The InVEST model reflects the macro-scale ecosystem services, which works as an ecosystem service reflecting the macro-scale. No matter from the research scale or the description of the object, it is difficult to accurately describe the process. It also brings challenges to the verification of the model results. In terms of watershed description, the dynamics of the river network are ignored in the algorithm of subbasin extraction in the plain river network area.

Meanwhile, this manuscript can only provide attribution explanations for NOP in terms of landscape composition and configuration, but their considerations are incomplete. Therefore, we hope to improve it from the following points in future research: (1) Natural climate factors such as temperature and rainfall, topographical factors, and socio-economic factors were not included in the correlation analysis (Su et al. [34]; Wan et al. [35]; Xu et al. [44]; Li et al. [26]; Duan et al. [7]; Wang et al., [39], [40], [41]). (2) In time scale, we need to consider the diversity in the driving effect of landscape pattern metrics under the long-term sequence with years as the time unit in the future. When months or seasons as the time unit are taken, the effect of changes in landscape pattern on NOP under natural landscape changes should be considered, such as vegetation replacement (Guo et al. [16]; Krausfeldt et al. [24]). (3) The diversity in pollutant transport under different landscape types and layouts are attempted to be answered, rather than simply defining the relationship as linear regression (Huang et al. [23]; Liu et al. [27]). (4) Based on the analysis of the landscape pattern metrics, a landscape amendatory plan should be proposed which could reach the water quality standard auxiliary. At the same time, managers and planners should consider the impact of the landscape metrics comprehensively and pay attention to the trend of the spatial distribution and propose measures to improve the ecosystem services.

Based on the LULC data of the TLB in 2020, our research calculated the typical landscape pattern metrics of the TLB. The InVEST model was used to spatially calculate the NOP. Based on the 125 subbasin units, the multiple stepwise regression model is used to quantitatively describe the relationship between landscape metric and NOP. The research results showed that the average value of NOP was 146.14 kg/km2, and the total output reached 23677.92t. Use observations to verify the simulation results, and the fitting coefficient is 0.76, the RMSE is 34.81, and the MAE is 68.39. At the subregional scale, the response effect of the landscape metric on NOP is better than that of the whole watershed scale. The main response land types include cropland, developed land, and forest. Among them, the PLAND and LPI of cropland show positive effects, the ED and NP of forest show negative effects, and the effect of developed land is diverse depending on the region.

Conceptualization: Ya'nan Wang and Guishan Yang. Methodology: Ya'nan Wang. Formal analysis and investigation: Ya'nan Wang. Writing—original draft preparation: Ya'nan Wang. Writing—review and editing: Ya'nan Wang, Guishan Yang, Bing Li, Chun Wang, and Weizhong Su. Funding acquisition, resources, and supervision: Guishan Yang.

This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA23020201).

Data can be provided upon request from the corresponding author.

Not applicable.

Not applicable.

The authors declare no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.