An improved source apportionment mixing model combined with a Bayesian approach for nonpoint source pollution load estimation

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.

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.

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_n, and p_i. Then, the Bayesian approach is used to calibrate the unknown parameters.

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).

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². 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).

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₄⁺, 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.

To obtain the prior distribution of targeted parameters k_j in the model, a database of k_j for NH₄⁺, COD, TN, and TP was compiled from the literature. Specifically, the prior distribution of k_j was determined as normal distribution (Liu et al. 2008; Chen et al. 2012), with a range of 0.005–0.70 d⁻¹ for NH₄⁺ (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⁻¹ for COD (Guo et al. 2008; Gao 2014; Zhang et al. 2015), 0.005–0.5 d⁻¹ for TN (Alexander et al. 2000; Feng et al. 2017), and 0.01–0.6 d⁻¹ 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³/s and the interval value of 10–30 m³/s by the measured discharge from points and H1. The prior distributions of the key parameters are summarized in Table 3.

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²) 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.

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_n decrease as the water increases upstream when the Erhe floodgate opens, indicating that the model is reasonable.

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₄⁺ was 20.47 g/s (0.259 g/m³ d), which was lower than the mean value of 0.544 g/m³ 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³ kg/d) and is in the range of 1.0 × 10³ – 36.4 × 10³ 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₄⁺, 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₄⁺, 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.

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.

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.

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).