## Abstract

A nonpoint source (NPS) loads evaluation method based on the Bayesian source apportionment mixing model was established in this research. The model assumed that (1) the pollutant concentration from each source mixed with the others in the monitoring section during transport, (2) transport only considered first-order attenuation, (3) point source pollution had relatively stable emissions, and (4) the measurement error was random, unrelated, and consistent with a normal distribution (mean of 0). All unknown parameters in the model were taken as random variables, and their posterior distributions were derived by Markov chain Monte Carlo procedures based on historical data, literature, and empirical information. The outflow system of the Huaihe River was adopted as a case study to verify the feasibility of the model. Gelman–Rubin, automatic frequency control statistics, and the determination coefficient (*R*^{2}) verified the reliability. The results showed that the total loads of ammonia nitrogen (NH_{4}^{+}), chemical oxygen demand, total nitrogen, and total phosphorus from NPSs accounted for 16.35–27.58%, 18.78–25.69%, 21.68–29.71%, and 42.11–52.09%, respectively. The parameter sensitivity analysis showed that prior distribution of NPS concentration was the most sensitive one, which should be determined reasonably based on the empirical or historical information.

## INTRODUCTION

Nonpoint source (NPS) pollution is often defined as diffuse pollution that occurs over a wide area and is associated with land use (Nikolaidis *et al.* 1998; Xu *et al.* 2013), and this form of pollution has had an important effect on water quality impairment worldwide (Chen *et al.* 2010; Zhuo *et al.* 2017). Traditional assessment methods for NPS pollution load estimation often demand a substantial amount of data consisting of basin data and river data, which are not always obtained during routine monitoring. Recently, an inverse modelling approach, which has been widely applied to environmental and hydrological problems (Woodbury & Ulrych 2000; Shen *et al*. 2006; Zou *et al.* 2007), has provided the possibility of estimating NPS loads based on limited data.

Currently, there are two simple and effective formulas used to estimate NPS loads with inverse models. The first one is the concentration formula (Liu *et al.* 2008; Gronewold *et al.* 2009; Shen & Zhao 2009, 2010), which estimates the NPS load via the concentration relationship between upstream and downstream regions. However, this formula lacks descriptions for the influences of point sources and discharge in the relationship. Therefore, the NPS load assessment results may be too high or too low. The second formula is the load formula (Chen *et al.* 2011, 2012; Zhao *et al.* 2014), which estimates the NPS load via load conservation and joins the point source and discharge descriptions. However, the requirements for upstream and downstream regions and point source flow monitoring are not easily met, and it is difficult to employ the formula for a wide range of applications. Thus, it would be beneficial to develop a robust method for load estimation by using incomplete routine monitoring data.

The source apportionment mixing model (SAMM) is often used to calculate the proportional contributions of pollution sources in a sample by using stable isotopes (Xue *et al.* 2012; Yang *et al.* 2013; Meghdadi & Javar 2017; Wang *et al.* 2017) or chemical elements (Massoudieh & Kayhanian 2013; Sharifi *et al.* 2014). The SAMM uses the characteristics of each source to identify and apportion the volumetric flow rate contributions of each source based on the mass conservation law. Based on this point of view, the SAMM can also be used to estimate the proportion of NPS pollution to recipient regions based on emission characteristics if the emission characteristics of the sources are not collinear. Furthermore, NPS loads can be estimated by using the contribution proportion and known point source load. This method can reduce the requirements for flow data and increase the applicability.

Due to the uncertainties of the model structure, parameters, and data (Lindenschmidt *et al.* 2007), it is obviously impossible to obtain the true values of the parameters. The Bayesian method, which treats all unknown variables as random variables, is usually applied to evaluate the parameters and uncertainties in an inverse model (Liu *et al.* 2008; Shen & Zhao 2010; Chen *et al.* 2011, 2012; Zhao *et al.* 2014). The Bayesian approach can make full use of empirical knowledge or historical data through prior distributions (Freni & Mannina 2010; Neuman *et al.* 2012); therefore, it can estimate parameters more accurately than traditional methods, especially when observation data are scarce. In addition, it addresses the uncertainty of the model parameters in a joint posterior distribution format (Massoudieh & Kayhanian 2013), which is convenient for analysing the risks of decision-making in water environment management (Malve & Qian 2006). Therefore, the Bayesian method is used to address the uncertainty of NPS loads in our model.

This study develops a new tool for estimating NPS loads from the point of view of source apportionment. This tool overcomes the disadvantages of needing more discharge monitoring requirements in traditional load formulas and can carry out NPS load assessments based on regular monitoring data to extensively support practical water quality management. The heavily polluted Huaihe River outflow system is invoked as a case study. The model results will help local decision makers quantify NPS loads and determine the key issue in water environment treatment.

## METHODS AND MATERIALS

### Flowchart of the Bayesian source apportionment mixing model

As shown in Figure 1, we simplify the model and turn it into a mathematical system by making a reasonable assumption of the physical conditions based on the conceptual model. The simulation error is introduced to transform the system into an expressed form of the random variable required in the Bayesian method, which generates the likelihood function of the model by replacing the observation data. The prior distribution of parameters is defined by historical data or experience-based information. Then, the Bayes’ theorem is used to express the posterior distribution of the parameters associated with the prior distribution and likelihood function, where the Markov chain Monte Carlo (MCMC) method is used to calculate the posterior distribution of parameters, which results in the final output of the NPS load evaluation.

### Source apportionment mixing model

*C*represents the concentration of the downstream receiver (mg/L);

_{d}*T*and

*T*represent the in-stream travel times (d), which are calculated from the mean velocity and flow-path length, respectively;

_{i}*p*represents the volumetric flow rate contribution of each source to the recipient sample (dimensionless), where the sum is 1;

*k*is the first-order decay rate of the pollution species (1/d);

*C*,

_{u}*C*, and

_{i}*C*represent the pollution concentrations at upstream, point, and NPSs, respectively (mg/L); and

_{n}*m*represents the number of point sources (dimensionless). Pollution from NPSs is assumed to be lost over half of the stream length on average (Alexander

*et al.*2006; Chen

*et al.*2011).

*L*represents the NPS pollution load of species

_{n,j}*j*,

*q*represents the point source discharge,

_{i}*C*represents the NPS concentration of species

_{n,j}*j*,

*p*represents the fractional concentration of pollution from the NPS, and

_{n}*p*represents the fraction of pollution from point source

_{i}*i*.The model in Equation (4) can be viewed as a forward model from

*p*to

_{n}*C*. The reverse model for Equation (5) is

_{j}Thus, the problem of load estimation has transformed into an inverse model with the purpose of fitting the set of *C _{n}*,

*p*,

_{u}*p*, and

_{n}*p*. Then, the Bayesian approach is used to calibrate the unknown parameters.

_{i}### Bayesian approach

*θ*represents the parameter that needs to be estimated;

*p*(

*θ|y*) represents the posterior probability of

*θ*, which also represents the strength of our belief in

*θ*when the data have been taken into account;

*p*(

*θ*) represents the prior distribution of

*θ*, which is defined by historical or experience-based information; and

*p*(

*y|θ*) represents the likelihood function, which defines the probability that the data could be generated by the model with a parameter value of

*θ*.

*ɛ*is added to the SAMM, as shown in Equation (7), with a zero mean and a variance of

_{j}*σ*

^{2}. Therefore, the likelihood function for all observed values is given in Equation (8), where

*z*represents the number of pollution species.

Because the posterior distribution is a complex, high-dimension integral that is difficult to evaluate analytically, the MCMC was utilized to solve this issue, as it is commonly applied in Bayesian analyses (Lunn *et al.* 2000; Gelman 2006; Malve & Qian 2006).

### Parameter sensitivity analysis

*et al.*1990; Yoshie

*et al.*2007). Sensitivity index (

*S*), the ratio of the fractional change in the model outputs to the fractional changes in the parameter value, was introduced to describe the influence of each parameter on the modelling outputs, which can be calculated by the following equation:where

_{p}*E*(

*p*) is the model outputs in the baseline run (i.e., all parameters at their calibrated values

_{b}*p*), and

_{b}*E*(

*p*

_{high}) and

*E*(

*p*

_{low}) are the model outputs for when the parameter is set to its high value

*p*

_{high}(for example, 10% higher than calibrated) and its low value

*p*

_{low}(for example, 10% lower than calibrated). Because our model is a Bayesian model where the parameter and outputs are distributions, so we use the mean of the distributions of the parameter and outputs to calculate the sensitivity index.

### Case study description

The Huaihe River outflow is located in the lower Huaihe River Basin in northern Jiangsu Province and is one of the five flood protection channels for Hongze Lake, as shown in Figure 2(a). The upper inflow from Hongze Lake is controlled by the Erhe floodgate. The outflow towards the East China Sea is limited by the Haikou floodgate. The Erhe floodgate has always been closed, except in July and August of 2003, July and August of 2007, and April and May of 2010, where only a small amount of leakage occurred when the gate closed. The river is 163.5 km in length, with a basin area of 1,710 km^{2}. The discharge of the Huaihe River mainly derives from rainfall runoff in the watershed and urban sewage water when the upper floodgate is closed. In recent years, the Huaihe River outflow has been heavily polluted by both point and NPSs. It has been listed as a key object for the prevention and control of water pollution in the Huaihe River Basin. Identifying the source of river pollutants is essential for prevention and control plans and provides a scientific basis for prevention and control measures.

Ten sets of hydrological sites (H1) and two water-quality monitoring sites (S1 and S2) from 30 November to 4 December 2016, and 11 April to 18 April 2017, were collected. Moreover, six sewage outlets with monitoring data during the same period were supplied by the Jiangsu Province Hydrology and Water Resources Investigation Bureau. Their spatial distributions are shown in Figure 2(a). Here, we generalized the river and point sources in the study area, as shown in Figure 2(b).

## RESULTS AND DISCUSSION

### Prior distribution determination and analysis

The observed data in Tables 1 and 2 for the variables considered in the model were initially analysed as references for the prior distribution of unknown parameters before the implementation of Bayesian statistics. As seen from Table 1, the water quality downstream (S2) is significantly worse than that upstream (S1), indicating that the interval flows into the river have a lower water quality than those upstream. The emission concentrations from point sources in Table 2 are far greater than those downstream (S2) in Table 1, indicating that the NPS inflow is cleaner than that of the point source. Therefore, the prior distribution of NPS pollutant concentrations is greater than that for upstream concentrations. From the above, the NPS pollutant concentrations were empirically determined as the normal distribution with the range of 0.09–8.14 mg/L for NH_{4}^{+}, 7.8–66.5 mg/L for chemical oxygen demand (COD), 2.09–12.4 mg/L for total nitrogen (TN), and 0.02–0.98 for total phosphorus (TP), in which the lowest values were taken from the upstream (S1) minimum and the highest values were taken from the downstream (S2) maximum.

Sample . | Term . | Median . | Minimum . | Maximum . | S.D. . |
---|---|---|---|---|---|

S1 | NH_{4}^{+} (mg/L) | 0.10 | 0.09 | 0.54 | 0.14 |

S1 | COD (mg/L) | 11.05 | 7.80 | 18.00 | 3.17 |

S1 | TN (mg/L) | 2.61 | 2.09 | 3.00 | 0.32 |

S1 | TP (mg/L) | 0.08 | 0.02 | 0.28 | 0.08 |

S2 | NH_{4}^{+} (mg/L) | 3.59 | 0.69 | 8.14 | 2.46 |

S2 | COD (mg/L) | 39.55 | 12.00 | 66.50 | 20.73 |

S2 | TN (mg/L) | 6.89 | 3.49 | 12.40 | 2.89 |

S2 | TP (mg/L) | 0.24 | 0.09 | 0.98 | 0.31 |

H1 | Velocity (m/s) | 0.75 | 0.21 | 1.15 | 0.38 |

H1 | Discharge (m^{3}/s) | 82.65 | 22.10 | 123.70 | 44.02 |

Sample . | Term . | Median . | Minimum . | Maximum . | S.D. . |
---|---|---|---|---|---|

S1 | NH_{4}^{+} (mg/L) | 0.10 | 0.09 | 0.54 | 0.14 |

S1 | COD (mg/L) | 11.05 | 7.80 | 18.00 | 3.17 |

S1 | TN (mg/L) | 2.61 | 2.09 | 3.00 | 0.32 |

S1 | TP (mg/L) | 0.08 | 0.02 | 0.28 | 0.08 |

S2 | NH_{4}^{+} (mg/L) | 3.59 | 0.69 | 8.14 | 2.46 |

S2 | COD (mg/L) | 39.55 | 12.00 | 66.50 | 20.73 |

S2 | TN (mg/L) | 6.89 | 3.49 | 12.40 | 2.89 |

S2 | TP (mg/L) | 0.24 | 0.09 | 0.98 | 0.31 |

H1 | Velocity (m/s) | 0.75 | 0.21 | 1.15 | 0.38 |

H1 | Discharge (m^{3}/s) | 82.65 | 22.10 | 123.70 | 44.02 |

Point source . | Discharge (m^{3}/s)
. | NH_{4}^{+} (mg/L)
. | COD (mg/L) . | TN (mg/L) . | TP (mg/L) . |
---|---|---|---|---|---|

PS1 | 0.82 | 1.24 | 86.52 | 10.45 | 0.17 |

PS2 | 1.06 | 2.10 | 61.67 | 4.23 | 1.22 |

PS3 | 4.04 | 4.95 | 62.15 | 7.38 | 0.85 |

PS4 | 2.20 | 44.85 | 471.10 | 66.40 | 1.16 |

PS5 | 0.79 | 1.97 | 136.59 | 17.03 | 1.77 |

PS6 | 0.26 | 4.39 | 97.30 | 8.03 | 3.08 |

Weighted average^{a} | 9.17^{b} | 124.59 | 1557.01 | 204.51 | 9.62 |

Point source . | Discharge (m^{3}/s)
. | NH_{4}^{+} (mg/L)
. | COD (mg/L) . | TN (mg/L) . | TP (mg/L) . |
---|---|---|---|---|---|

PS1 | 0.82 | 1.24 | 86.52 | 10.45 | 0.17 |

PS2 | 1.06 | 2.10 | 61.67 | 4.23 | 1.22 |

PS3 | 4.04 | 4.95 | 62.15 | 7.38 | 0.85 |

PS4 | 2.20 | 44.85 | 471.10 | 66.40 | 1.16 |

PS5 | 0.79 | 1.97 | 136.59 | 17.03 | 1.77 |

PS6 | 0.26 | 4.39 | 97.30 | 8.03 | 3.08 |

Weighted average^{a} | 9.17^{b} | 124.59 | 1557.01 | 204.51 | 9.62 |

^{a}Weighted by the flow volume.

^{b}Sum value.

To obtain the prior distribution of targeted parameters *k _{j}* in the model, a database of

*k*for NH

_{j}_{4}

^{+}, COD, TN, and TP was compiled from the literature. Specifically, the prior distribution of

*k*was determined as normal distribution (Liu

_{j}*et al.*2008; Chen

*et al.*2012), with a range of 0.005–0.70 d

^{−1}for NH

_{4}

^{+}(Alexander

*et al.*2000; Guo

*et al.*2008; Liu

*et al.*2008; Gao 2014; Zhang

*et al.*2015; Feng

*et al.*2017), 0.009–0.9 d

^{−1}for COD (Guo

*et al.*2008; Gao 2014; Zhang

*et al.*2015), 0.005–0.5 d

^{−1}for TN (Alexander

*et al.*2000; Feng

*et al.*2017), and 0.01–0.6 d

^{−1}for TP (Zhang

*et al.*2015; Feng

*et al.*2017).

The in-stream travel times (T or Ti) are determined through the length between the points or NPS location with the downstream section divided by the mean velocity (Chen *et al.* 2011). The length was measured in Google Earth, and the velocity was defined as normal distribution (Liu *et al.* 2008) with mean value measured from H1 and a standard deviation of 0.01–0.02 m/s.

The prior distribution of the volumetric flow rate (*p*) was determined according to the point source flow and monitored flow data for H1, in which the NPS discharge was defined as normal distribution with a mean of 15 m^{3}/s and the interval value of 10–30 m^{3}/s by the measured discharge from points and H1. The prior distributions of the key parameters are summarized in Table 3.

Symbol (unit) . | Parameter . | Variation range . | Prior distribution . |
---|---|---|---|

Cn1 (mg/L) | NPS concentration of NH_{4}^{+} | 0.09–8.14 | dnorm (2,1) |

Cn2 (mg/L) | NPS concentration of COD | 7.8–66.5 | dnorm (25,0.04) |

Cn3 (mg/L) | NPS concentration of TP | 2.09–12.4 | dnorm (4,1) |

Cn4 (mg/L) | NPS concentration of TN | 0.02–0.98 | dnorm (0.2,1) |

k1 (d^{−1}) | First-order attenuation coefficient of NH_{4}^{+} | 0.005–0.7 | dnorm (0.09,100) |

k2 (d^{−1}) | First-order attenuation coefficient of COD | 0.009–0.9 | dnorm (0.15,100) |

k3 (d^{−1}) | First-order attenuation coefficient of TN | 0.005–0.5 | dnorm (0.10,100) |

k4 (d^{−1}) | First-order attenuation coefficient of TP | 0.01–0.6 | dnorm (0.09,100) |

U.tau (s^{2}/m^{2}) | Reciprocal of velocity variance | 400–10,000 | dnorm (2500,0.01) |

Qn (m^{3}/s) | NPS discharge | 10–30 | dnorm (15,0.01) |

Symbol (unit) . | Parameter . | Variation range . | Prior distribution . |
---|---|---|---|

Cn1 (mg/L) | NPS concentration of NH_{4}^{+} | 0.09–8.14 | dnorm (2,1) |

Cn2 (mg/L) | NPS concentration of COD | 7.8–66.5 | dnorm (25,0.04) |

Cn3 (mg/L) | NPS concentration of TP | 2.09–12.4 | dnorm (4,1) |

Cn4 (mg/L) | NPS concentration of TN | 0.02–0.98 | dnorm (0.2,1) |

k1 (d^{−1}) | First-order attenuation coefficient of NH_{4}^{+} | 0.005–0.7 | dnorm (0.09,100) |

k2 (d^{−1}) | First-order attenuation coefficient of COD | 0.009–0.9 | dnorm (0.15,100) |

k3 (d^{−1}) | First-order attenuation coefficient of TN | 0.005–0.5 | dnorm (0.10,100) |

k4 (d^{−1}) | First-order attenuation coefficient of TP | 0.01–0.6 | dnorm (0.09,100) |

U.tau (s^{2}/m^{2}) | Reciprocal of velocity variance | 400–10,000 | dnorm (2500,0.01) |

Qn (m^{3}/s) | NPS discharge | 10–30 | dnorm (15,0.01) |

dnorm(*α*,*β*) denotes the norm distribution with the mean *α* and the reciprocal of variance *β.*

### Model calibration

The OpenBUGS software was used to develop and run this model. The observations of each pollutant from S2 were used to calibrate the model by adjusting the mean and variance of the prior distribution. The detailed data can be seen in Supplementary Material (Table S1) (available with the online version of this paper).

In order to be sure that the MCMC chain is a truly representative sample from the distribution, Gelman–Rubin Scale Reduction factor (SR) statistic (Gelman & Rubin 1992) and autocorrelation function (ACF) were used (Kruschke 2011). SR checks if the chains get stuck in unrepresentative regions of parameter space by the variability between MCMC chains relative to that within chains. If not, the SR should be around 1. ACF is to measure the correlation of the chain values with the values lag steps behind (or ahead). A highly auto-correlated chain is ‘clumpy’, which over-represents some values and under-represents other values. The way to mitigate the problem is to thin the chain until the ACF of lag (>0) steps close to zero. The determination coefficient (*R*^{2}) between the observed data and modelled values was calibrated, which should be close to 1.

It should be noted that the MCMC has a burn-in period due to the reduced influence of the initial value, and early simulation samples are often omitted by the Gelman–Rubin statistical diagnoses. In Figure 3, the Gelman–Rubin statistic for *k* chains shows that the model converges (SR close to 1) after several hundred iterations; therefore, the first 1,000 runs are discarded in our model to ensure that the chains converge. The chains were thinned to reduce their auto-correlation. The ACF for *k* chains with the thin step of 20 was shown in Figure 4, which shows that the chains are not ‘clumpy’ (ACF (lag > 0) close to 0).

The fit between the observed data and modelled values of the downstream section (S2) is calibrated by the deterministic coefficient, as shown in Figure 5. The results are acceptable when considering the complexities of pollutant transport and the measuring error. Therefore, the model can be used to assess NPS loads.

### NPS pollution loads using the Bayesian source apportionment mixing model

The posterior distributions of *k* varied significantly among the different pollution species, as described in Figure 6. The statistics of the 10-day *p _{n}* are presented in Table 4, including the mean value, median value, standard deviation (S.D.), Monte Carlo error (MC error), and two confidence levels (2.5% and 97.5%). As shown in Table 3, the NPS contributions of

*p*decrease as the water increases upstream when the Erhe floodgate opens, indicating that the model is reasonable.

_{n}P[_{n}i]
. | Mean . | S.D. . | MC error . | 2.50% . | Median . | 97.50% . |
---|---|---|---|---|---|---|

P[1] _{n} | 0.5857 | 0.05101 | 5.21 × 10^{−4} | 0.4989 | 0.5869 | 0.6679 |

P[2] _{n} | 0.4725 | 0.06489 | 6.76 × 10^{−4} | 0.3644 | 0.4727 | 0.5782 |

P[3] _{n} | 0.4006 | 0.1036 | 0.001163 | 0.2122 | 0.4046 | 0.5731 |

P[4] _{n} | 0.1911 | 0.05075 | 4.68 × 10^{−4} | 0.1007 | 0.192 | 0.2759 |

P[5] _{n} | 0.1525 | 0.043 | 4.12 × 10^{−4} | 0.08456 | 0.1494 | 0.2339 |

P[6] _{n} | 0.5949 | 0.04907 | 5.01 × 10^{−4} | 0.5166 | 0.593 | 0.6794 |

P[7] _{n} | 0.3339 | 0.09155 | 0.000952 | 0.1798 | 0.3318 | 0.4983 |

P[8] _{n} | 0.1562 | 0.0441 | 4.46 × 10^{−4} | 0.08713 | 0.1525 | 0.2409 |

P[9] _{n} | 0.1493 | 0.04265 | 4.08 × 10^{−4} | 0.0848 | 0.1452 | 0.2338 |

P[10] _{n} | 0.153 | 0.0437 | 4.58 × 10^{−4} | 0.08748 | 0.149 | 0.2413 |

P[_{n}i]
. | Mean . | S.D. . | MC error . | 2.50% . | Median . | 97.50% . |
---|---|---|---|---|---|---|

P[1] _{n} | 0.5857 | 0.05101 | 5.21 × 10^{−4} | 0.4989 | 0.5869 | 0.6679 |

P[2] _{n} | 0.4725 | 0.06489 | 6.76 × 10^{−4} | 0.3644 | 0.4727 | 0.5782 |

P[3] _{n} | 0.4006 | 0.1036 | 0.001163 | 0.2122 | 0.4046 | 0.5731 |

P[4] _{n} | 0.1911 | 0.05075 | 4.68 × 10^{−4} | 0.1007 | 0.192 | 0.2759 |

P[5] _{n} | 0.1525 | 0.043 | 4.12 × 10^{−4} | 0.08456 | 0.1494 | 0.2339 |

P[6] _{n} | 0.5949 | 0.04907 | 5.01 × 10^{−4} | 0.5166 | 0.593 | 0.6794 |

P[7] _{n} | 0.3339 | 0.09155 | 0.000952 | 0.1798 | 0.3318 | 0.4983 |

P[8] _{n} | 0.1562 | 0.0441 | 4.46 × 10^{−4} | 0.08713 | 0.1525 | 0.2409 |

P[9] _{n} | 0.1493 | 0.04265 | 4.08 × 10^{−4} | 0.0848 | 0.1452 | 0.2338 |

P[10] _{n} | 0.153 | 0.0437 | 4.58 × 10^{−4} | 0.08748 | 0.149 | 0.2413 |

*P _{n}*[

*j*] denotes the NPS volumetric flow rate proportion of the

*j*

_{t}_{h}monitored data set.

The posterior distributions of the NPS loads are described in Figure 7, including the mean value and the 50% and 95% posterior confidence intervals. The mean NPS load of NH_{4}^{+} was 20.47 g/s (0.259 g/m^{3} d), which was lower than the mean value of 0.544 g/m^{3} d in the study of the Hun-Taizi River system in China (Liu *et al.* 2008). The possible reason for this is that the load in the study by Liu includes point sources. The mean TN load is 39.81 g/s (3.44 × 10^{3} kg/d) and is in the range of 1.0 × 10^{3} – 36.4 × 10^{3} kg/d in the study of the ChangLe river system in China (Chen *et al.* 2012). This result may be due to our evaluation period being the non-flood season, while the period in Chen's study was annual. The NPS load is often higher in flood season than the non-flood season (Shen & Zhao 2010; Chen *et al.* 2011).

The posterior distributions of the NPS loads proportional to the nonpoint and point source loads (PLns) in the river are shown in Figure 8. The total loads of NH_{4}^{+}, COD, TN, and TP from NPSs in the region accounted for approximately 16.35–27.58%, 18.78–25.69%, 21.68–29.71%, and 42.11–52.09%, respectively. It is clear that the point source load is the main source of pollution in this channel for COD, NH_{4}^{+}, and TN, but TP accounts for approximately half of the total pollution. The first national pollution source census (MOEP *et al.* 2010) shows that the proportions of COD, TN, and TP in water environments from agricultural NPS pollution are 43.7%, 57.2%, and 67%, respectively. Obviously, it is reasonable that the NPS pollution ratios in our study, where point sources are the main source of pollution, are lower than that in China.

This method can be applied on rivers with high velocity, where the algae biomass is relatively low and the uptake amount of ammonium by algae was small. However, in rivers with algae blooming and eutrophication, algae uptaking and the benthic flux should be considered to better describe the attenuation process of ammonium, and corresponding likelihood function by the Bayesian method should be developed to cover the coupled attenuation process.

### Parameter sensitivity analysis

To analyse the degree to which an input parameter affects the model outputs, the sensitivity index of model parameters for the model outputs PLns is calculated by Equation (9), and the results are given in Table S1. It is clear that the NPS concentration of pollution has the greatest impact on corresponding model output PLns, and the NPS discharge is the second sensitive parameter. The velocity and the first-order attenuation have low sensitivities. All the sensitivity indices are below one, which means that the model output shrinks the uncertainty of the parameters.

So, the sensitivity ranking of the input parameters sorted by the amount of influence in our model is Cn, Qn, Uh1, *k*, and U.tau. In practice, the parameters are not included in routine monitoring data, so experience or historical data should be used to define a reasonable prior distribution in order to get a precise modelling result. Besides, the more detailed the data, less uncertain is the result.

## CONCLUSION

An improved inverse model based on the SAMM and coupled with the Bayesian approach for NPS load estimation was constructed. The model assumed that the pollutant concentrations in downstream samples were a mixture from various sources after attenuation, and the volumetric flow rate contributions were apportioned based on known information and pollution source emission characteristics. The model assessed NPS loads via a source analysis, which was carried out based on conventional and limited monitoring data without a sufficient number of discharge sites. The model not only considered the decrease in pollution transport but also addressed uncertainties in the estimation parameters, which are important for environmental risk management. This practical application shows that the model can serve as a simple and effective tool for researchers and managers to estimate NPS loads and uncertainties based on a small amount of water quality and hydrology data.

The Bayesian source apportionment mixing model (BSAMM) applied in this study roughly classified all pollution sources except for upstream and point source as non-point sources, which may include rural wastewater, agricultural wastewater, groundwater, and unmonitored industrial wastewater. If detailed pollution source data could be collected, more precise results regarding the source and volume of the pollution can be calculated.

Parameter sensitivity analysis shows the prior distribution of NPS concentration exert the greatest influence on modelling, and thereby sufficient empirical or historical information would facilitate a reasonable prior distribution and more precise result.

## ACKNOWLEDGEMENT

This work was supported by the Major Science and Technology Program for Water Pollution Control and Treatment in China (No. 2017ZX07602-003 and No. 2014ZX07204-005).

## REFERENCES

*The First National Pollution Source Census*. Available from: http://www.stats.gov.cn/tjsj/tjgb/qttjgb/qgqttjgb/201002/t20100211_30641.html (in Chinese)