Combined Exceedance Probability Assessment of Water Quality Indicators Based on Multivariate Joint Probability Distribution in Urban Rivers

Abstract

Discharge and water quality are two important attributes of rivers, although the joint response relationship between discharge and multiple water quality indicators is not clear. In this paper, the joint probability distributions are established by copula functions to reveal the statistical characteristics and occurrence probability of different combinations of discharge and multiple water quality indicators. Based on the data of discharge, ammonia nitrogen content index (NH₄⁺) and permanganate index (COD_Mn) in the Xiaoqing River in Jinan, we first tested the joint change-point with the data from 1980–2016, before we focused on analyzing the data after the change-point and established the multivariate joint probability distributions. The results show that the Gaussian copula is more suitable for describing the joint distribution of discharge and water quality, while the year of 2005 is a joint change-point of water quantity and quality. Furthermore, it is more reasonable to use the trivariate joint probability distribution as compared to the bivariate distributions to reflect the exceedance probability of water quality combination events under different discharge conditions. The research results can provide technical support for the water quality management of urban rivers.

1. Introduction

Water quantity and quality are two important attributes of rivers. In recent years, the water quality of rivers is deteriorating with the development of society and economy, with water dispatching having become one of the main technologies of water quality improvement [1,2,3,4]. However, due to the disturbance of urban rivers from human activities, the relationship between water quantity and water quality has become more complicated. The relationship between a single water quality indicator and the water quantity can hardly reflect the real situation of river pollution [5]. Thus, the methods for improvement of the river water quality should be based on a correct understanding of the joint response relationship between water quantity and multiple water quality indicators.

Many scholars have simulated and predicted the river water quality from the aspects of mechanism models and non-mechanism models. The mechanism models are based on the hydrodynamic model and water quality equations in order to simulate the diffusion and attenuation process of pollutants [6,7,8]. However, due to the wide variety of influencing factors and the cognitive limitations of the real water quality, this makes the model construction and parameter determination difficult [9,10,11]. Therefore, this indicates that the simulation effect of the nature river is better than that of the urban river, which is affected by human activities. The non-mechanism model can avoid this type of problems to some extent. Sun et al. [12] used the gray model group method and the exponent smoothing method to construct a probabilistic combination prediction framework of water quality to forecast the individual indicator of chemical oxygen demand in a complex environment. After this, some new models for water quality prediction have been proposed, which integrates different non-mechanism models, such as the neural networks model, fuzzy model, and wavelet model [13,14,15]. In addition, there are also some joint models that combine hydrodynamics and water quality simulation tests [16,17]. However, most non-mechanism models require that the variables obey the same marginal distribution or that these variables need to be transformed to the same distribution, which not only leads to the loss of the original data information, but it also does not conform to the fact that variables, such as discharge and water quality, obey different marginal distributions. Meanwhile, because less emphasis is placed on the joint distribution between the variables, non-mechanism models are more suitable to use a simulation of the discharge and individual water quality indicator when compared to a simulation of the discharge and multiple water quality indicators. Copula can construct the joint probability distribution of multiple variables [18] and it does not require each variable to follow the same marginal distribution. Thus, it can effectively preserve the nature of the original data. Copula has the advantages of great adaptability and flexibility, so it has been widely used in the field of hydrology, especially in determining flood risks [19,20,21,22] and drought characteristics [23,24,25,26,27,28]. However, copula has been used less frequently in the joint analysis of discharge and water quality, with most of the uses focusing on bivariate joint distribution. As such, there is a distinct lack of studies that use copula for a multivariable joint distribution.

Therefore, the purpose of this paper is to use copulas to establish the multivariate joint probability distribution of the discharge and water quality indicators, before analyzing the exceedance probability of multiple water quality indicators under different conditions of discharge.

2. Methods and Materials

2.1. Methodology

2.1.1. Establishment of Copula Function

Copula is a connection function, which is defined in R with marginal distribution functions $F_1,F_2,\cdots ,F_n$ for n-dimensional random variables $X_1,X_2,\cdots ,X_n$ [29]. If H is the n-dimensional distribution function of n-dimensional random variables, for any $x_i\in R(i=1,2,\cdots ,n)$, there exists the following copula function:

$$\begin{array}{cc}\hfill H_{x_1,x_2,\dots ,x_d}(x_1,x_2,\dots ,x_d)& =P(X_1\le x_1,X_2\le x_2,\cdots ,X_n\le x_n)\hfill \\ & =C(F_1(x_1)，F_2(x_2)，\dots ，F_n(x_n))\hfill \end{array}$$

(1)

where n is the number of variables; C is the associated dependence copula function; and, $F_i(x_i)$ is the marginal distribution functions.

Archimedean copulas [30] and Elliptical copulas [31] are commonly used in hydrology. The typical Archimedean copulas include Clayton, Gumbel-Hougaard (G-H), and Frank copula, which require the variables to fall within a symmetric structure. Elliptical copulas are neither limited by the correlation nor restricted by the symmetric structure between variables. However, the multiple hydrological variables may be asymmetric, with independent, positive, or negative relationships occurring among the variables. Thus, it is more reasonable to use the Elliptical copulas to describe the joint distribution of variables when compared to the Archimedean copulas. The Gaussian copula in Elliptical copula is expressed as:

$$\begin{array}{cc}\hfill C(u_1,u_2,\dots ,u_d；{\displaystyle \sum )}& ={\mathsf{\Phi }}_{{\displaystyle \sum }}({\mathsf{\Phi }}^{-1}(u_1),{\mathsf{\Phi }}^{-1}(u_2),\cdots ,{\mathsf{\Phi }}^{-1}(u_d))\hfill \\ & ={\displaystyle {\int }_{-\infty }^{{\mathsf{\Phi }}^{-1}(u_1)}\cdots }{\displaystyle {\int }_{-\infty }^{{\mathsf{\Phi }}^{-1}(u_d)}\frac{1}{{(2\pi )}^{\frac{d}{2}}{\displaystyle \sum {}^{\frac{1}{2}}}}}\exp (-\frac{1}{2}{w}^T{\displaystyle \sum {}^{-1}}w)dw\hfill \end{array}$$

(2)

where ${\mathsf{\Phi }}^{-1}(.)$ is the inverse function of the standard normal distribution; ${\mathsf{\Phi }}_{{\displaystyle \sum }}({\mathsf{\Phi }}^{-1}(u_1),{\mathsf{\Phi }}^{-1}(u_2),\cdots ,{\mathsf{\Phi }}^{-1}(u_d))$ is the multivariate normal distribution; ${\displaystyle \sum =\left[\begin{array}{ccc}1& \cdots & {\rho }_{1d}\\ ⋮& \ddots & ⋮\\ {\rho }_{d1}& \cdots & 1\end{array}\right]}$, ${\rho }_{ij}=\{\begin{array}{c}1\begin{array}{cc},\begin{array}{cc}& \end{array}& \end{array}i=j\\ {\rho }_{ji},\begin{array}{cc}& \end{array}i=j\end{array}$, $-1<{\rho }_{ij}<1$; d is the dimensions of the random variable; and, $\mathbf{w}={[w_1,w_2\cdots w_d]}^T$.

2.1.2. Test of Joint Change-Point

The purpose of testing the change-point of the multivariable dependence is to make the data more representative in different periods. At the same time, we can utilize the collected data as much as possible to improve the accuracy of the simulation.

Assuming that there is one change-point in the copula function constructed based on the observed data $(x_1,y_1,z_1),\cdots (x_n,y_n,z_n)$, the original hypothesis and the alternative hypothesis can be, respectively, assumed as:

$$H_0:{\theta }_1={\theta }_2=\cdots ={\theta }_n{,}_{}^{}\hspace{1em}{H}_1:{\theta }_1=\cdots {\theta }_n{}_{*}\ne {\theta }_{n^{*}+1}=\cdots {\theta }_n$$

If the original hypothesis is rejected, then $k^{*}$ is the moment of change. If $k^{*}=k$, the log likelihood ratio statistic of copula can be established based on the maximum likelihood estimation method, which is shown in Equation (4):

$$\begin{array}{cc}\hfill -2\ln ({\mathsf{\Lambda }}_k)& =2[{\displaystyle \sum _{i=1}^k{C}_{uvw}({\theta }_k;U(x_i),V(y_i),W(z_i))}+{\displaystyle \sum _{i=k+1}^n{C}_{uvw}({\theta }_{k_{}^{*}};U(x_i),V(y_i),W(z_i))}\hfill \\ \hfill & \hspace{1em}+{\displaystyle \sum _{i=1}^n{C}_{uvw}({\theta }_n;U(x_i),V(y_i),W(z_i))}]\hfill \end{array}$$

(3)

where ${\theta }_k,{\theta }_{k_{}^{*}},{\theta }_n$ are the maximum likelihood estimations of the corresponding data; and, $C_{uvw}(U(x),V(y),W(z))=\frac{\partial }{\partial u}\frac{\partial }{\partial v}\frac{\partial }{\partial w}C(u,v,w)$.

The statistic $Z_n=\underset{1\le k<n}{\max }(-2\ln ({\mathsf{\Lambda }}_k))$ is constructed. If $k^{*}$ is unknown, the original hypothesis will be rejected when $Z_n$ is greater than the critical value, which essentially demonstrates that the copula function has change-points [32].

2.1.3. Exceedance Probability of Water Quality

In order to quantitatively analyze the response relationship between the discharge and water quality, according to the full probability formula and the copula functions [33,34], we can calculate the probability that the water quality indicators exceed specific values. However, discharge needs to be less than a specific value, which addresses the exceedance probability of water quality. Taking the joint probability distribution of two and three dimensions as an example, the expressions are as follows:

$$\begin{array}{ccc}\hfill P_{U_{1}U{}_2}& =P(U_1\le u_1,U_2\ge u_2)\hfill & =F_{U_1}(u_1)-F_{U_{1}U{}_2}(u_1,u_2)\hfill \\ \hfill & \hfill & =F_{U_1}(u_1)-C(u_1,u_2)\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill P_{U_1{U}_2{U}_3}& =P(U_1\le u_1,U_2\ge u_2,U_3\ge u_3)\hfill \\ \hfill & \hspace{1em}=F_{U_1}(u_1)-F_{U_{1}U{}_2}(u_1,u_2)-F_{U_{1}U3}(u_1,u_3)+F_{U_{1}U{}_{2}U{}_3}(u_1,u_2,u_3)\hfill \\ \hfill & \hspace{1em}=F_{U_1}(u_1)-C(u_1,u_2)-C(u_1,u_3)+C_{}(u_1,u_2,u_3)\hfill \end{array}$$

(5)

where $P_{U_1{U}_2}$ is the value of the bivariate joint probability distribution; $P_{U_1{U}_2{U}_3}$ is the value of the trivariate joint probability distribution; $U_1$ is the discharge; $U_2$ and $U_3$ are the water quality indicators; $F_{U_i}(u_i)$ is the marginal distribution function of $U_1$; and, $C(u_1,\cdots ,u_i)$ is the copula function of pair $(U_1,\cdots ,U_i)$.

2.2. Study Area and Data

Xiaoqing River, which is located in Jinan, is a manually excavated city river. Furthermore, it is the main flood discharge channel in Jinan. The Huangtaiqiao hydrologic station is located in the main stream of the river, which has monitored the discharge from the river since 1953 with complete water quality monitoring having begun in 1970. The distance between the Muli Gate to Huangtaiqiao hydrologic station is 21.0 km, with a drainage area of 321.0 km² and an average annual discharge of 10.4 m³/s. The river basin is located in the main urban area of Jinan. Due to the rapid development of urbanization, the water pollution in this area is serious. In order to build the ecological civilization, there has been water diversion from the Eastern Route of South-to-North Water Transfer Project to improve the water quality. However, the combined effect of urbanization and water diversion result in the response relationship between the water quantity and quality in Xiaoqing River has been severely disturbed, which has seriously affected its water quality simulation. Therefore, the statistical analysis of the occurrence probability of the combined events of water quantity and water quality can provide a basis for the scientific dispatch. The drainage map of Xiaoqing River in Jinan can be seen in Figure 1.

This study selected the data of discharge and water quality of Huangtaiqiao station during 1980–2016, with the data of water quality having been recorded every other month and the discharge having been recorded at the same time when the water quality samples were taken (using the main pollutants, nitrogen content index (NH₄⁺), which includes the ionised ammonia nitrogen and the free ammonia nitrogen, and permanganate index (COD_Mn) as the calculation indicators). To remove some unreasonable data, 161 sets of water quantity and water quality data were applied in this paper. The sample data series of NH₄⁺, COD_Mn, and discharge of the Huangtai Station from 1980 to 2016 are shown in Figure 2, and the characteristics and trend results of the generalized basic data about the water quality are given in

Table 1. The trend results of NH₄⁺, COD_Mn, and discharge.

	NH4+ (mg/L)	CODMn (mg/L)	Discharge (m3/s)	Water Temperature (°C)	pH
Minimum	0.43	2.60	3.99	3	5.50
Maximum	18.50	48.10	36.10	31	11.60
Mean	7.23	16.86	12.43	18.94	7.51
Trend	↓	↓ +	↑ +	↑ −	↑

Note: ↓ means a decreasing trend; ↑ means a raising trend; ⁺ means significant trend; ⁻ means non-significant trend.

. The trends of the generalized basic data are tested by Mann-Kendall (M-K). It can be seen that the concentrations of NH₄⁺ and COD_Mn are basically larger than the grade III value (Environmental Quality Standards for Surface Water (GB3838-2002)). Moreover, the COD_Mn has a significant decreasing trend, and the discharge has a significant increasing trend, while the trend of NH₄⁺ descending was not obvious.

In order to better reflect the relationship between the discharge and the multiple water quality indicators in the Xiaoqing River Basin, which is a rapidly developing area, we carried out a test for the change-point of the long sequence data (1980–2016), before focusing on the analysis of the data after the change-point.

3. Results and Discussion

3.1. Marginal Distributions and Correlation Coefficient of Q, NH₄⁺ and COD_Mn

The marginal distributions of the Q, NH₄⁺, and COD_Mn of Xiaoqing River are established by the lognormal, gamma, and Pearson-III distributions to search for the proper cumulative distribution function (CDF). The proper CDF of the parameters was chosen using two goodness fit tests: the Kolmogorov–Smirov (K–S) test and root mean squared error (RMSE) test. The results of the proper marginal distribution of each variable are shown in

Table 2. The parameters of marginal distributions of Q, NH₄⁺, and COD_Mn and results of Kolmogorov–Smirov (K–S) and root mean squared error (RMSE).

	Lognormal			Gamma			Pearson III
	$(\mathit{\mu },\mathit{\sigma }$)	KS	RMSE	$(\mathit{\mu },\mathit{\sigma }$)	KS	RMSE	$(\mathit{\mu },\mathit{\sigma },{\mathit{\alpha }}_0$)	KS	RMSE
Q	2.398, 0.494	0.057	0.028	4.255, 2.921	0.051	0.021	0.553, 23.406, −0.508	0.431	0.228
NH4+	1.738, 0.751	0.152	0.026	2.238, 3.230	0.063	0.022	2.367, 3.149, −0.222	0.066	0.028
CODMn	2.588, 0.690	0.072	0.053	2.266, 7.440	0.121	0.067	1.470, 10.291, 1.737	0.097	0.055

Note: The number in bold means that RSME is the minimum and the K–S is 0.10172 according to the 5% significance level when n = 161.

. From Table 2, we can see that the marginal distribution of Q and NH₄⁺ sequence is more consistent with a Gamma distribution, while COD_Mn is more consistent with a lognormal distribution.

The Kendall correlation coefficient $\tau$ is used to analyze the correlation of Q, NH₄⁺, and COD_Mn, respectively, which is shown in

Table 3. The correlation coefficient $\tau$ of Q, NH₄⁺, and COD_Mn.

Pair	Q, NH4+	Q, CODMn	NH4+, CODMn
$\tau$	−0.2712	−0.3992	0.1933

. We can see that the pair of (Q, NH₄⁺) and the pair of (Q, COD_Mn) are negatively correlated, while the pair (NH₄⁺, COD_Mn) has a positive correlation. This provides the conditions for the application of Gaussian Copula.

3.2. Establishment of Joint Probability Distribution

3.2.1. Establishment of Multivariate Joint Probability Distribution of Q, NH₄⁺, and COD_Mn with the Data from 1980–2016

The Frank, GH, and Gaussian copulas are used to construct the bivariate and trivariate joint probability distributions with the data of water quantity and water quality of Xiaoqing River in Jinan that were obtained during 1980–2016. The proper distribution was chosen using K–S and the RMSE tests, with the parameters of copulas and the results of the goodness fit being shown in

Table 4. The parameters of copulas of Q, NH₄⁺, and COD_Mn, and results of K–S and RMSE.

	Pairs	Parameters	K–S	RMSE
Frank	Q, NH4+	$\theta =-0.6114$	0.1032	0.0541
	Q, CODMn	$\theta =-0.5890$	0.0867	0.0353
	Q, NH4+, CODMn	$\theta =-0.9315$	0.1078	0.0331
GH	Q, NH4+	$\theta =0.7866$	0.0842	0.0251
	Q, CODMn	$\theta =0.7147$	0.0892	0.0337
	Q, NH4+, CODMn	$\theta =0.9236$	0.1003	0.0343
Gaussian	Q, NH4+	$\rho =-0.4036$	0.0812	0.0251
	Q, CODMn	$\rho =-0.5319$	0.0752	0.0323
	Q, NH4+, CODMn	$$\begin{gathered} {\rho }_{12}=-0.4036 \\ {\rho }_{13}=-0.5319 \\ {\rho }_{23}=0.2645 \end{gathered}$$	0.0906	0.0222

Note: ${\rho }_{12}$ is the correlation coefficient of (Q, NH₄⁺); ${\rho }_{13}$ is the correlation coefficient of (Q, COD_Mn); and ${\rho }_{23}$ is the correlation coefficient of (NH₄⁺, COD_Mn). The K–S is 0.1626 according to the 5% significance level when n = 70.

. Figure 3 compares the sample empirical distribution with the joint distribution of the pairs of (Q, NH₄⁺), (Q, COD_Mn), and (Q, NH₄⁺, COD_Mn). From these figures, we can see that the correlation coefficients R² of Gaussian copula are higher. This illustrates that the Gaussian copula is more suitable in establishing both the bivariate joint probability distributions and the trivariate joint probability distribution.

3.2.2. Test of Joint Change-Point of Q, NH₄⁺, and COD_Mn

According to the optimal fitting copula functions of Q, NH₄⁺, and COD_Mn using the data from 1980–2016, the joint change-points can be identified. When the number of data is 161, the critical value of $Z_n$ is 13.35 [35]. The results are shown in Figure 4, which show that $Z_n$ is greater than its critical value in No. 94 of the data, which is July 2005. Therefore, this is the change-point of the combination of (Q, NH₄⁺, COD_Mn). Furthermore, the result also matches the actual situation: the continuous drought for several years before 2005 in Jinan resulted in the management department transferring about 1 × 10⁸ m³ of water into the Xiaoqing River, which might cause a change in the relationship between the water quantity and quality.

This indicates that the dependency relationship of the data was inconsistent between 1980 and 2016, so a change-point test was needed to identify the data sequences that are consistent with the current situation.

3.2.3. Establishment of Multivariate Joint Probability Distribution after the Joint Change-Point

The data of Q, NH₄⁺, and COD_Mn after the joint change-point is more representative of the relationship between water quantity and quality in recent years. Therefore, we re-established the multivariate joint probability distribution using the data after the change-point. The parameters of the copulas and the results of the goodness fit are shown in

Table 5. The parameters of copulas of Q, NH₄⁺, and COD_Mn and the results of K–S and RMSE.

	Pairs	Parameters	K–S	RMSE
Frank	Q, NH4+	$\theta {=}_{}^{}-0.5869$	0.1667	0.0654
	Q, CODMn	$\theta {=}_{}^{}-0.5971$	0.1546	0.0621
	Q, NH4+, CODMn	$\theta {=}_{}^{}-0.2998$	0.1816	0.0651
GH	Q, NH4+	$\theta {=}_{}^{}0.8295$	0.2136	0.0711
	Q, CODMn	$\theta {=}_{}0.9509$	0.2095	0.0604
	Q, NH4+, CODMn	$\theta {=}_{}2.8152$	0.1791	0.0743
Gaussian	Q, NH4+	$\rho =-0.4073$	0.1431	0.0462
	Q, CODMn	$\rho =-0.4037$	0.1315	0.0493
	Q, NH4+, CODMn	$$\begin{gathered} {\rho }_{12}=-0.4073 \\ {\rho }_{13}=-0.4037 \\ {\rho }_{23}=0.4151 \end{gathered}$$	0.1163	0.0524

. From Table 5, we can see that the joint probability distributions are still consistent with the Gaussian copula function. Meanwhile, their corresponding contour plots and contour surface are presented in Figure 5.

3.3. Analysis of Multivariate Joint Probability Distribution of Q, NH₄⁺ and COD_Mn after the Joint Change-Point

According to the real situation of the river water quality and the existing water quality standards “Environmental Quality Standards for Surface Water” (GB3838-2002), we conducted a cluster analysis of the water quality to obtain the new classes of water quality, which are shown in

Table 6. Classes of water quality for NH₄⁺ and COD_Mn.

	Cluster Center				Boundary Values
	Center 1	Center 2	Center 3	Center 4	Class 1	Class 2	Class 3	Class 4	Class 5
NH4+	4.2	6.3	8.5	11.2	<4.2	4.2~6.3	6.3~8.5	8.5~11.2	>11.2
CODMn	8.0	12.0	24.0	35.0	<8.0	8.0~12.0	12.0~24.0	24.0~35.0	>35.0

.

Using the constructed joint probability distribution of Q, NH₄⁺, and COD_Mn, the probability of the water quality exceeding a specific value under the different conditions of discharge can be analyzed by Equations (3) and (4). $P_{U_1{U}_2}$ and $P_{U_1{U}_3}$ are the bivariate joint probability distributions of the pair of (Q, NH₄⁺) and the pair of (Q, COD_Mn), respectively, which represent the exceedance probability of NH₄⁺ or COD_Mn in different discharge situations. $P_{U_{1}U{}_{2}U{}_3}$ is the trivariate joint probability distribution of the pair of (Q, NH₄⁺, COD_Mn), which represents the combined exceedance probability of (NH₄⁺, COD_Mn) in different discharge situations. The contour plots of $P_{U_1{U}_2}$ and $P_{U_1{U}_3}$ are given in Figure 6a,b, respectively, while the contour surface of $P_{U_{1}U{}_{2}U{}_3}$ is shown in Figure 6c.

From Figure 6a,b, we can see that when $P_{U_1{U}_2}$ or $P_{U_1{U}_3}$ is in the range of 0.1–0.4, the contour plot is relatively sparse, especially in the upper right panel of the figures. This indicates that when the concentration of NH₄⁺ or COD_Mn is high, the joint probability is insensitive to a change in discharge. For the combination of NH₄⁺ and COD_Mn, Figure 6c shows that, when $P_{U_{1}U{}_{2}U{}_3}$ is between 0.1–0.3, the contour surface is relatively sparse. This is mainly due to the large amount of pollutants in the water body of Xiaoqing River at this moment. At the same time, the tributaries also carry pollutants so the joint probability is insensitive to a change in the discharge.

The typical values of the joint probability of Q, NH₄⁺, and COD_Mn are given in

Table 7. The typical values of exceedance probability of NH₄⁺.

NH4+	Q(m3/s)
(mg/L)	3	6	12	24	36
4.2	0.0133	0.0653	0.2589	0.4079	0.4177
6.3	0.0120	0.0552	0.1985	0.2916	0.2964
8.5	---	0.0426	0.1379	0.1900	0.1921
11.2	---	0.0287	0.0827	0.1070	0.1078

Note: “---” denotes impossible events.

,

Table 8. The typical values of exceedance probability of COD_Mn.

CODMn	Q(m3/s)
(mg/L)	3	6	12	24	36
8.0	0.0135	0.0645	0.2638	0.4147	0.4231
12.0	0.0134	0.0587	0.2135	0.3043	0.3077
24.0	---	0.0341	0.0896	0.1077	0.1080
35.0	---	0.0184	0.0395	0.0443	0.0443

Note: “---” denotes impossible events.

and

Table 9. The typical values of combined exceedance probability of NH₄⁺ and COD_Mn.

NH4+	CODMn	Q(m3/s)
(mg/L)	(mg/L)	3	6	12	24	36
4.2	8	0.0138	0.1056	0.3897	0.5596	0.5658
	12	0.0128	0.0964	0.3194	0.4244	0.4270
	24	0.0096	0.0566	0.1387	0.1607	0.1609
	35	0.0064	0.0308	0.0625	0.0685	0.0685
6.3	8	0.0120	0.0894	0.3011	0.4094	0.4125
	12	0.0116	0.0819	0.2493	0.3173	0.3187
	24	0.0088	0.0487	0.1116	0.1263	0.1264
	35	0.0059	0.0268	0.0514	0.0554	0.0554
8.5	8	0.0101	0.0691	0.2108	0.2725	0.2740
	12	0.0098	0.0635	0.1766	0.2159	0.2165
	24	---	0.0384	0.0818	0.0905	0.0906
	35	---	0.0214	0.0386	0.0411	0.0411
11.2	8	---	0.0467	0.1275	0.1569	0.1575
	12	---	0.0432	0.1082	0.1272	0.1274
	24	---	---	0.0522	0.0565	0.0566
	35	---	---	0.0254	0.0266	0.0266

Note: “---” denotes impossible events.

, which describe the probability of the water quality indicators exceeding the specific values with different combinations of discharge and water quality indicators. For example, under the condition that only NH₄⁺ is considered, $P_{U_1{U}_2}=P(U_1\le 12,U_2\ge 6.3)=0.1985$. Under the condition that only COD_Mn is considered, $P_{U_1{U}_3}=P(U_1\le 12,U_3\ge 12.0)=0.2135$. Finally, under the condition that NH₄⁺ and COD_Mn are synergistically considered, $P_{U_{1}U{}_{2}U{}_3}=P(U_1\le 12,U_2\ge 6.3,U_3\ge 12.0)=0.2493$. The results show that the trivariate joint probability of the combination of (Q, NH₄⁺, COD_Mn) is larger than the bivariate joint probability of the pairs of (Q, NH₄⁺) and (Q, COD_Mn). This indicates that the probability of exceeding the standards of multiple indicators is greater than that of exceeding individual indicators, which is more consistent with the actual situation. Therefore, only considering the discharge and individual water quality indicator cannot reasonably reflect the relationship between the discharge and water quality. The trivariate joint probability distribution considering Q, NH₄⁺, and COD_Mn can more comprehensively express the characteristics of discharge and multiple water quality indicators.

3.4. Discussion

There is not always a positive correlation between the discharge and water quality indicators, which means that the correlation between variables is asymmetric. In this paper, although Frank and G-H Copulas, which are symmetric copula functions, can fit the probability distribution between the discharge and water quality, the fitting accuracy, and the simulation of higher dimensions are not sufficient when compared to the asymmetric Gaussian copulas. Based on the statistical analysis of the historical data of discharge and water quality, this paper provides a method to quantitatively evaluate the exceedance probability of water quality by considering the joint probability of Q, NH₄⁺, and COD_Mn. When compared with the previous non-mechanical methods, this method adopts a multi-dimensional joint distribution that is based on a copula to more accurately describe the response relationship between water quality and discharge under the condition of the coexistence of multiple pollutants in rivers.

The trivariate joint probability distribution of the combination of (Q, NH₄⁺, COD_Mn) can reflect the joint probability for the different combinations of Q, NH₄⁺, and COD_Mn. Thus, when the combination events (NH₄⁺, COD_Mn) exceed their given design values in different discharge situations, the exceedance probability can be calculated. Generally, regardless of whether the pollutant is NH₄⁺, COD_Mn, or a combination of NH₄⁺ and COD_Mn, the joint probability tends to increase with an increase in the discharge and a decrease in the concentration of pollutants discharged into rivers. Thus, the exceedance probability of the water quality can be reduced by an increase in discharge (Q) and a decrease in concentration of pollutants (NH₄⁺ and COD_Mn). However, with an increase in discharge, the counter lines or the counter surface tend to level, which means that there is a limit imposed on the improvement of water quality. Therefore, further research is needed to determine the limit value. In addition, we suggest that the administration should strengthen the monitoring of water quality at each entry section of Xiaoqing River so that we can obtain a more comprehensive relationship between the discharge and water quality.

4. Conclusions

In this paper, the joint probability distributions of discharge and the main pollutants of NH₄⁺ and COD_Mn of Xiaoqing River in Jinan were constructed by copulas, before the exceedance probability of water quality combination events under different discharge conditions was analyzed. The specific conclusions are as follows:

(1)

The Gaussian copula is more suitable for describing the multivariate joint probability distribution of discharge and water quality. As the relationship between the discharge and water quality indicators is not always positive, the Gaussian copula is more suitable than the Archimedean copulas in the simulation of trivariate or the above joint distribution.

(2)

Based on the copula, the joint change-point can be identified. For urban rivers, the dependence of water quantity and quality is often affected by human activities, which leads to the emergence of change-point. Therefore, it is necessary to identify the mutation points.

(3)

The trivariate joint probability distribution of the combination of (Q, NH₄⁺, COD_Mn) is more suitable for estimating the various exceedance probability of the water quality effectively under different discharge conditions. Thus, when the combination events (NH₄⁺, COD_Mn) exceed their given values in the specific discharge situation, the exceedance probability can be calculated.