Low Flow Regimes of the Tarim River Basin, China: Probabilistic Behavior, Causes and Implications

Abstract

Droughts are a frequent occurrence in Xinjiang, China, and therefore fundamental to determining their hydrologic characteristics is low flow analysis. To that end, 11 probability distribution functions and 26 copulas functions were employed to analyze the changing characteristics of low flow regime (defined as seven-day low flow) of the Tarim River Basin. Results indicated that: (1) The Wakeby distribution satisfactorily described the probabilistic behavior of the low flow regime. According to Akaike Information Criterion (AIC), Bayesian Information Criterions (BIC), maximum likelihood, and other residual-based metrics, Tawn copula, Farlie–Gumbel–Morgenstern copula and Frank copula were the best choice and used in this current study. (2) After 1987, hydrological droughts of longer return periods were prone to higher occurrence frequency. (3) The low flow volume has been increasing in recent years due to the temperature-induced increase of snowmelt and increasing precipitation. However, hydrological droughts can be expected to occur due to the massive increase in water demand from the development of irrigated agriculture, increasing arable land and livestock farming. As a result, the water shortage in the lower Tarim River Basin will be increasingly severe under the influence of climate change and human activities. To alleviate the shortage would call for the development of water-saving agricultural irrigation, water-saving technology, conservation of eco-environment and sustainable development of local socio-economy.

1. Introduction

Traditionally, hydrologists have paid greater attention to floods that are normally associated with more visible and dramatic hazards, damages and economic losses [1]. This has been especially true in China, partially due to the recurrence of severe floods in major river basins in the past few decades [2,3,4]. However, low-flow information for both gauged and ungauged sites is critical for a wide range of applications, including the planning and design of water supply systems, analysis of environmental and economic impacts, modeling of stream water quality, and determination of optimum maintenance flow for instream use and hydropower operation [5,6]. Low flow regime is also a major component in the drought conception [7,8]. Gumbel (1963) defined a drought as the smallest annual value of daily streamflow [9]. Palmer (1965) described a drought as a significant deviation from the normal hydrologic condition of an area [10]. Therefore, an understanding and estimation of low-flow characteristics are vital to the sustainable development and management of water resources and also for monitoring practices of drought hazards [11,12,13].

Northwest China is characterized by arid and semi-arid climate. Availability and variability of water resources have a direct influence on local eco-environmental conservation and sustainable socio-economic development. The Tarim River (Figure 1), with an annual flow of 4–6 billion cubic meters, is the longest inland river in China. The basin has a population of about 10 million, including ethnic minorities, such as Uyghurs and Mongolians. The climate of the basin is characterized by low precipitation and high evaporation characteristic of arid climate. The temporal distribution of precipitation throughout the year is strongly heterogeneous. More than 80% of the total annual precipitation falls between May and September in the high flow season and less than 20% of the total falls from November to the following April [14]. Precipitation in the mountainous regions can exceed 300 mm per year in some areas and is mostly in the form of snowfall [15]. Water resources play a key role in the conservation of ecological environment and sustainable development of social economy and agriculture is heavily dependent on irrigation. Sustainable water distribution depends on average multi-year discharges and water requirements for ecological and economic development in the Tarim River Basin [16]. This requires an investigation of the hydrological processes and water resources availability in the basin where the supply of water resources is a major constraint for further socio-economic development and ecological protection.

There have been several studies dealing with climate change and water resources in the Tarim River Basin. Using 50 years of data on hydrology, temperature, and precipitation, Chen et al. (2007) investigated the effect of climate change on water resources [14]. Xu et al. (2008) analyzed the relationship between ecological change and agricultural development and discussed water resources management and its ecological significance [16]. Analyzing the characteristics of hydro-climatic changes, Tao et al. (2011) found that local human activities since the 1970s led to a decrease of the water volume diverted into the main stream of the basin, which was aggravated in the 2000s [15]. Based on daily precipitation and temperature dataset of 23 gauging stations and annual streamflow data of five hydrological stations covering 1960–2005, Zhang et al. (2010) analyzed spatio-temporal variations of climatic change and associated impacts on water resource and found that annual streamflows of headstream rivers, except Aksu River, did not exhibit significant trends [17,18]. Gu et al. (2016) analyzed the intra-annual, seasonal and inter-seasonal clustering of floods and investigated possible impacts of climate indices on the occurrence rates of floods [19,20]. Li et al. (2016) developed a new agricultural drought index, which is based on the Variable Infiltration Capacity (VIC) model coupled with an irrigation scheme and a reservoir module. The new drought index was derived from the simulated soil moisture data from a retrospective VIC simulation from 1961 to 2007 over the irrigated area [21].

Sources of streamflow in the Tarim River Basin are precipitation and glacial melting, and hence an increase in streamflow is the result of increasing precipitation and/or accelerating ice melting. Xu et al. [22] detected trends in major hydrometeor-hydrological variables during the period of 1960–2007. Results showed that both mean annual air temperature and precipitation experienced an increasing trend, while annual streamflow demonstrated a mixed trend of decreasing and increasing: The mountainous region upstream showed an increasing trend and the region downstream exhibited a decreasing trend.

However, despite its importance to the development of agriculture and conservation of eco-environment, only a few reports are available addressing the changing properties of low flow regimes. This constitutes the motivation of this study and therefore the objective of this study was to investigate the probabilistic behavior of hydrological droughts represented by low flows and to discuss the possible causes behind the changing properties of low flows.

2. Data

The hydrological data analyzed in this study were daily streamflow for the period of 1962–2008 from eight major hydrological stations, as shown in

Table 1. Information of seven-day low flow at eight stations of the Tarim River Basin.

No.	Stations	Abbreviation	Streamflow Series	Basin
1	Tongguziluoke	TG	1962–13 December 2008	Hotan river
2	Yuzimenluoke	YZ	1 January 1962–13 December 2008	Yarkand river
3	Kaqun	KQ	1 January 1962–13 December 2008	Yarkand river
4	Shaliguilanke	SL	1 January 1962–13 December 2008	Aksu river
5	Xiehela	XH	1 January 1962–13 December 2008	Aksu river
6	Dashankou	DS	1 January 1972–13 December 2008	Kaidu river
7	Huangshuigou	HS	1 January 1972–13 December 2008	Kaidu river
8	Alaer	AL	1 January 1962–13 December 2008	Mainstrean of Tarim river

: Tongguziluoke (TG), Yuzimenleke (YZ), Kaqun (KQ), Shaliguilanke (SL), Xiehela (XH), Huangshuigou (HS), Dashankou (DS) and Alaer (AL). The streamflow data were obtained from the Management Bureau of the Tarim River Basin. The missing data were processed based on the method by Zhang et al. [17,18]. Specifically, the missing values of 1–2 days were replaced by the average of neighboring days. Adequate consecutive days with missing data were filled by the long-term average of the same days of other years. In this study, the low flow regime was defined as the seven-day low flow [17,23,24]. In addition, information on water reservoirs and irrigated areas was also collected. Locations of the hydrological stations, water reservoirs and irrigated areas are shown in Figure 1.

3. Methods

3.1. Frequency Analysis

First, an autocorrelation analysis was conducted to confirm the independence of observations. Then, the probabilistic behavior of hydrological extremes was analyzed using the following steps: (1) Eleven probability distribution functions (PDFs), commonly used in extreme value analysis, were used [25,26]. (2) For each individual station, these 11 probability distribution functions were fitted to the seven-day low flow series with their parameters estimated using the L-moment estimation [27]. (3) Goodness-of-fit was judged by the Kolmogorov–Smirnov statistic D (K-S D). Here, 95% confidence level was accepted to reject or accept a fit. (4) Based the on K-S D value, the PDF better describing the probabilistic behavior of seven-day low flow was selected for each station. The PDF that fitted the seven-day low flow series of the most of the stations was adopted as the best choice.

3.2. Copula Function

The copula functions have been used in the analysis of hydro-meteorological extremes [28,29,30,31,32]. In some studies, the Archimedean copula family was used to analyze the joint probabilistic behavior of hydrological extremes [30]. The reason for choosing this copula family is that it can be easily constructed, a large variety of copulas belong to this family, and it can be applied when the correlation amongst hydrologic variables is positive or negative [30,33]. The copula functions considered in this study were the 26 Copulas in

Table 2. Copula Families and Their Closed-Form Mathematical Description [34].

Name	Abbreviation	Mathematical Descriptiona	Parameter Range	Reference
Gaussian	Gau	${\displaystyle {\int }_{-\infty }^{{\varphi }^{-1}(u)}{\displaystyle {\int }_{-\infty }^{{\varphi }^{-1}(v)}\frac{1}{2\pi \sqrt{1-{\theta }^2}}}}\exp (\frac{2\theta xy-x^2-y^2}{2(1-{\theta }^2)})dxd{y}^b$	$\theta \in \left[-1,1\right]$	[38]
t	t	${\displaystyle {\int }_{-\infty }^{t_{{\theta }_2}^{-1}(u)}{\displaystyle {\int }_{-\infty }^{t_{{\theta }_2}^{-1}(v)}\frac{\Gamma \left(\left({\theta }_2+2\right)/2\right)}{\Gamma ({\theta }_2/2)\pi {\theta }_2\sqrt{1-{\theta }_1^2}}}}{(1+\frac{x^2-2{\theta }_{1}xy+y^2}{{\theta }_2})}^{({\theta }_2+2)/2}dxd{y}^c$	${\theta }_1\in \left[-1,1\right]and{\theta }_2\in \left(0,\infty \right)$	[38]
Clayton	Cla	$\max {(u^{-\theta }+v^{-\theta }-1,0)}^{-1/\theta }$	$\theta \in \left[-1,\infty \right]\backslash 0$	[39]
Frank	F	$-\frac{1}{\theta }\ln \left[1+\frac{(\exp (-\theta u)-1)(\exp (-\theta v)-1)}{\exp (-\theta )-1}\right]$	$\theta \in R\backslash 0$	[38]
Gumbel	Gum	$\exp \left\{-{\left[{(-\ln (u))}^{\theta }+{(-\ln (v))}^{\theta }\right]}^{1/\theta }\right\}$	$\theta \in \left[1,\infty \right)$	[38]
Independence	I	$uv$		[40]
Ali–Mikhail–Haq	AMH	$\frac{uv}{1-\theta (1-u)(1-v)}$	$\theta \in \left[-1,1\right)$	[41]
Joe	Joe	$1-{\left[{(1-u)}^{\theta }+{(1-v)}^{\theta }-{(1-u)}^{\theta }{(1-v)}^{\theta }\right]}^{1/\theta }$	$\theta \in \left[1,\infty \right)$	[38]
Farlie–Gumbel–Morgenstern	FGM	$uv\left[1+\theta (1-u)(1-v)\right]$	${\theta }_1\in \left(0,\infty \right),{\theta }_2\in \left(1,\infty \right)$	[40]
Gumbel–Barnett	GB	$\mathrm{u}+v-1+(1-u)(1-v)\exp \left[-\theta \ln (1-u)\ln (1-v)\right]$	$\theta \in \left[0,1\right]$	[42,43]
Plackett	P	$\frac{1+(\theta -1)(u+v)-\sqrt{{\left[1+(\theta -1)(u+v)\right]}^2-4\theta (\theta -1)uv}}{2(\theta -1)}$	$\theta \in \left(0,\infty \right)$	[44]
Cuadras-Auge	CA	${\left[\min (u,v)\right]}^{\theta }{(uv)}^{(1-\theta )}$	$\theta \in \left[0,1\right]$	[45]
Raftery	R	$\{\begin{array}{c}u-\frac{1-\theta }{1+\theta }{u}^{\frac{1}{1-\theta }}(v^{\frac{-\theta }{1-\theta }}-v^{\frac{1}{1-\theta }}),if\mathrm{u}\le \mathrm{v}\\ v-\frac{1-\theta }{1+\theta }{v}^{\frac{1}{1-\theta }}(u^{\frac{-\theta }{1-\theta }}-u^{\frac{1}{1-\theta }}),if\mathrm{v}\le \mathrm{u}\end{array}$	$\theta \in \left[0,1\right)$	[40]
Shih-Louis	SL	$\begin{array}{l}\{\begin{array}{l}(1-\theta )uv+\theta \min (u,v),if\theta \in \left(0,\infty \right)\\ (1+\theta )uv+\theta (u+v-1)\Psi (u+v-1),if\theta \in \left(-\infty ,0\right]\end{array}\\ \Psi (a)=1\mathrm{if}\mathrm{a}\ge 0\mathrm{and}\Psi (a)=0\mathrm{if}\mathrm{a}\text{<}0\end{array}$		[46]
Linear-Spearman	LS	$\begin{array}{c}\left[u+\theta (1-u)\right]v,\mathrm{if}\mathrm{v}\le \mathrm{u}\mathrm{and}\theta \in \left[0,1\right]\\ \left[v+\theta (1-v)\right]u,\mathrm{if}\mathrm{u}\text{<}1\mathrm{and}\theta \in \left[0,1\right]\\ (1+\theta )uv,\mathrm{if}\mathrm{u}+\mathrm{v}\text{<}1\mathrm{and}\theta \in \left[-1,0\right]\\ uv+\theta (1-u)(1-v),\mathrm{if}\mathrm{u}+\mathrm{v}\ge 1\mathrm{and}\theta \in \left[-1,0\right]\end{array}$	$\theta \in \left[-1,1\right]$	[47]
Cubic	Cu	$uv\left[1+\theta (u-1)(v-1)(2u-1)(2v-1)\right]$	$\theta \in \left[-1,2\right]$	[48]
Burr	B	$u+v-1+{\left[{(1-u)}^{-1/\theta }+{(1-v)}^{-1/\theta }-1\right]}^{-\theta }$	$\theta \in \left(0,\infty \right)$	[49]
Nelsen	N	$\frac{-1}{\theta }\log \left\{1+\frac{\left[\exp (-\theta u)-1\right]\left[\exp (-\theta v)-1\right]}{\exp (-\theta )-1}\right\}$	$\theta \in \left(0,\infty \right)$	[40]
Galambos	Ga	$uv\exp {\left\{{(-\ln (u))}^{-\theta }+{(-\ln (v))}^{-\theta }\right\}}^{-1/\theta }$	$\theta \in \left[0,\infty \right)$	[50]
Marshall–Olkin	MO	$\min \left[u^{(1-{\theta }_1)}v,u{v}^{(1-{\theta }^2)}\right]$	${\theta }_1,{\theta }_2\in \left[0,\infty \right)$	[50]
Fischer–Hinzmann	FH	${\left\{{\theta }_1{\left[\min (u,v)\right]}^{{\theta }_2}+(1-{\theta }_1){\left[uv\right]}^{{\theta }_2}\right\}}^{1/{\theta }_2}$	${\theta }_1\in \left[0,1\right],{\theta }_2\in R$	[51]
Roch–Alegre	RA	$\exp \left\{1-{\left[{({({(1-\ln (u))}^{{\theta }_1}-1)}^{{\theta }_2}+{({(1-\ln (v))}^{{\theta }_1}-1)}^{{\theta }_2})}^{1/{\theta }_2}+1\right]}^{1/{\theta }_1}\right\}$	${\theta }_1\in \left(0,\infty \right),{\theta }_2\in \left[1,\infty \right)$	[52]
Fischer–Kock	FK	$uv{\left[1+{\theta }_2(1-u^{1/{\theta }_1})(1-v^{1/{\theta }_1})\right]}^{{\theta }_1}$	${\theta }_1\in \left[1,\infty \right),{\theta }_2\in \left[-1,1\right]$	[34]
BB1	BB1	${\left\{1+{\left[{(u^{-{\theta }_1}-1)}^{{\theta }_2}+{(v^{-{\theta }_1}-1)}^{{\theta }_2}\right]}^{1/{\theta }_2}\right\}}^{-1/{\theta }_1}$	${\theta }_1\in \left(0,\infty \right),{\theta }_2\in \left(1,\infty \right)$	[35]
BB5	BB5	$\exp \left\{-{\left[{(-\ln (u))}^{{\theta }_1}+{(-\ln (v))}^{{\theta }_1}-{({(-\ln (u))}^{-{\theta }_1{\theta }_2}+{(-\ln (v))}^{-{\theta }_1{\theta }_2})}^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{${\theta }_2$}\right.}\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\theta }_1$}\right.}\right\}$	${\theta }_1\in \left[1,\infty \right),{\theta }_2\in \left(0,\infty \right)$	[35]
Tawn	Tawn	$\exp \left\{\ln (u^{(1-{\theta }_1)})+\ln (v^{(1-{\theta }_2)})-{\left[{(-{\theta }_1\ln (u))}^{{\theta }_3}+{(-{\theta }_2\ln (v))}^{{\theta }_3}\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\theta }_3$}\right.}\right\}$	${\theta }_1,{\theta }_2\in \left[0,1\right],{\theta }_3\in \left[1,\infty \right]$	[50]

[34]. The Ali–Mikhail–Haq family is similar to the Frank family but is limited by the correlation structure of two hydrological variables [35,36]. Sadegh et al. (2017) empirically estimated the marginal distribution of each hydrological variable and constructed the joint distribution of precipitation and soil moisture anomalies using the 26 copulas. They found that only the Tawn and Marshall–Olkin copulas are capable of characterizing the asymmetric dependence structure of this dataset. Neither is among the commonly used copulas in the hydrological literature. The Tawn copula provides a very good fit to the data [34,37]. Therefore, in this study, we construct the joint distribution of low flow based on 26 copulas.

3.3. Determination of the Generating Function and the Resulting Copula

The first step in determining a copula is to obtain its generating function from bivariate observations. The procedure to obtain the generating function and the resulting copula was introduced by Genest and Rivest [53] and it was followed in this study. The two variables here are (X, Y): (x₁, y₁), (x₂, y₂), …, (x_n, y_n). X and Y are the seven-day low flow series analyzed in this study. The procedure can be described here:

(1) Determination of Kendall’s τ based on the observed series:

$${\tau }_N={\left(\begin{array}{c}N\\ 2\end{array}\right)}^{-1}{\displaystyle \sum _{i<j}sign[(x_i-x_j)(y_i-y_j)]}$$

(1)

where N is the length of the series; sign = 1 when x_i ≤ x_j and y_i ≤ y_j, otherwise, sign = −1; i, j = 1, 2, …, N; and τ_N is the estimate of τ from the observed series.

(2) The copula parameter, θ, can be estimated based on Markov Chain Monte Carlo (MCMC) simulation within a Bayesian framework [54,55]. MCMC simulation estimates the posterior distribution of parameter values, which are then translated into uncertainty ranges for the copula probability isolines. The MCMC simulation searches for the region of interest with multiple chains are running in parallel. Chains share information on the fly, characterize the posterior region (even in the presence of multimodality), and estimate the global optimum [34,56,57].

3.4. Selection of Copula Family

The computation procedure for copula-based analysis of hydrological extremes involves the following steps [30]: (a) determination of marginal distributions based on the conventional statistical approach; (b) identification of generator and parameter of copulas; (c) determination of the joint probability distribution; and (d) application to real data. The marginal probability distribution function was evaluated using the Kolmogorov–Smirnov goodness-of-fit statistics (K-S D statistics). We use several goodness of fit measures to evaluate the performance of different copula models, including likelihood value, AIC, BIC, RMSE, and NSE [34]. Higher model complexity (more degrees of freedom) provides the advantage of greater model flexibility and hence usually results in a better fit to the observed data. However, this might stimulate over-conditioning of the model. AIC, in contrast to the ad hoc likelihood value, considers both complexity of the model and minimization of error residuals and provides a more robust measure of quality of model predictions. AIC avoids the problem of over-conditioning by adding a penalty term based on the number of parameters. AIC is formulated as [34,58,59,60]:

$$AIC=2D-2\ell$$

(2)

in which D is the number of parameters of the statistical model and $\ell$ is the log-likelihood value of the best parameter set. This equation can be simplified to:

$$AIC=2D+n\ln \left\{\frac{{\displaystyle \sum _{i=1}^n{\left[\overline{y_i}-y_i(\theta )\right]}^2}}{n}\right\}-2cs$$

(3)

$$\ell =-\frac{n}{2}\ln (2\pi )-\frac{n}{2}\ln {\sigma }^2-\frac{1}{2}\ln {\sigma }^{-2}{\displaystyle \sum _{i=1}^n{\left[{\tilde{y}}_i-y_i(\theta )\right]}^2}$$

(4)

Given the Gaussian assumption of error residuals, ${\sigma }^2=\frac{{\displaystyle \sum _{i=1}^n{\left[{\tilde{\mathrm{y}}}_i-y_i(\theta )\right]}^2}}{n}$, $\theta$ is Copula parameter and cs is a constant. A lower AIC value associates with a better model fit.

Similar to AIC, BIC is presented as [34,61]:

$$BIC=D\ln n-2\ell$$

(5)

which similarly simplifies to:

$$BIC=D\ln n+n\ln \left\{\frac{{\displaystyle \sum _{i=1}^n{\left[\overline{y_i}-y_i(\theta )\right]}^2}}{n}\right\}-2cs$$

(6)

If residuals are independent and identically distributed following a Gaussian distribution centered around zero, then, similar to AIC, a lower BIC value associates with a better model fit.

NSE and RMSE are also two widely used measures of goodness of fit, which only focus on minimization of residuals,

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left[{\tilde{\mathrm{y}}}_i-y_i(\theta )\right]}^2}}{n}}$$

(7)

$$NSE=1-\frac{{\displaystyle \sum _{i=1}^n{\left[{\tilde{\mathrm{y}}}_i-y_i(\theta )\right]}^2}}{{\displaystyle \sum _{i=1}^n{\left[{\tilde{\mathrm{y}}}_i-{\overline{\tilde{y}}}_i\right]}^2}}.$$

(84)

A perfect model fit is associated with $RMSE=0$, $RMSE\in \left[0,\infty \right)$, and $NSE=1$. All these metrics evaluate, in different ways, the performance of copulas in terms of how close modeled bivariate probabilities ($y$) are to their empirical observed counterparts ($\tilde{y}$). Although number of copula parameters (model complexity) impacts some of the evaluation metrics, parameter ranges have zero impact on them.

4. Results

4.1. Selection of the Marginal Distributions and Copula Functions

Parameters of the PDFs were estimated using the L-moment technique and the goodness-of-fit was tested by the K-S method (

Table 3. K-S statistic D for 11 probability distribution functions describing the statistical properties of extreme streamflow.

	Hydrological Stations	TG	YZ	KQ	SL	XH	DS	HS	AL
Probability Functions
Wakeby (5P)		0.043	0.058	0.059	0.046	0.057	0.089	0.060	0.071
Weibull (3P)		0.042	0.096	0.055	0.063	0.094	0.144	0.076	0.162
Gamma (3P)		0.041	0.083	0.054	0.058	0.085	0.130	0.070	0.128
Lognormal (3P)		0.041	0.075	0.049	0.057	0.081	0.106	0.075	0.182
Log-logistic (3P)		0.046	0.067	0.055	0.054	0.070	0.102	0.075	0.193
General Pareto (3P)		0.081	0.088	0.097	0.056	0.119	0.119	0.064	0.111
Gen. Extreme Value (3P)		0.045	0.062	0.055	0.057	0.083	0.094	0.077	0.074
Maximum Extreme Value (3P)		0.062	0.070	0.105	0.048	0.098	0.107	0.079	0.153
Beta (4P)		0.072	0.084	0.052	0.063	0.086	0.132	0.081	0.401
Gumbel Max Distribution (2P)		0.072	0.078	0.082	0.073	0.092	0.085	0.081	0.099
Gumbel Min Distribution (2P)		0.123	0.194	0.111	0.181	0.158	0.224	0.204	0.122

Note: The threshold values of K-S D are: 0.198 (n = 47, 1-α = 0.95); 0.200 (n = 46, 1-α = 0.95); 0.224 (n = 37, 1-α = 0.95). Smaller K-S D indicates better goodness-of-fit of probability functions describing the seven-day low flow regimes. Bold numbers in Table 3 denote higher goodness-of-fit degree when compared to alternative functions.

). It can be observed in Table 3 that, at the 95% confidence level, the Wakeby, gamma, lognormal, log-logistic and general extreme value distribution functions well described the probabilistic behavior of the seven-day low flow regime. However, comparatively, the Wakeby function had the smallest K-S statistic D values (Table 3) and

Table 4. Estimated parameters of the Wakeby function describing the probabilistic behavior of seven-day low flow regime using the L-moment technique.

Stations	α	β	γ	δ	ξ
TG	18.90	9.08	3.02	−0.41	2.63
YZ	41.28	31.35	1.02	−0.05	0.53
KQ	70.26	9.52	6.20	−0.35	29.67
SL	935.35	131.39	5.82	−0.24	0
XH	16.82	5.44	1.48	0.09	18.55
DS	1836.5	61.10	14.15	−0.22	0
HS	47.92	58.27	1.20	10.37	1.06
AL	21.78	1.99	2.57	0.19	−1.35

lists the estimated Wakeby distribution parameters. Figure 2 reinforces this point by the distribution and cumulative distribution curves. Thus, the Wakeby functions were used to describe the probabilistic behavior of the seven-day low flow regime in the Tarim River Basin. It should also be noted that the Wakeby distribution is a more flexible distribution than other distributions and is widely used in extreme value analysis practice [62]. Therefore, the Wakeby distribution is usually regarded as a first option in extreme value analysis [63].

During the 1980s, climatic experienced abrupt changes [64,65], as characterized by increasing temperature and precipitation and also by the decrease in dust storm days [64]. Besides, after the 1990s, the net increase of arable land reached 68.75 × 10⁴ hm² with an increasing magnitude of 3.8% [66]. Therefore, 1987 was taken as the demarcation line and the probabilistic behavior of the seven-day low flow regime was studied before and after 1987 (Figure 3 and

Table 5. Designed seven-day low flow (m³/s) and related return periods before 1987.

Stations	T = 2	T = 3	T = 5	T = 7	T = 10	T = 20	T = 30	T = 50	T = 70	T = 100
TG	5.53	4.94	4.33	3.98	3.67	3.23	3.06	2.91	2.85	2.80
YZ	2.26	2.09	1.92	1.78	1.63	1.33	1.20	1.07	1.01	0.97
KQ	39.5	37.8	36.5	35.9	35.5	35.0	34.8	34.7	34.5	33.8
SL	10.3	9.05	7.87	7.28	6.80	6.20	5.98	5.81	5.73	5.67
XH	22.1	21.3	20.5	20.1	19.8	19.4	19.3	19.2	19.14	19.11
DS	36.7	34.3	31.9	30.6	29.5	28.0	27.5	27.0	26.8	26.6
HS	2.36	2.15	1.99	1.92	1.87	1.81	1.79	1.77	1.73	1.66
AL	4.95	2.83	1.20	0.51	0	0	0	0	0	0

and

Table 6. Designed seven-day low flow (m³/s) and related return periods after 1987.

Stations	T = 2	T = 3	T = 5	T = 7	T = 10	T = 20	T = 30	T = 50	T = 70	T = 100
TG	7.93	7.18	6.62	6.38	6.12	5.93	5.78	5.62	5.53	5.46
YZ	3.14	2.71	2.42	2.30	2.22	2.12	2.09	2.06	2.03	1.94
KQ	43.6	41.8	39.4	37.5	35.3	31.4	29.7	28.1	27.4	26.8
SL	13.0	10.8	9.16	8.50	8.02	7.47	7.29	7.14	7.08	7.04
XH	23.6	22.9	22.3	21.9	21.2	19.7	18.8	17.9	17.5	16.2
DS	44.9	39.1	34.7	32.8	31.4	28.1	24.8	19.4	15.7	12.2
HS	3.34	2.95	2.57	2.39	2.24	2.07	2.01	1.96	1.94	1.92
AL	13.0	11.2	10.0	9.31	8.41	6.10	4.69	3.18	2.39	1.75

). It can be observed in Figure 3 and Table 5 and Table 6 that no evident differences could be found in the seven-day low flow magnitudes with return periods of <20 years. The seven-day low flow increased after the end of the 1980s at the Shaliguilanke, Alaer, Huangshuigou and Yuzimenleke stations, and decreased at the Kaqun, Xiehela, Dashankou, and Tongguziluoke stations. Besides, the magnitude of the seven-day low flow with return periods of <20 years after the end of the 1980s was larger than that after the 1980s at the Kaqun, Dashankou, Xiehela stations, and the reverse was expected for the seven-day low flow with return periods of >20 years. These changes imply a decreasing magnitude of the seven-day low flow or increasing probability of hydrological droughts at the Kaidu and Yarkand River Basins (Figure 1). However, increased magnitude of the seven-day low flow could be detected at the Tongguziluoke, Alaer, Huangshuigou, and Shaliguilanke stations, implying a decreasing probability of hydrological droughts in the river basins these hydrological stations are located in (Figure 1).

4.2. Joint Probabilistic Behavior of Seven-Day Low Flow Regimes at Different Stations

The joint distributions of seven-day low flow at all the stations are determined by copulas. In total, 26 copulas (Table 2) were selected as candidate copulas for the joint distribution of seven-day low flow in Tarim River Basin. The association between seven-day low flow is measured using Kendall’s τ coefficient, while the parameters of the corresponding copulas distributions are estimated based on Markov Chain Monte Carlo (MCMC) simulation within a Bayesian framework. Based on the goodness-of fit test introduced in Section 3.4, selection of copula functions was based on the AIC, BIC, and Maxi-Likelihood, and results are illustrated in Figure 4. It is worth noting that a visual inspection of the 26 copulas fitted to the seven-day low flow shows that the Tawn copula provides a very good fit to the data with a NSE = 0.9978 (NSE = 1 is associated with a perfect fit) and is selected as the best copula according to AIC, BIC, maximum likelihood, and other residual-based metrics in Yarkand River Basin. Farlie–Gumbel–Morgenstern copula is a very good fit to the data with minimum value of AIC, BIC and Maximum likelihood in Aksu River Basin. Moreover, Frank copula, with minimum value of AIC, BIC and Maximum likelihood, is the best copula and will be used for the further analysis in Kaidu River Basin. The parameters of best copulas in Tarim River Basin based on MCMC are shown in

Table 7. Estimated parameters of the copula functions with the highest modelling performance.

Joint Distribution of Seven-Day Low Flow	Copula	θ1	θ2	θ3
YZ and KQ stations	Tawn	0.394	0.931	3.014
SL and XH stations	Farlie–Gumbel–Morgenstern	1.874	-	-
DS and HS stations	Frank	6.715	-	-

.

The joint probabilistic behavior of seven-day low flow regimes was analyzed for Yuzimenleke vs. Kaqun stations (with time interval of 1962–2008), Shaliguilanke vs. Xiehela stations (with the time interval of 1962–2008), and Dashankou vs. Huangshuigou stations (with time interval of 1972–2008). The marginal distribution function was the Wakeby function, with parameters estimated using the L-moment technique (Table 4). Figure 5, Figure 6 and Figure 7 show the joint probabilistic behavior of the seven-day low flow in the Yarkand River Basin, i.e., the Kaqun and the Yuzimenleke stations. The results of joint distribution, joint return periods and concurrent return periods are illustrated in Figure 5, Figure 6 and Figure 7, respectively. Results for other river basins considered in this study showed similar features and are not be given here. It can be seen in Figure 5, Figure 6 and Figure 7 that the probability was relatively small when low flow regimes with small streamflow magnitudes occurred simultaneously in both tributaries, so did the case when a large magnitude of low flow occurred in one tributary and a small magnitude of low flow in another tributary. However, the probability was increasingly higher when the magnitude of low flow increased in both tributaries.

Table 8. Joint return periods (JRP) and concurrence return periods (CRP) with different groups of hydrological stations and different groups of return periods.

Designed Return Periods T	YZ and KQ Stations		SL and XH Station		DS and HS Stations
	JRP T0	CRP Ta	JRP T0	CRP Ta	JRP T0	CRP Ta
2	1.5	3.1	1.6	2.7	1.4	3.6
5	3.1	13.7	3.4	9.6	2.9	19.9
10	5.7	42.5	6.2	25.5	5.4	74.7
20	10.8	131.4	11.7	67.6	10.4	288.8
30	15.9	254.3	17.2	119.5	15.4	641.9
50	26.1	584.3	27.8	245.1	25.4	1765.3
70	36.3	1010.6	38.4	393.4	35.4	3445.1
90	46.4	1521.6	48.9	560.2	45.4	5681.2
100	51.4	1806.4	54.2	649.6	50.4	7007.9

shows joint return periods (JRP) and concurrence return periods (CRP) with different groups of hydrological stations and different groups of return periods. It can be observed from the table that the JRP was smaller than the designed return period, and CRP was larger than the designed return period. Changes in JRP and CRP of low flows at the Yarkand and Aksu Rivers were largely in agreement. On the premium of the same occurrence frequency, JRP of low flows at the Kaidu River was larger than that at the Yarkand and Aksu Rivers; however, CRP of low flows at the Kaidu River was smaller than that at the Yarkand and Aksu Rivers. The concurrent occurrence probability was higher for the hydrological droughts with return periods of < 2 years at the Yarkand and Aksu Rivers, and for the hydrological droughts with return periods of < 5 years at the Kaidu River. In general, the Joint return periods (JRP) of hydrological droughts at the Kaidu River was smaller than that at the Yarkand and Aksu Rivers, implying a higher frequency of hydrological droughts in the Kaidu River Basin when compared to the Yarkand and Aksu River Basins. However, the concurrence return periods (CRP) of hydrological droughts at the Kaidu River was larger than that at the Yarkand and Aksu Rivers, implying a lower risk of hydrological droughts in the all Kaidu River Basin when compared to the Yarkand and Aksu River Basins.

5. Discussion

The low flow of the same return period at the Tongguziluoke, Yuzimenleke, Shaliguilanke, Huangshuigou and Alaer stations was larger in magnitude after 1987 than that before 1987. The low flow regime at the Kaqun station with return periods of larger than 10 years was smaller in magnitude after 1987 than that before 1987. At the Xiehela and Dashankou stations, however, the low flow regime with return periods of larger than 30 years was larger before 1987 than that after 1987. This observation is in agreement with the published results [18,64,67] that shifts occurred from warm-dry to warm-wet climatic conditions which may trigger increased volume of low flow regimes at most of the hydrological stations within the Tarim River Basin, except the Kaqun, the Xiehela and Dashankou stations which are dominated by decreased low flow volumes after about 1987. Streamflow of the Yarkand River is mainly from snow melt, which accounts for about 64% of the total streamflow measured at the Kaqun station; and 22.6% of streamflow is from ground water supply. Streamflow at the Kaqun station is mainly from snow melt and streamflow at the Yuzimenleke station is mainly due to snow melt and rainfall [68]. The temperature and precipitation in Xinjiang are increasing in general. In the Yarkand River Basin, the increase in temperature is evident in autumn when compared to other seasons. No observed changes in temperature can be found in winter. However, the increasing tendency of precipitation is obvious [68]. In this case, the increase in streamflow due to snowmelt is not evident. Thus, hydrological droughts are not alleviated, although precipitation is increasing because the streamflow of the Yarkand River Basin is mainly from snowmelt. The irrigation demand in the Yarkand River Basin is the largest within the Tarim River Basin (Figure 1). The water consumption by livestock farming and agriculture can reach 2.17 billion cubic meters per year, accounting for 28.4% of the total streamflow availability [69]. Although precipitation in the Yarkand River Basin is slightly on the increase, the increase in irrigation demand is even more significant. Thus, increase of precipitation cannot satisfy the water demand from agriculture and livestock farming.

The Shaliguilanke and Xiehela stations control the hydrological processes of the Aksu River Basin. The streamflow of the Aksu River Basin is mainly from rainfall and snowmelt. After 1987, precipitation in spring and winter was increasing and temperature in winter was evidently increasing. Thus, an increase of the low flow volume can be expected after the 1970s. The water demand within the Aksu River Basin is also large, although it is smaller than that in the Yarkand River Basin. The water consumption by livestock farming and agriculture irrigation can reach 1.49 billion m³, accounting for 17.6% of the total streamflow of this river basin. The streamflow supply of the Kaidu River depends heavily on rainfall and snow melting. The increase of streamflow due to snow melting accounts for 23.2% of the total streamflow [70]. The mean temperature during December, January, February, and March is −20.4 °C with weak evaporation and water supply of this river is mainly from ground water. However, precipitation changes in summer and autumn will directly impact the low flow amount in winter [71]. The increase of precipitation may alleviate the water deficit; however, increasing water demand due to massive agricultural irrigation may call for the mitigation of hydrological droughts to face new challenges if water-saving agriculture technology cannot be developed in due time (Figure 1).

It should be noted here that the low flow volume increased after 1987, and the occurrence frequency after 1987 also decreased when compared to that before 1987. However, Figure 8 shows that the drought-affected area was increasing at a rate of 2.3 thousand ha per year. The low flow changes are not in agreement with those of the drought-affected areas. It is attributed to the significant improvement in irrigation facilities and seepage prevention of the irrigation canal (Figure 1 and

Table 9. Detailed information of water reservoirs within the Tarim River Basin.

River	Number of Water Reservoirs	Total Storage Capacity/108m3	Utilizable Capacity/108m3	Cumulative Irrigated Areas/hm2	Effective Irrigated Areas/hm2	Designed Flood Volum/108m3
Hetian River Basin	20	2.35	2.05	36,263.04	33,268.3	2.20
Yarkand River Basin	37	14.20	11.57	303,816.3	201,283.4	19.97
Aksu River Basin	6	4.90	4.20	105,116.2	80,870.7	4.14
Kaidu River Basin	5	0.77	0.52	-	-	0.48
Mainstem of the Tarim River Basin	8	5.86	4.76	66,326.7	49,969.2	9.20
Total	76	28.08	23.10	511,522.2	365,391.6	35.99

and

Table 10. Detailed information of water supply facilities within the Tarim River Basin.

River	Number of Water Supply Facilities	Designed Irrigated Areas/hm2	Effective Irrigated Areas/hm2	Designed Water Supply Capacity/m3/s	Standing Water Supply Capacity/m3/s	Standing Water Supply Rate/%
Hetian River Basin	27	62,903.15	48,402.42	81.40	62.60	76.90
Yarkand River Basin	26	1,155,391	329,349.8	220.30	169.50	76.94
Aksu River Basin	63	702,035.1	570,028.5	198.60	165.50	83.33
Kaidu-Kongque River Basin	32	244,945.6	201,343.4	89.12	74.26	83.33
Mainstem of the Tarim River Basin	138 a	85,370.94	79,937.33	293.00	293.00	100.00
Total	286	2,250,646	1,229,061	882.42	764.86	86.68

^a: Provisional water withdrawal place.

). The agriculture of Xinjiang is heavily dependent on irrigation. The arable land for agriculture has increased within three time intervals: the increase of the arable land is about 0.449 million ha during 1949–1960; 0.265 million ha during 1963–1978; and 0.689 ha during 1990–2008 [66,72]. The irrigated arable land in the upper Tarim River Basin has increased from 0.348 million ha in 1950 to 1.645 million ha in 2007. The water consumption due to the significant increase of the irrigated arable land cannot be satisfied by the increase of streamflow as a result of increased precipitation in recent years. Thus, streamflow within the main stem of the Tarim River Basin is increasing correspondingly due to increased streamflow as a result of temperature-induced increase of snowmelt and increase of precipitation. Precipitation is increasing and snowmelt is also increasing due to increased temperature; however, water demand is significantly increasing due to the rapid development of agricultural irrigation. Therefore, droughts in Xinjiang are still increasingly serious and severe and have no signal of alleviation under the influence of climate change and human activities, particularly the increasing water demand due to the booming development of irrigated agriculture and livestock farming.

6. Conclusions

In summary, some important and interesting conclusions can be drawn as follows:

(1)

The Wakeby distribution can be used to describe the probabilistic behavior of the seven-day low flow regimes within the Tarim River Basin. Tawn copula provides a very good fit to the data with a NSE = 0.9978 and is selected as the best copula according to AIC, BIC, maximum likelihood, and other residual-based metrics in Yarkand River Basin. Farlie–Gumbel–Morgenstern copula and Frank copula is the best copula and should be used for further analysis in Aksu River Basin and Kaidu River Basin.

(2)

After 1987, the increase of temperature and precipitation enable the low flow volume to increase to a certain degree. However, the climate change within the Tarim River Basin is uneven in both space and time. Water supply sources in different tributary basins of the Tarim River are different. Hence, the increasing magnitude of the low flow regime is different in different tributaries. However, the massive increase of water demand due to increased irrigated agriculture and livestock farming greatly reduces the streamflow input into the main stem of the Tarim River Basin. Thus, water shortage in the lower Tarim River will not be alleviated.

(3)

Hydrological droughts of longer return periods are prone to increasing occurrence frequency. The water supply cannot satisfy the increasing water demand due to significantly increased irrigated arable land and growing population, although the precipitation is increasing and the snowmelt is also increasing due to increased temperature in recent years. The water shortage in Xinjiang is still a challenge for the sustainable development of the local socio-economy. In this case, the development of water-saving technology for irrigated agriculture and effective water resources management is necessary for the sustainable development of regional social economy and conservation of the eco-environment.