Combined effect of climate change and human activities on the hydrological processes of an urban river: case study of the Xiaoqing River, China

ABSTRACT

The aim of this study is to reveal statistical characteristics and exceedance probability of discharge under the combined effect of climate change and human activities. The study is conducted in the Xiaoqing River in Jinan, China, based on data of discharge, land-use types and precipitation from the period 1970–2016. A multivariate joint probability distribution of the data is established to test the univariable, bivariable and trivariable change points. These are then used to calculate and analyse the risk probability of discharge exceeding the specific values under different conditions of precipitation and land-use type. The results show that the change point calculated by trivariate joint distribution can reduce the disturbance of the change point obtained with the univariable or bivariable approach and reflect the changes of various factors in the hydrological processes more objectively. When the land-use type is taken into consideration, the trivariate distribution can reflect the variation of hydrological processes more reasonably.

KEYWORDS:

• human activities
• copulas
• multivariate joint change point
• joint probability distribution
• hydrological processes
• urban rivers

Editor A. CastellarinGuest editor Y. Chen

1 Introduction

Human activities and climate change are two major factors affecting the hydrological processes of rivers, especially in rapidly urbanizing regions where hydrological processes are being seriously disturbed. Many rivers in urban regions have experienced an increasing frequency of floods over the past few years, which have led to ecological problems and affected livelihoods (Lian et al. 2012, Sun et al. 2012, Suriya and Mudgal 2012, Chen et al. 2015, Mcgrane 2015, Zhao et al. 2016). Thus, more and more attention has been paid to changes in river hydrological processes (e.g. Jay and Naik 2011, Zhang et al. 2015), especially to excessive discharge. Many researchers have employed distributed hydrological models to investigate the changes in hydrological processes (e.g. Paturel et al. 2003, Piniewski et al. 2013, Khoi and Suetsugi 2014). However, compared with natural rivers, the combined effect of climate change and human activities is more evident in urban regions. Moreover, with the intensification of urbanization, it not only has a direct impact on the hydrological processes of rivers, but also affects the urban climate where the rivers are located (Shepherd 2006, Chen et al. 2016, Cristiano et al. 2017). This then has an indirect effect on the hydrological processes of urban rivers in particular. With more interdependent factors affecting the hydrological processes of urban rivers, the calibration of parameters and the construction of the models are complex (Huang et al. 2015). As an alternative, non-mechanism models based on probability and statistics may be employed to avoid such problems to some extent. The elasticity method (Zhao et al. 2015, Li et al. 2016), the Budyko hypothesis and support vector machine (SVM) (Huang et al. 2016), as well as artificial neural networks (Guimaraes Santos and Silva 2014, Zemzami and Benaabidate 2015) are common non-mechanism models based on probability and statistics. These methods can be used to predict and simulate changes in hydrological processes by calculating and analysing the probability of hydrological events. However, most of them are based on the assumption that climate change and the human activities are independent factors (Peng et al. 2014, Xia et al. 2014, Ashofteh et al. 2016). Both hydrological models and non-mechanism models place emphasis on evaluating the contribution of either climate change or human activities to changes in hydrological processes, but it is relatively rare for the analysis to be of the joint effect of these two factors. Models ignore the combined effect of climate change and human activities (Liu et al. 2014), which means that their results cannot clearly reveal the change law of hydrological processes of rivers subject to a greatly changed environment. This suggests that a univariate or bivariate joint distribution of, for example, precipitation and discharge can no longer accurately reflect the changes in hydrological processes of urban rivers. Thus, to study highly urbanized regions, it is necessary to construct a multivariate joint probability distribution to reveal the occurrence probability of discharge under different conditions of precipitation, with the third variable of the land-use type.

Traditional joint probability distributions require all variables to have the same marginal distribution. However, in general, different variables may be subject to different marginal distributions, which leads to the information loss of data in the process of transformation (Zhang and Singh 2007). The copula probability distribution is one function that can overcome the drawbacks, and can connect with different univariate marginal distributions without changing the dependence structure of variables (Zegpi and Fernandez 2010, Dai et al. 2014). Thus, one can construct a joint probability distribution of any univariate marginal distribution, while the information included in initial data will be retained in the transformation process. Due to these advantages, the copula probability distribution has been applied in various hydrological fields, such as analysis of the probability of flood risk (e.g. Kao and Chang 2012, Bezak et al. 2014, Li and Zheng 2016, Debele et al. 2017), study of the drought characteristics of watersheds (e.g. Hao and Aghakouchak 2013, Yusof et al. 2013, Bazrafshan et al. 2015) and tests of the joint change point of precipitation and discharge (Guo et al. 2015, Lv et al. 2015). These studies are based on the two-dimensional copula. The land-use type has not been directly incorporated into the copula, and the application of three-dimensional copula probability distribution to analyse the combined effect of climate change and human activities on the hydrological processes of urban rivers has not been done to the authors’ knowledge.

Therefore, this paper introduces a novel approach to identify a joint change point of hydrological processes, integrating precipitation and land-use types, which are closely related to urbanization and river discharge. To this end, the exceedance probability of discharge after the joint change point under different combinations of precipitation and land-use type is calculated, and then the relationship between precipitation, land-use type and discharge is systematically analysed. The method can serve as a guide for urban managers to understand the change law of hydrological processes under different combined situations of precipitation and urbanization. The objectives of this paper are: (1) to construct a multivariate joint distribution of the main influence indexes (precipitation and land-use type) and the response index (discharge) and (2) to calculate the hydrological change points and the exceedance probability of discharge under the joint effect of precipitation and land-use type to evaluate the hydrological processes of an urban river in China. In the next section, the general characteristics of the study area and the collected data are presented, and the methods are introduced in Section 3. The results and discussion are given in Section 4, and the conclusions, and future research directions are summarized in Section 5.

2 Study area and data

The city of Jinan is located in the midwestern part of Shandong Province, eastern China (36°40′N, 116°57′E), and is in a monsoon climate area. Jinan is the capital city of Shandong Province with more than 8.8 million people and a GDP of US$86 billion. The Xiaoqing River is the main flood discharge channel that rises in Muli Gate; it spans five cities (Jinan, Binzhou, Zibo, Dongying and Weifang, in order from upstream to downstream) in Shandong Province and then discharges into the Bohai Sea. The study focuses on the 21 km reach between the Muli Gate in Jinan and the Huangtai hydrological station, which has a drainage area of 321.0 km². The Xiaoqing River basin covers the main urban area of Jinan. As the main flood discharge channel for Jinan, its water sources include precipitation in the flood season, groundwater and some treated wastewater. In recent years, with the rapid economic and social development of Jinan, the urban area has increased rapidly. This has resulted in dramatic changes in the underlying surface conditions in the study area, and the urban rain island effect has also been significantly enhanced, resulting in an increase in precipitation and significant changes in the hydrological process of the Xiaoqing River (Huang and He 2011).

The distribution of rivers, precipitation stations and hydrological stations in the study area is shown in Fig. 1. Figure 2 shows the average annual discharge for the period 1970–2016 monitored by the Huangtai hydrological station and the average annual precipitation monitored by seven precipitation stations. In addition, a total of 13 remote sensing images were collected from Landsat 3–8,¹ all of them in May to October. Details of the remote sensing images are shown in Table 1.

Table 1. Basic information of remote sensing images.

Number	Year	Data source	Acquisition time* (mm.dd)	Mean cloudiness (%)	Centre coordinates
1	1975	Landsat3	08.21	0	115.96E, 36.06N
2	1981	Landsat3	09.06	0	116.39E, 37.31N
3	1984	Landsat5	06.19	2.01	116.97E, 36.05N
4	1992	Landsat5	10.18	0	116.88E, 36.06N
5	1995	Landsat5	06.21	5.87	116.88E, 36.06N
6	2002	Landsat5	07.10	0	116.86E, 36.06N
7	2006	Landsat5	05.02	0.11	116.96E, 36.05N
8	2009	Landsat5	08.30	0.15	116.95E, 36.05N
9	2011	Landsat5	06.01	0.2	116.94E, 36.05N
10	2012	Landsat8	06.11	0	116.96E, 36.04N
11	2013	Landsat8	09.26	0.9	116.96E, 36.04N
12	2014	Landsat8	07.11	2.23	116.96E, 36.04N
13	2015	Landsat8	06.12	2.84	116.96E, 36.04N

Figure 1. Map of the Xiaoqing River Basin above Huangtai Station with the location of rainfall and hydrological stations.

Figure 2. Variation of the annual average (a) precipitation and (b) discharge of the Xiaoqing River basin above the Huangtai hydrological station from 1970 to 2016.

3 Research method

3.1 Extraction of land-use types

It is difficult to obtain a long data series of land-use type in most areas of China. At present, the commonly used method is to interpret remote sensing images of the study area. In this study, a method of combining ENVI (Environment for Visualising Images v5.5) software for the analysis and visualization of scientific data and imagery) and a CA (Cellular Automata)-Markov model is proposed to deduce the impervious area (C_A) value through the forward and reverse deduction of the probability matrix of state transition. The specific steps are as follows:

Firstly, for remote sensing images meeting the accuracy requirements, ENVI was used to perform radiometric calibration and quick atmospheric correction, and then supervised classification was used to interpret different land-use types. Secondly, the data of land-use types of the control years are used to calibrate the parameters of the CA-Markov model to simulate the area of different land-use types in the years lacking the remote sensing images or not meeting the accuracy requirements.

The transition probability matrix and transition area matrix of the land-use types are expressed, respectively, as follows:

(1) ${\mathit{F}}_{t+1}={\mathit{P}}_{ij}\times {\mathit{F}}_t$(1)

(2) $S_{t+1}=\mathrm{f}(S_t,N)$(2)

where F_t and F_t+1 are the area matrix of the land-use types at time t and t+ 1, respectively; P_ij is the probability matrix of state transition; S_t and S_t+1 denote the state set of the cellular automation at time t and t+ 1, respectively; N is the neighbourhood of the cellular automation; and f is the transfer function that defines the transfer of the cellular automation from time t to time t+ 1.

3.2 Identification of the joint change point based on copula

3.2.1 Selection of the multivariate copula functions

In order to identify the change points of a multivariable long sequence, we first need to establish a multivariate joint distribution. The copula theory proposed by Sklar in 1959 can connect any d-dimensional cumulative distributions with a copula function. The Archimedean copula function is commonly used in the hydrological field (Guo et al. 2008, Zhang et al. 2011). In general, the d-dimensional symmetric Archimedean copula function is defined as follows:

(3) $C(u_1,u_2,...,u_d)={\phi }^{-1}[\phi (u_1)+\phi (u_2)+...+\phi (u_d)]$(3)

where φ(~) is the generation function of the Archimedean copula function; and φ⁻¹ is the inverse function of φ.

Typical symmetrical Archimedean copula functions include Clayton, Gumbel-Hougaard (G-H), Frank and Ali-Mikhail-Haq (AMH). The construction of the copula joint distribution mainly involves the following steps:

Step 1. The distribution of hydrological variables is fitted to determine the optimal type of each marginal distribution.

Step 2. According to the distribution type of each marginal function, the parameters θ of the four symmetric Archimedean copula functions are calculated and the copula function with the optimal fitting is determined.

3.2.2 Identification of the joint change points

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

${\mathrm{H}}_0:{\theta }_1={\theta }_2=\cdots ={\theta }_n,{\mathrm{H}}_1:{\theta }_1=\cdots ={\theta }_n{}^{\ast }\ne {\theta }_{n_{\hspace{0.17em}}^{\ast }+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 the copula can be established based on the maximum likelihood estimation method, as shown by:

(4) $\begin{array}{rl}& -2ln({\mathrm{\Lambda }}_k)=2\left[\sum _{i=1}^k{C}_{uvw}({\theta }_k;U(x_i),V(y_i),W(z_i))\right.\\ & +\sum _{i=k+1}^n{C}_{uvw}({\theta }_{k_{\hspace{0.17em}}^{\ast }};U(x_i),V(y_i),W(z_i))\\ & \left.\begin{array}{ccc}\begin{array}{cc}\hspace{0.17em}& \hspace{0.17em}\end{array}& \hspace{0.17em}& \hspace{0.17em}\end{array}+\sum _{i=1}^n{C}_{uvw}({\theta }_n;U(x_i),V(y_i),W(z_i))\right]\end{array}$(4)

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

The statistic $Z_n=\underset{1\le k<n}{max}(-2ln({\mathrm{\Lambda }}_k))$ is constructed and, if k* is unknown, the original hypothesis will be rejected when z_n is greater than the critical value; that is, the Archimedean copula function has a change point. The binary segmentation method will be applied in the test if there is more than one change point in the sequence (Jarušková 2015).

3.3 Combined exceedance probability of discharge

According to the full probability formula and the copula functions (Song and Nie 2011), we can quantitatively analyse the exceedance probability, i.e. that the discharge exceeds specific values in the combined conditions of precipitation and land-use types. Taking the joint probability distributions of two and three dimensions as an example, the expressions are as follows:

(5) $P_{VW}=P(V\le v,W\ge w)=F_V(v)-F_{VW}(v,w)=F_V(v)-C_{\hspace{0.17em}}(v,w)$(5)

(6) $\begin{array}{rl}& P_{UVW}=P(U\le u,V\le v,W\ge w)\\ & \begin{array}{cc}\hspace{0.17em}& \hspace{0.17em}\end{array}=F_{UV}(u,v)-F_{UVW}(u,v,w)\\ & =C_{\hspace{0.17em}}(u,v)-C_{\hspace{0.17em}}(u,v,w)\end{array}$(6)

where V and U denote the affecting factors; v and u are their specific values; W denotes discharge and w is its specific value; P_VW is the bivariate joint probability distribution of the pair (P,Q), P_UVW is the trivariate joint probability distribution of set (C_A,P,Q); F_U(u), F_V(u) and F_W(u) are the marginal distributions of V, U and W, respectively; and C(u,v,w) is the copula function of the variables.

4 Results and discussion

4.1 Change points and trend of univariate distribution

The collected remote sensing images are preprocessed, and supervised classification done by ENVI, and then they are clipped by the boundary of the study area, which is extracted by ArcGIS, to obtain the areas of land-use types. Taking 2 years as the calculation step of land transfer, the land-use types in the year 2033 can be forecast and C_A of 0.64 is obtained. Some typical land-use types are shown in Fig. 3 and the ratio of impervious area (C_A) from 1970 to 2016 is shown in Fig. 4.

Figure 3. Land-use types of the Xiaoqing River basin above the Huangtai hydrological station in typical years.

Figure 4. The ratio of impervious area of the Xiaoqing River basin above the Huangtai hydrological station (1970–2016).

As can be seen from Fig. 4, the value of C_A has been increasing in the study basin. The Mann-Kendall (M-K) test (Yue and Wang 2004) is used to analyse the trend and identify the change points of the data series of annual precipitation (P), discharge (Q) and impervious area (C_A). The results are shown in Table 2. The annual discharge and C_A show a clear trend of growth, while the trend in annual precipitation is not obvious. In addition, the land-use types changed first in 1992, and the first change point of the discharge occurred in 1994, while the precipitation is the last one to change in 2003.

Table 2. Identification of the trend and change points of univariate distribution of the Xiaoqing River basin above the Huangtai hydrological station in the period 1970–2016. ↑ and ↓denote an increasing or decreasing trend; + and – represent a significant and insignificant trend, respectively.

	P	Q	CA
Trend	↑–	↑+	↑+
Change point	2003	1994,2004	1992,2009

4.2 Construction of multivariate joint distribution and analysis of joint change points

The Clayton, G-H and Frank copula functions are used to construct the bivariate and trivariate joint probability distributions with the data of P, C_A and Q of Xiaoqing River in Jinan for the period 1970–2016. The proper distribution was chosen to use the Kolmogorov-Smirnov (K-S) test and root mean square error (RMSE) (Justel et al. 1997, Willmott and Matsuura 2005). The parameters of the copulas and the results of the goodness-of-fit tests are shown in Table 3. It can be seen from Table 3 that the Clayton copula is proper for sets of (P,Q) and (C_A,Q), and the G-H copula is proper for sets of (C_A,P) and (C_A,P,Q). The corresponding contour plots and contour surface are presented in Fig. 5.

Table 3. Parameter values and test results of multivariable copula functions, Clayton, Gumbel-Hougaard (GH) and Frank. RMSE: root mean square error; K-S: Kolmogorov-Smirnov test. Bold denotes the minimum RMSE and K-S, and the best fitness parameter.

	Clayton			GH			Frank
	K-S	RMSE	θ	K-S	RMSE	θ	K-S	RMSE	θ
P,Q	0.0785	0.0318	1.6955	0.0953	0.0345	1.8477	0.0854	0.0635	4.2661
CA,Q	0.1014	0.0504	0.5256	0.1144	0.0572	1.2628	0.1302	0.0944	2.0648
CA,P	0.3771	0.2889	0.1772	0.1224	0.0591	1.3359	0.2106	0.0701	1.0124
CA,P,Q	0.1099	0.0541	0.5212	0.1077	0.0507	1.3231	0.1086	0.0589	2.4097

Figure 5. Contour plots and surface of bivariate and trivariate joint probability distribution for the period 1970–2016. Bivariate joint probability of: (a) set (P,Q); (b) set (C_A,Q); (c) set (C_A,P); and (d) trivariate joint probability of set (C_A,P,Q).

Combining the optimal copula functions of the two- and three-variable sets and Equation (4), their joint change points can be identified. When the number of data is 46, the critical value of Z_n is 8.82. The corresponding results are shown in Fig. 6 and Table 4, which indicate that the pairs (P,Q), (C_A,P) and set (C_A,P,Q) have only one change point in 2000, 2002 and 2007, respectively, while there are two change points for the pair (C_A,Q), in 1992 and 2003.

Table 4. Change points and tail-dependence of two- and three-variable parameter sets.

	P,Q	CA,Q	CA,P	CA,P,Q
Change point	2000	1992,2003	2002	2007

Figure 6. Analysis of the combined change points based on two- or three-variable sets: (a) pair (P,Q); (b) pair (C_A,Q); (c) pair (C_A,P); and (d) set (C_A,P,Q).

4.3 Analysis of multivariate joint distribution of C_A, P and Q before and after the joint change point

The data of C_A, P and Q after the joint change point is more suitable for representing hydrological processes in recent years. Therefore, we re-established the multivariate joint probability distribution using the data after the joint change point. In addition, in order to better reflect the change law of hydrological processes, we also re-established the copula before the joint change point. The parameters of the copulas are shown in Table 5.

Table 5. Selected copulas and their parameters before and after the combined change point. G-H: Gumbel-Hougaard, AMH: Ali-Mikhail-Haq.

Parameter set	Before the joint change point (1970–2000)		After the joint change point (2000–2016)
	Copula	θ	Copula	θ
(P,Q)	Clayton	1.109	AMH	0.041
(CA,P,Q)	G-H	1.343	G-H	1.112

With the constructed joint probability distributions of (P,Q) and (P,C_A,Q) before and after the joint change point, the combined exceedance probability of discharge of the Xiaoqing River in Jinan can be analysed. According to Equations (5) and (6), P_VW is the bivariate joint probability distribution of the pair of (P,Q), which represents the exceedance probability of discharge with different precipitations, while P_UVW is the trivariate joint probability distribution of the set (P,C_A,Q), which represents the combined exceedance probability of discharge in combination with precipitation and land-use types. The contour plots of P_VW and contour surfaces of P_UVW before and after the joint change point are shown in Figs. 7 and 8, respectively.

Figure 7. Contour plots of P_VW(v,w) probability: (a) with all the data of 1970–2016; (b) before the change point (1970–2000); and (c) after the change point (2000–2016).

Figure 8. Contour surface plots of P_UVW(u,v,w) probability: (a) with all the data of 1970–2016; (b) before the change point (1970–2007); and (c) after the change point (2007–2016).

4.4 Discussion

By comparing the change points of the univariate, bivariate and trivariate distribution, we can see that their occurrence time does not coincide. The change point of the single index occurred the earliest and that of the trivariate indexes occurred the latest. The reason is the single index only considers the changes of its own sequence without taking the effects of other factors into account. For example, for the Q series, although the discharge suddenly increased in 1994, the precipitation also significantly went up from 541.1 to 834.2 mm in 1993, so the dependency relationship of pair (P,Q) may still remain stable. Therefore, the change point of the single discharge sequence will not lead to a mutation in the dependency relationship of (P,Q). On the contrary, even if the P series has not mutated, the dependency relationship of (P,Q) may also be mutated. Similarly, when it comes to (P,Q) and (C_A,P,Q), the dependency relationship of (P,Q) mutated in 2000, yet if the land-use type is taken into consideration, the dependency relationship between the three-dimensional variables would not have mutated until 2007. Therefore, compared with the change point tests of univariate and bivariate distribution, the test of trivariate distribution considers the joint effect of different factors more comprehensively, which reduces the impact of individual variable and describes the changes of discharge more objectively under different levels of land development and precipitation.

From Fig. 7, we can see that the contour plot spacing of Fig. 8(c) is smaller than that of Fig. 8(a) and (b). This indicates that the discharge is more sensitive to changes in precipitation when the dependence relationship of (P,Q) has mutated. Besides, from Fig. 7, it can be seen that the contour plot of (P,Q) after the change point has an obvious defect: with the increase of P, especially when P is larger than 800 mm, the exceedance probability of Q basically remains unchanged, which indicates that the two-dimensional joint distribution cannot well reveal the change law of discharge in Xiaoqing River in Jinan.

From Fig. 8, we can see that the value of joint probability tends to be large with more precipitation and greater C_A. When the dependence relationship of (P,Q,C_A) passes through the joint change point, the contour surface spacing of Fig. 8(c) becomes smaller than that of Fig. 8(a) and (b). Meanwhile, Fig. 8 does not show the same defect as Fig. 7, which suggests that the established multivariate probability distribution can quantitatively evaluate the exceedance probability of discharge for a wider range of conditions.

We can further obtain some typical values of the joint probability of P, C_A and Q to evaluate the combined exceedance probability of discharge and to provide urban managers with more concise guidance. Typical values of joint probability are shown in Tables 6 and 7. They describe the risk probability of discharge exceeding a specific value under the conditions that only P is considered or that P and C_A are considered in combination. For example, under the condition that only P is considered, P_VW = P(V ≤ 800, W ≥ 12) = 0.1506, while under the condition that P and C_A are both considered, P_UVW = P(U ≤ 800, U ≤ 0.5, W ≥ 12) = 0.1727. The results show that the trivariate joint probability of the parameter set (C_A,P,Q) is larger than the bivariate joint probability of pair (P,Q). Therefore, if only the discharge and precipitation are considered, it cannot reasonably reflect the relationship between river water and water quality. The trivariate joint probability distribution considering Q, P and C_A can more comprehensively express the characteristics of the change law of hydrological processes.

Table 6. Joint probability of (P,Q).

P (mm)	Q (m^3/s)
	8	10	12	14	16
600	0.2003	0.1165	0.0574	0.0255	0.0107
700	0.3561	0.2118	0.1062	0.0477	0.0201
800	0.4847	0.2943	0.1506	0.0686	0.0291
900	0.5677	0.3504	0.183	0.085	0.0366
1000	0.6129	0.3827	0.2033	0.0964	0.0422

Table 7. Joint probability of (P,C_A,Q).

P (mm)	CA	Q (m^3/s)
		8	10	12	14	16
600	0.5	0.2298	0.0189	-	-	-
	0.55	0.2437	0.0731	0.057	0.0132	-
	0.6	0.2514	0.1163	0.0641	0.0197	0.0024
	0.65	0.2554	0.135	0.068	0.0234	0.006
700	0.5	0.4286	0.0571	-	-	-
	0.55	0.4545	0.1568	0.1232	0.0343	-
	0.6	0.4689	0.2349	0.1354	0.0449	0.0085
	0.65	0.4763	0.2689	0.1422	0.0511	0.0144
800	0.5	0.5998	0.0977	0.1727	0.0452	-
	0.55	0.6369	0.2408	0.1989	0.0643	0.0075
	0.6	0.6575	0.3515	0.2155	0.0775	0.0191
	0.65	0.6682	0.4003	0.225	0.0855	0.0263
900	0.5	0.7112	0.127	0.2334	0.0765	0.0079
	0.55	0.7568	0.3016	0.265	0.0968	0.0217
	0.6	0.7823	0.4369	0.2854	0.1116	0.0329
	0.65	0.7956	0.4977	0.2971	0.1208	0.0406
1000	0.5	0.7708	0.1436	0.2746	0.1026	0.0252
	0.55	0.8217	0.3364	0.3104	0.1243	0.0372
	0.6	0.8506	0.4868	0.3338	0.1403	0.0477
	0.65	0.8658	0.5559	0.3474	0.1506	0.0553

5 Conclusions

This study analyses the main influence and response indicators of the hydrological processes of the Xiaoqing River in Jinan. By building a three-dimensional joint distribution of the set (C_A,P,Q), this paper compares the occurrence time of change points of one-, two- and three-variable sets as well as the variation in the joint probability; thus urban managers can evaluate the given design value risk exceeded by a combination of events (C_A,P,Q), as well as the changes in risk values before and after the joint change point. The outcome of the study may be summarized as follows:

1. Due to taking the relationship between precipitation, land-use types and discharge into account, the test of multivariable joint change point is less affected by an individual index. Thus, this test can more accurately identify the change point in hydrological processes of urban rivers.

2. Urbanization has an obvious impact on the hydrological processes of rivers. Under the condition of little change in annual precipitation in the two change periods, no matter the bivariate and trivariate joint distributions, the joint probability distribution of discharge, and other factors is obviously different in the two periods, which indicates that human activities are having an impact on hydrological processes.

3. In the same case of precipitation, the exceedance probability of discharge calculated by three-dimensional joint distribution is greater than the value calculated by two-dimensional joint distribution. The three-dimensional joint distribution considers the factors more comprehensively, so it is more objective than the two-dimensional joint distribution to reflect the occurrence probability of hydrological variation of urban rivers.

4. Due to the limitations of data acquisition, especially that of land-use types, the sequence length in this paper is only 47 years, which results in the data sequence length after change points being short and the joint probability value calculated may have some deviation. Future work should continue to supplement the data in the future to improve the research on hydrological changes of urban rivers.

Acknowledgements

The authors are grateful for the insightful suggestions made by the editors and the anonymous reviewer, which helped to improve the paper. We are grateful for the basic data provided for this study by the Water Resources Department of Shandong Province in China.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was funded by the Public Welfare Industry Research Special Fund Project of Ministry of Water Resources of China “Key Technology Research and Demonstration Project of Water Ecological Civilization Construction” [Grant no. 201401003].

Notes

¹ http://www.gscloud.cn/.