Combining the multivariate statistics and dual stable isotopes methods for nitrogen source identification in coastal rivers of Hangzhou Bay, China.

Coastal rivers contributed the majority of anthropogenic nitrogen (N) loads to coastal waters, often resulting in eutrophication and hypoxia zones. Accurate N source identification is critical for optimizing coastal river N pollution control strategies. Based on a 2-year seasonal record of dual stable isotopes ([Formula: see text] and [Formula: see text]) and water quality parameters, this study combined the dual stable isotope-based MixSIAR model and the absolute principal component score-multiple linear regression (APCS-MLR) model to elucidate N dynamics and sources in two coastal rivers of Hangzhou Bay. Water quality/trophic level indices indicated light-to-moderate eutrophication status for the studied rivers. Spatio-temporal variability of water quality was associated with seasonal agricultural, aquaculture, and domestic activities, as well as the seasonal precipitation pattern. The APCS-MLR model identified soil + domestic wastewater (69.5%) and aquaculture tailwater (22.2%) as the major nitrogen pollution sources. The dual stable isotope-based MixSIAR model identified soil N, aquaculture tailwater, domestic wastewater, and atmospheric deposition N contributions of 35.3 ±21.1%, 29.7 ±17.2%, 27.9 ±14.5%, and 7.2 ±11.4% to riverine [Formula: see text] in the Cao'e River (CER) and 34.4 ±21.3%, 29.5 ±17.2%, 27.4 ±14.7%, and 8.7 ±12.8% in the Jiantang River (JTR), respectively. The APCS-MLR model and the dual stable isotope-based MixSIAR model showed consistent results for riverine N source identification. Combining these two methods for riverine N source identifications effectively distinguished the mix-source components from the APCS-MLR method and alleviated the high cost of stable isotope analysis, thereby providing reliable N source apportionment results with low requirements for water quality sampling and isotope analysis costs. This study highlights the importance of soil N management and aquaculture tailwater treatment in coastal river N pollution control.

Introduction

Increasing anthropogenic nitrogen (N) inputs have significantly elevated riverine N loading to coastal zones over the past decades, resulting in an expansion of eutrophication and hypoxia zones (dissolved oxygen <2 mg L−1) and red tides worldwide (Breitburg et al. 2018). As the most important conduit for anthropogenic N delivery from land to coastal waters, coastal rivers contributed more than 43.2 Tg N to the ocean in 2000, which is expected to further increase in the future (Izett and Fennel 2018; Breitburg et al. 2018). Given the considerable N exported to coastal waters, reduction of coastal river N pollution is essential for mitigating coastal water eutrophication (Jiang et al. 2020). Therefore, it is warranted to quantitatively assess and identify the main N pollution sources in coastal rivers for optimizing pollution control strategies.

In China, 16.6–37.5% of coastal waters experienced eutrophication and ~700 red tide events occurred due to excessive nutrients in 2011–2020, with 60% of these events occurring in the East China Sea (Ministry of Ecology and Environment of China, 2020). Large population density coupled with intensive agricultural, aquaculture, and industrial activities in coastal areas of southeast China cause excessive N loading to coastal rivers, estuaries, and oceans (Bu et al. 2011; Liu et al. 2018; Luo et al. 2018). The aggregation of multiple N sources in the coastal area of southeast China (Herbeck and Unger 2013) makes it difficult to quantitatively identify N sources in coastal rivers. Thus far, there have been only limited attempts to quantify the N sources in coastal rivers of China (Kang and Xu 2016; Wu et al. 2020; Herbeck et al. 2021).

Watershed-scale mathematical models (including lumped models and distributed models) have been widely developed to identify riverine N sources (Borah & Bera 2004; Chen et al. 2014; Wellen et al. 2015; Imani et al., 2019). However, distributed models usually require large amounts of data (including hydrological, hydrochemical records, meteorological, land use, topography) for calibration and validation, resulting in many challenges for their applications in many regions (Chen et al. 2014). For lumped models, although they can provide a simplified understanding of N sources, requirements for watershed attributes and N discharge/loads from wastewater treatment facilities and industries limit their widespread application (Wellen et al. 2015; Hu et al., 2019). To overcome several of the limitations associated with these mathematical models, the absolute principal component score–multiple linear regression model (APCS-MLR) was adopted to quantify contributions from different pollution sources based on river water quality data. For example, the APCS-MLR method was used to identify the major contributors of riverine N in urban (Meng et al. 2018; Shen et al. 2021), suburban (Yang et al. 2013b), agricultural (Yuan et al. 2020), and mixed land-use catchments (Su et al. 2011). However, the APCS-MLR method cannot effectively distinguish contributions from sources with similar pollution characteristics (e.g., domestic wastewater and soil N are often identified as a mixed source) (Pasten-Zapata et al. 2014).

The dual nitrate () stable isotope ) approach can be an efficient tool for assessing N transformations (e.g., nitrification and denitrification, Xia et al. 2017; Yi et al. 2017; Hu et al. 2019) and qualitatively identifying/tracing N sources (Xue et al. 2009; Ji et al. 2017; Hu et al. 2019; Jiang et al. 2020). The dual stable isotopic information coupled with Bayesian statistics (i.e., MixSIAR, Liu &Han 2021) can effectively quantify N sources with unique isotopic signatures (e.g., atmospheric deposition, soil N pools, chemical fertilizers, manure/sewage, Xue et al. 2009). However, the high analytical cost has limited the widespread application of the dual stable isotope-based MixSIAR model (irregular and sparse data sets), thereby hindering its rigorous adoption for N source apportionment (Yang et al. 2013a; Kim et al. 2015). Combining the multivariate statistics methods and dual nitrate () stable isotope approach or the dual stable isotope-based model might provide complementary information to enhance the accuracy of N source apportionment results while lowering the requirements/costs for water quality monitoring (Yuan et al. 2020). Several studies adopted the multivariate statistics approaches (e.g., PCA, CCA, PMF) to qualitatively identify potential pollution sources of groundwater or precipitation N, and then used the dual stable isotope-based source apportionment method to quantify contributions of different pollution sources (Cui et al. 2018; Yuan et al. 2020). However, these studies only qualitatively identified water N source by the multivariate statistics approaches, lacking of attempts for river water N source apportionments. Therefore, it is necessary to combine the quantitative multivariate statistics methods (i.e., APCS-MLR method) and the dual nitrate () stable isotope-based MixSIAR model to quantify riverine N sources.

This study collected a 2-year, seasonal record of dual isotopes ( and ) and associated water quality data that was subsequently used for a synergistic analysis employing APCS-MLR and the dual stable isotope-based MixSIAR model to quantify the primary N sources in two coastal rivers draining into Hangzhou Bay, southeast China. The major objectives of this study were to (i) assess the spatial and temporal variations of eutrophication potential of coastal rivers in Hangzhou Bay (Cao’e River and Jiantang River); (ii) determine and validate the primary N sources in two rivers by combining the APCS-MLR and the dual stable isotope-based MixSIAR model; (iii) provide specific strategies for efficient N pollution control. This study provides an integrated approach for quantifying N contributions from multiple sources in coastal rivers, thus providing theoretical support for developing effective N pollution control strategies.

Materials and methods

Study area

As an important inlet to the East China Sea, Hangzhou Bay comprises 11 county-level administrative divisions with a total area of 4800 km2 (Sun et al. 2016). Hangzhou Bay is one of the most severely polluted coastal waters in China and was the location of 60% of the extreme red tide events (max area >100 km2) in China during 2020 (Ministry of Ecology and Environment, 2020). Climate is subtropical monsoon with an average annual temperature and precipitation of 15–25°C and 1300–1600 mm, respectively (Zhou et al. 2020). The water year was separated into four hydrologically contrasting seasons: dry season (spring: March–May; winter: December–February) and wet season (summer: June–August, autumn: September–November). Over the past 60 years, tideland reclamation projects have spurred the development of aquaculture ponds in the area by 138% (Zhou et al. 2020; Qiu et al. 2021). The southern coast of Hangzhou Bay (area with red boundary, Fig. 1) has an aquaculture zone comprised of 572.9 km2 of ponds (accounting for ~12.0% of the total area). Major aquaculture species are fish, shrimp and crabs, with breeding areas of 9.1%, 26.2%, and 42.2%, respectively (Ni et al. 2020).

Fig. 1

Location of study area and the distribution of sampling sites in Cao’E River (CER) and Jiantang River (JTR) in Hangzhou Bay

The Cao’e River (CER) and Jiantang River (JTR) are two typical coastal rivers flowing into southern Hangzhou Bay. Aquaculture ponds averaged ~34% of the CER nearshore area (within 3 km of river channel, 72.77 km2), with developed lands and arable land contributing ~33% and ~24%, respectively. Primary land-use types in the JTR nearshore (accounting for 39.61 km2) were aquaculture ponds, arable land, and developed lands, comprising 77%, 3%, and 2%, respectively (Table 1). The major aquaculture species in CER and JTR were shrimp. Population density in the CER region was ~560 person km−2 and the JTR region has ~1300 person km−2 (Zhejiang Statistics Bureau 2020). In the CER and JTR regions, multiple-year average precipitation is ~1600 mm and ~1800 mm, respectively, with 40–45% of rainfall occurring in the summer season due to typhoon events (Fig. 2; Zhou et al. 2020).

Table 1

Land-use distributions for Cao’E River (CER) and Jiantang River (JTR) in Hangzhou Bay

	Total area(km2)	A (%)	D (%)	G (%)	P (%)	Others (%)
CER	72.8	24.1	32.9	9.4	33.6	<0.1
JTR	39.6	3.2	1.6	<0.1	77.4	17.9

A, arable lands; D, developed lands; G, grasslands; P, aquaculture ponds

Fig. 2

Daily precipitation in Cao’E River (CER) and Jiantang River (JTR) catchments in Hangzhou Bay. The dotted line indicated a daily precipitation of 10 mm (the threshold of moderate rain)

Field sampling and water quality and stable isotope analysis

We collected water samples from four sites (XS, CXF, CXW, and PH, Fig. 1) around Hangzhou Bay for a total of 4 seasonal sampling events in 2019. Basic water quality indicators were measured, including total nitrogen (TN), total phosphorus (TP), ammonia (), nitrate (), nitrite (), dissolved organic carbon (DOC), chemical oxygen demand (CODMn), chlorophyll-a (Chl-a), pH, dissolved oxygen (DO), and transparency (Trans). Water samples for and analyses were collected seasonally in the CER and JTR from 2019 to 2020. Notably, water samples were only collected in summer (June) and autumn (September) in 2020 due to COVID-19 restrictions (i.e., water samples were collected in 6 seasons, 47 water samples including surface water, aquaculture tailwater, and rainfall). Samples were collected uniformly at 3–5 points across the river channel and mixed to obtain a single depth-width integrated sample. All water samples were immediately sealed and stored in the dark at 4°C.

The , , and concentrations were measured after filtering through a sterile 0.45-μm millipore polycarbonate membrane, whereas DOC was determined after filtering through a sterile 0.45-μm millipore organic phase membrane. and were determined using the spectrophotometric salicylic acid and phenol disulfonic acid methods, respectively. DOC was measured with a Jena multi-N/C 3100 analyzer. TN concentrations were analyzed as following alkaline potassium persulfate digestion. TP concentrations were measured by potassium persulfate digestion method and CODMn was determined by potassium permanganate oxidation method (State Environment Protection Bureau of China 2002). Transparency and chlorophyll-a (Chl-a) were measured by Secchi disk and chlorophyll-a fluorescence detector (FuloroQuik, AMI, USA), respectively. Temperature, pH, and DO were measured in the field using YSI sensors (Xylem, New York, USA). All water samples were analyzed within 24 h of collection.

Water samples for and analyses were stored at −20°C and subsequently transported to the Environmental Stable Isotope Lab (Chinese Academy of Agricultural Sciences, Beijing, China) for analysis. was converted to N2O by the denitrifier method (Christensen and Tiedje 1988), and then detected by a continuous-flow isotope ratio mass spectrometer (IRAM, Isoprime 100). Isotope ratio values are reported in parts per thousand (‰) relative to atmospheric N2 and Vienna Standard Mean Ocean Water (VSMOW) for  and , respectively (Xue et al., 2009). Sample analysis had an average precision of 0.31‰ for and 0.55‰ for .

All stable isotope results were expressed as delta values by representing the deviations in per mil (‰) from the respective reference standard as follow:

where Rsample and Rstandard are the isotopic ratios (15N/14N, 18O/16O) for the sample and standard, respectively.

Trophic level index (TLI)

To assess the river eutrophication status, we calculated TLI values based on measured water quality parameters (China Environmental Monitoring Station, 2001). TLI is a weighted sum based on correlations between Chl-a and various water quality parameters, and was calculated based on the concentrations of Chl-a, TP, TN, Trans, and CODMn (China Environmental Monitoring Station 2001):

where Wj is the weight of each index parameter and TLI(j) is the calculated index based on each pollutant concentration; r is the correlation coefficients between the reference Chl-a and each parameter j (Chl-a: 1, TP: 0.84, TN: 0.82, Trans: −0.83 and CODMn: 0.83), which were obtained from water survey data sets from China (Jin et al. 1995; Wang et al. 2019). TLI values range from 0 to 100 and characterize the eutrophication status based on five categories: oligotrophic (TLI<30), mesotrophic (3070) (Wang et al. 2019).

Nitrogen source identification methods

Absolute principal component score–multiple linear regression method

The absolute principal component score–multiple linear regression (APCS-MLR) method is derived from principal component analysis (PCA), and can be applied to identify pollutant sources (Yang et al. 2013a; Yuan et al. 2020). The overall PCA can be expressed as follow:

where Z is the normalized value of element i concentration, g represents the concentration of i in pollution source k (i.e., factor load), and h is the contribution of pollution source k to the sample j (i.e., factor score).

Based on the rotation factor loadings and principal eigenvectors obtained by the PCA, the multiple linear regression model (MLR) was employed to estimate the contribution of each source to . To quantitative source contributions, the normalized factor scores estimated by the PCA should be rescaled and converted to un-normalized component scores (i.e., APCS). The principle of the APCS-MLR model can be expressed as follows:

where C is the concentration of , r represents a constant, indicating the contribution from the pollution sources undetermined by PCA, r is coefficient of regression of source k, APCS is the absolute principal component scores of source k. The ratio of standard coefficient (r × APCSk) for each source and mean C represents the contribution to . The determination coefficient (R2) is obtained from the regression between modeled values (C) and observed nitrate concentration to test the efficacy of APCS-MLR model.

The dual stable isotope-based MixSIAR model

The Bayesian tracer mixing model (e.g., SIAR) has been applied to identify contributions of multiple N sources in individual water bodies (Stock et al. 2018; Hu et al., 2019). In recent years, an inclusive, rich, and flexible Bayesian tracer mixing model (i.e., MixSIAR) was further developed to include fixed and random effects as covariates explaining variability in mixture proportions (Stock et al. 2018; Liu and Han 2021). The MixSIAR (Ver. 3.0.2) was applied to apportion the contribution rates using R software (Ver. 4.0.3).

As MixSIAR needs pre-determination information for pollution sources, we adopted the results from APCS-MLR as the preliminary sources. Combining the field survey, land-use distribution and results from APCS-MLR (Table 1), this study apportioned riverine N to atmospheric deposition, soil N, aquaculture tailwater, and domestic wastewater. Aquaculture tailwater and atmospheric deposition N isotope information were directly measured, while isotopic values for soil and domestic wastewater N were obtained from published values for a nearby watershed (Table 4; Ji et al. 2017).

Table 4

, and isotope characteristic values of different pollution sources as well as observed values in water samples (‰) in coastal rivers of Hangzhou Bay

	Mean δ15N	SD of δ15N	Mean δ18O	SD of δ18O	Sources
Soil N	2.18	2.59	0.63	2.01	Ji et al. 2017
Domestic wastewater	10.49	4.53	3.45	2.63
Atmospheric deposition	−1.59	0.56	54.3	6.40	Measured
Aquaculture wastewater	2.83	3.07	7.21	2.60
River water in CER	5.34	3.35	4.02	3.98
River water in JTR	5.11	4.20	4.22	3.52

SD, standard deviation

The principle of MixSIAR can be expressed by the following formula:

where X is the tracer of species j of sample i, p represents the number of assumed sources, f is the tracer value for species j in source k, g denotes the contribution of source k to sample i, and e is an error term.

The PCA and APCS-MLR analyses as well as other relevant regression and correlation analyses were performed with SPSS software (Ver. 17.0, SPSS, IL, USA). One-way ANOVA with an LSD multiple comparisons test was used to examine the spatio-temporal variability of water quality parameters.

Results and discussion

Spatio-temporal variability of water quality

According to China’s national water quality standards (GB 3838-2002, Ministry of Environmental Protection of China, 2002), most of the coastal rivers discharging to Hangzhou Bay are considered severely polluted. Several pollutants concentrations exceeded the standard of “Grade IV” (including TP and CODMn, Table 2), which means severely polluted and can only be used for industry. Although no national river water quality standard is established for TN, TN concentrations (3.12±2.14 mg N L−1) were much higher than the 2.0 mg L−1 value classified as “Grade V” for reservoirs/lakes (Table 2), which is unable to support ecological functions. Riverine DO concentrations (7.13±2.34 mg L−1) were high due to natural reaeration by running water, as well as the artificial aeration of aquaculture discharge waters before releasing to rivers. Nutrient pollutant (TP, TN, , and ) concentrations in aquaculture ponds were all much higher (1–2 times) than those measured in river water (Table 2), implying that aquaculture discharge waters were a potential pollution source for coastal rivers.

Table 2

Mean and standard deviation of water quality data for all sampling sites in coastal rivers of Hangzhou Bay

	TP(mg L−1)	TN(mg L−1)	NH4+(mg L−1)	NO3-(mg L−1)	NO2-(mg L−1)	DOC(mg L−1)	COD(mg L−1)	Chl-a(μg L−1)	pH	DO(mg L−1)	Transparency(cm)	TLI
PH n=4	0.38±0.31c (0.10–0.89)	3.27±2.06b (1.64–6.80)	0.53±0.64b (0.11–1.63)	1.07±0.74b (0.16–2.07)	0.06±0.06b (0.03–0.17)	11.5±9.6b (2.8–27.8)	6.72±1.39c (4.44–7.92)	15.5±15.3c (1.3–40.8)	7.76±0.41a (7.44–8.45)	6.30±2.99a (1.65–8.92)	30.0±0.0c (30.0–30.0)	62.3±5.7b (53.8–69.9)
CXW n=4	0.23±0.18c (0.04–0.46)	2.32±1.19b (0.81–3.92)	0.93±0.74b (0.20–2.16)	0.80±0.41b (0.17–1.26)	0.10±0.14a (0.01–0.34)	10.4±6.8b (3.9–21.1)	8.78±1.82c (5.70–10.40)	57.7±42.6a (25.0–131.0)	8.05±0.40a (7.45–8.59)	6.67±1.73a (4.67–8.51)	31.3±2.2bc (30.0–35.0)	65.5±4.6b (59.3–72.3)
CXF n=4	1.21±0.27a (0.90–1.64)	6.58±5.11a (3.11–15.4)	1.09±1.42a (0.22–3.57)	3.88±3.77a (0.28–10.20)	0.18±0.22a (0.01–0.56)	15.8±16.9a (4.7–44.9)	8.33±1.89c (6.34–11.10)	33.3±11.1b (15.8–46.7)	8.46±0.62a (7.49–9.21)	6.66±2.05a (3.41–8.69)	37.5±5.6b (30.0–45.0)	72.9±1.4a (70.6–74.2)
XS n=4	0.11±0.07d (0.02–0.20)	0.85±0.23c (0.56–1.14)	0.37±0.20c (0.06–0.60)	0.35±0.09c (0.18–0.40)	0.01±0.01b (0.00–0.03)	7.7±1.3c (6.4–9.5)	9.31±2.24b (6.97–11.70)	32.7±15.3b (12.4–49.4)	8.50±0.19a (8.36–8.83)	8.89±1.22a (7.45–10.54)	65.0±5.0a (60.0–70.0)	57.5±4.0c (52.0–61.9)
CER n=16	0.49±0.23c (0.15–1.13)	3.08±0.96b (1.50–4.80)	0.72±0.24b (0.30–1.20)	1.92±0.61b (1.10–3.78)	0.02±0.02b (0.00–0.04)	–	16.14±8.77a (0.80–28.90)	–	–	–	–	–
JTR n=17	0.74±0.41b (0.11–1.34)	3.04±0.84b (2.10–5.30)	0.56±0.27b (0.25–1.50)	1.82±0.39b (1.49–3.30)	0.01±0.01b (0.00–0.05)	–	15.07±3.24a (8.00–19.80)	–	–	–	–	–
River water n=49	0.58±0.41 (0.02–1.34)	3.12±2.14 (0.56–15.40)	0.67±0.57 (0.06–3.57)	1.75±1.44 (0.16–10.20)	0.04±0.09 (0.00–0.56)	–	13.2±6.49 (0.80–39.30)	–	–	–	–	–
Aquaculture ponds n=6	1.32±0.30** (0.97–1.82)	5.88±0.41** (4.20–8.20)	1.36±0.22** (1.10–1.80)	3.20±0.87** (2.3–4.9)	0.11±0.03* (0.08–0.17)	–	27.15±3.28** (21.3–31.6)	–	–	–	–	–

The numbers in the brackets represent the range of each parameter

Superscript letters after numbers denote significant differences (threshold is p<0.05, among PH, CXW, CXF, XS, CER, and JTR)

*Denotes significant differences between river waters and aquaculture pond waters (*p<0.05, **p<0.01)

=16

=49

=6

Generally, water quality across rivers showed considerable spatio-temporal variability (Table 2, Fig. 3). Pollutants concentrations (including TP, TN, , , , DOC, and chlorophyll-a) at site CXF were 2.0–15.3 fold higher than the other sites. Site XS showed the lowest concentrations for most pollutants (including TP, TN, , , , and DOC, Table 2). Concentrations of  and exhibited similar spatial variations (i.e., highest in CXF and lowest in XS). The higher and concentrations observed at CXF are attributed to the higher population density (1300 person km−2, ~150% higher than other sites) and higher precipitation (Fig. 2). This implies that soil N and domestic wastewater might be important N sources to CXF. Most of the sampling sites were characterized as moderately-eutrophic (6070), followed by CXW, PH, and XS (Table 2).

Fig. 3

Seasonal dynamics of water quality parameters in river sampling sites PH, CXW, CXF, and XS in coastal rivers of Hangzhou Bay

In terms of temporal variability, TN (p<0.01), (p<0.01), and (p<0.05) concentrations were significantly higher in the winter/spring dry season than in the summer/autumn wet season (Fig. 3). Moreover, the higher N concentrations in winter/spring may be associated with higher aquaculture tailwater discharge following the aquaculture harvest (Ni et al., 2021), as well as lower denitrification due to lower temperatures (Jiang et al. 2020). Conversely, no significant seasonal variations were observed in TP and CODMn concentrations across all sites (Fig. 3), implying differential processes affecting the various pollution sources (i.e., aquaculture tailwater, domestic wastewater, and agricultural runoff) and their transformations. TLI values for XS, PH, CXW, and CXF were higher in summer or autumn (Fig. 3), which was mainly associated with seasonal variations of Chl-a as regulated by temperature and light conditions (Herb and Stefan 2003; Paerl and Paul 2012; Jiang et al. 2018).

Riverine N source apportionments

Source apportionment by the APCS-MLR method

Estimated sphericity values based on Kaiser-Meyer-Olkin (KMO) and Bartlett’s test for the PCA incorporating 9 hydrochemistry parameters (TP, TN, CODMn, pH, DOC, , DO, , ) were 0.580 and 212 (p<0.05), respectively. These values support the efficacy of the PCA approach for pollution source identification (Mei et al. 2014). Three principal factors were identified by PCA explaining a total variance of 84.5% (Table 3). In detail, factor 1 explained 50.3% of total variance with strong positive loadings on TN (0.943), TP (0.282), and (0.937) and a moderate negative loading on DO (−0.164). Strong positive loadings on and moderate negative loading on DO usually indicated pollution contribution of domestic wastewater discharge (Yang et al. 2013a; Mei et al. 2014), while strong positive loadings on TN and TP can be interpreted as pollution contribution of agricultural runoff (Mei et al. 2014; Yuan et al. 2020). Therefore, factor 1 implied mixed impacts from multiple pollution sources of domestic wastewater and agricultural runoff.

Table 3

Loadings of the water quality parameters on the principal components in coastal rivers of Hangzhou Bay

	Factor 1	Factor 2	Factor 3
TP	0.282	−0.139	0.913
CODMn	0.316	0.697	0.078
DO	−0.164	0.786	−0.302
TN	0.943	−0.008	0.06
pH	−0.324	0.652	0.621
DOC	0.915	0.162	−0.168
NH4+-N	0.937	0.209	−0.195
NO2--N	0.954	−0.071	−0.039
NO3--N	0.839	0.132	0.222
% of variance	50.3	18.2	16.0
Cumulative %	50.3	68.5	84.5
Contributions to NO− 3	0.695	0.222	0.083

The absolute values of bold values exceeding 0.50 are the dominant components of each factor

Factor 2 explained 18.5% of the total variance with strong positive loadings on DO (0.786), pH (0.652), and CODMn (0.697) and moderate correlation with (0.209). Factor 2 characteristics are consistent with a pollution source from aquaculture tailwaters. The strong positive loading on , pH, and DO is attributed to the usage of fertilizers (e.g., organic fertilizer, ammonia fertilizer, compound fertilizer) (Herbeck et al. 2021) and artificial aeration activities in aquaculture management practices (Ni et al. 2020). The high loading on CODMn was associated with large organic food inputs and excreta of aquatic animals, as well as organic-rich sediments in the tailwater (Wang et al. 2020; Kim et al. 2021).

Factor 3 explained 16% of the total variance with strong positive loadings on TP (0.913) and pH (0.621), as well as negative loadings on DO (−0.302) and  (−0.195). We ascribe these loading characteristics to point source pollution (Mei et al. 2014). The low loading on may result from efficient domestic sewage treatment and high nitrification potential during transport (Yang et al. 2013a; Yuan et al. 2020). In summary, the PCA results identified pollutant sources in the rivers originating from soil, domestic wastewater (including untreated and well treated sewage), and aquaculture tailwater.

The APCS-MLR method (i.e., Eq. 5) was further applied to quantify the contribution of each N source to riverine N loads. The reasonably strong coefficient of determination (R2 = 0.845) between the APCS-MLR model simulated and observed concentrations supports the efficacy of the APCS-MLR model. We interpret the results to indicate that a mixed source (soil N and untreated domestic wastewater) was the primary N source to coastal rivers (69.5%). Additional sources were assigned to aquaculture tailwaters (22.2%) and residual sources (8.3%). This quantitative information highlights the importance of controlling pollution from wastewater and soil erosion/leaching in coastal watersheds. It is important to note that limitations associated with the APCS-MLR method prevent further characterization of the mixed source into higher resolution components, which limits the use of the results for the development of effective pollution control strategies.

Source apportionment by the dual stable isotope-based MixSIAR model

The overall and values were similar (p > 0.05) between JTR (4.22 ± 3.52‰ and 5.11 ± 4.20‰) and CER (4.02 ± 3.98‰ and 5.34 ± 3.35‰), respectively (Fig. 4, Table 4). The ranges of and signatures were overlapping between soil N, domestic wastewater, aquaculture tailwater, and river water, suggesting that all of these sources might contribute to riverine N loads (Hu et al. 2019).

Fig. 4

Scatter plot of and values in rivers (2019 and 2020) and aquaculture ponds in Hangzhou Bay. Dotted boxes show the typical ranges of and values from Xue et al. (2009), Yi et al. (2017), Ji et al. (2017), and Hu et al. (2019)

More than 70% of the values fell within the range of formed as a product of nitrification (−1.92 to 7.44‰; calculated based on the δ18O values of atmospheric O2 and precipitation; Ji et al. 2017). The slope between and for 2019 and 2020 samples were 1.11 (n=19, R2=0.54, p<0.01) and 1.45 (n=12, R2=0.48, p<0.01), respectively. These values exceeded the range of slopes expected for N experiencing considerable denitrification (i.e., slopes between 0.43 and 1.0, Fig. 4, Xue et al. 2009; Hu et al. 2019). The absence of a strong denitrification signal is consistent with the high DO concentrations in river waters (7.12 ± 2.34 mg L−1, Table 2). These high DO concentrations would be expected to hinder the denitrification process (<2 mg L−1 favors denitrification, Rivett et al. 2008). Overall, nitrification rather than denitrification was the likely dominant N transformation process in CER and JTR.

Similar results were observed between aquaculture pond and river waters, with more than 60% of   values falling within the range indicative of  formation as a product of nitrification (−1.92 to 7.44‰). No significant relationship (p=0.90, n=15) was observed between and for aquaculture pond waters (Fig. 4). Collectively, these results infer that artificial aeration of aquaculture ponds maintains adequate DO to promote nitrification and inhibit denitrification (Jiang et al. 2020; Ni et al. 2020). Thus, aquaculture tailwaters may be an important source of pollution in both CER and JTR.

According to the field investigations (Table 1), PCA results (Table 3) and observed ranges of δ18O-NO3− and δ15N-NO3− values (Table 4 and Fig. 4), soil N, domestic wastewater, aquaculture tailwater, and atmospheric deposition were identified as four major N pollution sources for river. Using the original and values of these sources (Table 4), the dual stable isotope-based MixSIAR model was adopted to quantify contribution of each pollution source. For JTR, the dual stable isotope-based MixSIAR model estimated that soil N (34.4 ± 21.3%) was the dominant riverine N source, followed by aquaculture tailwaters (29.5 ± 17.2%) and domestic wastewater N (27.4 ± 14.7%). The contribution from atmospheric deposition was deemed relatively low (8.7 ± 12.8%, Fig. 5c). Similar results were determined for CER, with riverine N source contributions from soil N, aquaculture tailwater, domestic wastewater, and atmospheric deposition contributing 35.3 ± 21.1%, 29.7 ± 17.2%, 27.9 ± 14.5%, and 7.2 ± 11.4%, respectively (Fig. 5b).

Fig. 5

The contribution of atmospheric deposition N (AD), soil N (SN), aquaculture tailwater (AQ), and domestic wastewater N (DW) in a coastal rivers of Hangzhou Bay; b Cao’E River (CER) and c Jiantang River (JTR) across different seasons based on the dual nitrate stable isotope approach (the MixSIAR model ran 3,000 simulations, i.e., n = 3000 for each box). Box plots in (a), (b), and (c) illustrate the 25th, 50th, and 75th percentile; the whiskers indicate the 2.5th and 97.5th percentiles

High N contributions from soil N pools are derived from long-term excess N applications (i.e., legacy N) that have accumulated in the soil. The average soil N contents were 2.53 g kg−1 in JTR and 2.40 g kg−1 in CER, values much higher than the China national average (1.45 g kg−1) for cropland soils (Zhou et al. 2020; China Soil Database). A large contribution from aquaculture is consistent with the high density of aquaculture ponds and the high N concentrations found in pond waters. Whereas pond waters are artificially aerated to avoid low DO concentrations in tailwaters, there are no specific management practices to address the release of N in tailwaters (Ni et al. 2021). High N contribution from domestic wastewater may be derived from high population intensity in both CER (~560 person km−2) and JTR (~1300 person km−2). Atmospheric deposition also contributes a considerable proportion of N due to excess N fertilizer application and fossil fuels burning (Zhou et al. 2020).

In terms of temporal variability, the contribution of N from aquaculture tailwaters was significantly higher in winter in both CER and JTR (p<0.01, Fig. 5b and c), which was consistent with the shrimp breeding cycle. Large amounts of N accumulate in the ponds during the breeding period (summer to autumn) and are subsequently discharged into the rivers during the winter following harvest (Ni et al. 2020). Similarly, the N contribution from domestic wastewater was highest in winter (p<0.01, Fig. 5), owing to point source pollution being most strongly expressed during the winter dry season due to the lack of dilution (Hu et al. 2018). Seasonal variation of soil N contributions in CER and JTR were similar, with the lowest contribution in winter and the highest contribution in summer/autumn (Fig. 5b and c). This seasonal pattern can be attributed to the seasonal distribution of precipitation (Fig. 2). Greater soil erosion/leaching occurs in the summer when ~45% of annual precipitation occurs (CER = 7.85 and JTR = 8.89 mm day−1), compared to only ~13% of precipitation in winter (CER = 2.48 and JTR = 2.38 mm day−1). Higher rainfall amounts and intensities were shown to result in greater soil N losses to rivers (Deng et al. 2019; Zhang et al. 2021). The contribution of atmospheric deposition (wetfall + dryfall) was relatively stable with no significant seasonal fluctuations (p>0.05).

Comparison between the two methods

The combined utilization of the APCS-MLR and the dual stable isotope–based MixSIAR model provided consistent results and a limited verification step for estimating riverine N sources to rivers (Table 3 and Fig. 5). The estimated N contribution from soil N and domestic wastewater from the dual stable isotope with the MixSIAR method (63.8%) was remarkably similar to the mixed source contribution calculated from the APCS-MLR method (69.5%). Similarly, the estimated N contributed from aquaculture tailwaters was consistent between the two methods (27.5% vs. 22.2%). These results provide compelling evidence that soil N, aquaculture tailwaters, and domestic wastewater were the major riverine N pollution sources.

Notably, the APCS-MLR method failed to distinguish between the contributions from soil N pool and domestic wastewater (Table 3). While the dual stable isotope-based MixSIAR model was able to identify soil N and domestic wastewater as distinct sources, the estimated results incorporate uncertainty due to lack of independent verification (Xue et al. 2009; Hu et al. 2019) and are often constrained by the high cost of isotope analysis. The APCS-MLR results provide the necessary source discrimination information for subsequent the dual stable isotope–based MixSIAR model applications, and reducing the cost of source screening for the stable isotope method. Furthermore, the two methods provide a measure of verification, which reduces the uncertainty associated with major source identification. Integration of the two methods provides an opportunity to achieve reliable N pollution source identification for rivers having limited data records for dual isotopes and water quality. Future studies incorporating relevant data from higher monitoring frequencies should be adopted to further evaluate the advantages/limits of the combined APCS-MLR model and the dual stable isotope–based MixSIAR model by testing performance as a function of temporal data resolution.

Implications for N pollution control in coastal rivers

Results of this study confirmed the high eutrophication status of coastal rivers in Hangzhou Bay (Table 2), which contributes to persistent red tide events in coastal waters of eastern China (Breitburg et al. 2018; Jiang et al. 2020). Quantitative N source identification/apportionment provides critical information for adopting specific strategies to mitigate riverine N pollution in receiving waterbodies. Considering the high pollution contribution from cropland soil N runoff to Hangzhou Bay (Fig. 5a), sustainable nutrient management strategies are necessary to optimize utilization of soil N by crops and to reduce excessive N fertilizer inputs and subsequent N losses to runoff/leaching (Hu et al. 2019). Moreover, soil and water conservation practices (e.g., conservative tillage, plant/residue soil cover, and terraces) are important approaches for mitigating N pollution from soil erosion (Chen et al. 2018). For aquaculture tailwaters (Table 3, Fig. 5), it is necessary to reduce N and other pollutant discharge loads by improving aquaculture management practices (Turcios and Papenbrock 2014). For example, an integrated multi-trophic aquaculture (IMTA) system can greatly reduce aquaculture tailwater pollutant discharges by recycling by-products of waste from one aquaculture species to fertilize/feed other aquatic species (Neori et al. 1998; Naylor et al. 2000; Rakocy et al. 2006; Martan 2008). Zhang and Kitazawa (2016) estimated that IMTA can achieve maximum reductions of 30% for phytoplankton and 35% for N in tailwaters. Additionally, flow-through wetlands consisting of oxidation and reduction ponds connected in series can effectively enhance nitrification and subsequent denitrification to remove inorganic N species. A large quantity of untreated domestic wastewater is often discharged directly into coastal rivers in the Hangzhou Bay region due to the lack of efficient collection and treatment facilities (~46% in 2019, Environmental Monitoring Center of Zhejiang, 2020). Our analysis indicated that domestic wastewater contributed considerable N to river waters, comprising a larger proportion of the riverine N load during the winter/spring dry period (Table 3, Fig. 5). Therefore, it is necessary to improve domestic waste collection for treatment.

Given the changing seasonal dynamics of riverine N sources (Fig. 5), management strategies must consider how specific sources change across different seasons. As eutrophication occurs disproportionately in summer and autumn due to higher temperature and light conditions (Paerl & Paul 2012; Jiang et al. 2018), management strategies should target N control measures to reduce riverine N concentrations especially during these two seasons. For example, effective soil and water conservation practices, as well as cropland nutrient management strategies, are warranted to reduce soil N losses during this summer/fall critical period for eutrophication. Correspondingly, more attention should be paid to domestic wastewater and aquaculture tailwater N reduction strategies during low flow periods (e.g., winter and spring) to reduce eutrophication potential.

Conclusion

This study combined the APCS-MLR model and the dual stable isotope–based MixSIAR model to quantitatively assess N dynamics and source contributions in coastal rivers discharging into Hangzhou Bay. This complementary and integrated approach achieved robust N source information with limited sampling requirements. Our water quality measurements confirmed that the studied rivers were severely polluted, with more than half experiencing moderate eutrophication (TLI>60). Spatio-temporal variability in water quality was associated with seasonal agricultural, aquaculture, and domestic activities, coupled with the seasonal rainfall pattern. The APCS-MLR model identified a mixed source (soil N + domestic wastewater) and aquaculture tailwater as the main N sources, with contributions of 69.5% and 22.2%, respectively. The dual stable isotope method identified considerable nitrification and minimal denitrification as N transformation processes. Furthermore, the dual stable isotope–based MixSIAR model apportioned riverine N contributions from soil N, aquaculture tailwater, domestic wastewater, and atmospheric deposition N, with contributions of 35.3 ± 21.1%, 29.7 ± 17.2%, 27.9 ± 14.5%, and 7.2 ± 11.4% in the Cao’e River (CER) and 34.4 ± 21.3%, 29.5 ± 17.2%, 27.4 ± 14.7%, and 8.7 ± 12.8% in the Jiantang River (JTR), respectively. Moreover, the APCS-MLR and the dual stable isotope–based MixSIAR model provided complementary/synergist information and a measure of verification for the quantitative apportionment results. Overall, the results highlight that coastal river N pollution control must focus on soil N management, aquaculture tailwater, and domestic wastewater treatment in the watersheds surrounding Hangzhou Bay.