Multiple causes of nonstationarity in the Weihe annual low-flow series

Under the background of global climate change and local anthropogenic activities, multiple driving forces have introduced various nonstationary components into low-flow series. This has led to a high demand on low-flow frequency analysis that considers nonstationary conditions for modeling. In this study, through a nonstationary frequency analysis framework with the generalized linear model (GLM) to consider time-varying distribution parameters, the multiple explanatory variables were incorporated to explain the variation in low-flow distribution parameters. These variables are comprised of the three indices of human activities (HAs; i.e., population, POP; irrigation area, IAR; and gross domestic product, GDP) and the eight measuring indices of the climate and catchment conditions (i.e., total precipitation P, mean frequency of precipitation events λ, temperature T, potential evapotranspiration (EP), climate aridity index AI_EP, base-flow index (BFI), recession constant K and the recession-related aridity index AI_K). This framework was applied to model the annual minimum flow series of both Huaxian and Xianyang gauging stations in the Weihe River, China (also known as the Wei He River). The results from stepwise regression for the optimal explanatory variables show that the variables related to irrigation, recession, temperature and precipitation play an important role in modeling. Specifically, analysis of annual minimum 30-day flow in Huaxian shows that the nonstationary distribution model with any one of all explanatory variables is better than the one without explanatory variables, the nonstationary gamma distribution model with four optimal variables is the best model and AI_K is of the highest relative importance among these four variables, followed by IAR, BFI and AI_EP. We conclude that the incorporation of multiple indices related to low-flow generation permits tracing various driving forces. The established link in nonstationary analysis will be beneficial to analyze future occurrences of low-flow extremes in similar areas.

1 Introduction

Low flow is defined as the flow of water in a stream during prolonged dry weather (WMO, 2009). Yu et al. (2014) quantitatively described a low-flow event as a segment of hydrograph during a period of dry weather with discharge values below a preset (relatively small) threshold. According to WMO (2009), annual minimum flows averaged over several days can be used to measure low flows. During low-flow periods, the magnitude of river flow will greatly restrict its various functions (e.g., providing water supply for production and living, diluting waste water, ensuring navigation, meeting ecological water requirement). Therefore, the investigation of the magnitude and frequency of low flows is of primary importance for engineering design and water resources management (Smakhtin, 2001). In recent years, low flows, as an important part of river flow regime, have been attracting an increasing attention of hydrologists and ecologists in the context of the significant impacts of climate change and human activities (HAs; Bradford and Heinonen, 2008; Du et al., 2015; Kam and Sheffield, 2015; Kormos et al., 2016; Liu et al., 2015; Sadri et al., 2016). In general, under the impact of a changing environment, combinations of multiple factors, such as precipitation change, temperature change, irrigation area (IAR) change and construction of reservoirs, can drive various patterns of streamflow changes (Liu et al., 2017; Tang et al., 2016). Unfortunately, when subjected to a variety of influencing forces, low flow is more vulnerable than high flow or mean flow. Therefore, it is a pretty important issue in hydrology to identify low-flow changes, track multiple driving factors and quantify their contributions from the perspective of hydrological frequency analysis.

In hydrological analysis and design, conventional frequency analysis estimates the statistics of a hydrological time series based on recorded data with the stationary hypothesis which means that this series is “free of trends, shifts or periodicity (cyclicity)” (Salas, 1993). However, global warming and human forces have changed climate and catchment conditions in some regions. Time-varying climate and catchment conditions (TCCCs) can affect all aspects of the flow regime, i.e., changing the frequency and magnitude of floods, altering flow seasonality and modifying the characteristics of low flows. The hypothesis of stationarity has been suspected (Milly et al., 2008). If this problematic method is still used, the frequency analysis may lead to high estimation error in hydrological design. Therefore, considerable literature has introduced the concept of hydrologic nonstationarity into analysis of various hydrological variables, such as annual runoff (Arora, 2002; Jiang et al., 2015a; Xiong et al., 2014; Yang and Yang, 2013), flood (Gilroy and Mccuen, 2012; Kwon et al., 2008; Yan et al., 2017; Zhang et al., 2015), low flow (Du et al., 2015; Jiang et al., 2015b; Liu et al., 2015), precipitation (Gu et al., 2017; Mondal and Mujumdar, 2015; Villarini et al., 2010) and so on. Compared with the literature on annual runoff, floods and precipitation, the literature on the nonstationary analysis of low flow is relatively limited.

Previous hydrological literature on frequency analysis of nonstationary hydrological series mainly focuses on two aspects: development of the nonstationary method and exploration of covariates reflecting changing environments. Strupczewski et al. (2001) presented the method of time-varying moment which assumes that the hydrological variable of interest obeys a certain distribution type, but its moments change over time. The method of time-varying moment was modified to be the method of time-varying parameter values for the distribution representative of hydrologic data (Richard et al., 2002). Villarini et al. (2009) presented this method using the generalized additive models for location, scale and shape parameters (GAMLSS; Rigby and Stasinopoulos, 2005), a flexible framework to assess nonstationary time series. The time-varying parameter method can be extended to the physical covariate analysis by replacing time with any other physical covariates (Jiang et al., 2015b; Kwon et al., 2008; López and Francés, 2013; Liu et al., 2015; Villarini and Strong, 2014). For example, Jiang et al. (2015b) used reservoir index as an explanatory variable based on the time-varying copula method for bivariate frequency analysis of nonstationary low-flow series in Hanjiang River, China. Du et al. (2015) took precipitation and air temperature as the explanatory variables to explain the inter-annual variability in low flows of the Weihe River, China (also known as the Wei He River). Liu et al. (2015) took the sea surface temperature in the Nino3 region, the Pacific Decadal Oscillation, the sunspot number (3 years ahead), the winter areal temperature and precipitation as the candidate explanatory variables to explain the inter-annual variability in low flows of Yichang station, China. Kam and Sheffield (2015) ascribed the increasing inter-annual variability of low flows over the eastern United States to the North Atlantic Oscillation and Pacific North America.

To our knowledge, compared with the nonstationary flood frequency analysis, the studies on the nonstationary frequency analysis of low-flow series are not very extensive because of incomplete knowledge of low-flow generation (Smakhtin, 2001). Most of these studies explain nonstationarity of low-flow series only by using climatic indicators or a single indicator of human activity. However, the indicators of catchment conditions (e.g., recession rate) related to physical hydrological processes have seldom been attached in nonstationary modeling of low-flow series. This lack of linking with hydrological processes makes it impossible to accurately quantify the contributions of influencing factors for the nonstationarity of low-flow series, and such a scientific demand for tracing the sources of nonstationarity of low-flow series and qualifying their contributions motivated the present study. The knowledge of low-flow generation has been increased by efforts of hydrologists, which can help develop physical covariates to address nonstationarity. Low flows generally originate from groundwater or other delayed outflows (Smakhtin, 2001; Tallaksen, 1995). Their generation relates to both an extended dry weather period (leading to a climatic water deficit) and complex hydrological processes which determine how these deficits propagate through the vegetation, soil and groundwater system to streamflow (WMO, 2009). Thus, not only climate condition drivers (e.g., potential evaporation exceeds precipitation), but also catchment condition drivers (e.g., the faster hydrologic response rate to precipitation) can cause low flows.

The significant factors such as precipitation, temperature, evapotranspiration (EP), streamflow recession, large-scale teleconnections and human forces may play important roles in influencing low-flow generation (Botter et al., 2013; Giuntoli et al., 2013; Gottschalk et al., 2013; Kormos et al., 2016; Sadri et al., 2016). Gottschalk et al. (2013) presented a derived low-flow probability distribution function with climate and catchment characteristics parameters (i.e., the mean length of dry spells λ⁻¹ and recession constant of streamflow K ) as its distribution parameters. Botter et al. (2013) derived a measurable index (${\mathit{\lambda }}^{-\mathrm{1}}/K)$ which can be used for discriminating erratic river flow regimes from persistent river flow regimes. Recently, Van Loon and Laaha (2015) used climate and catchment characteristics (e.g., the duration of dry spells in precipitation and the base-flow index, BFI) to explain the duration and deficit of the hydrological drought event and offered a further understanding of low-flow generation. These studies indicated that climate and catchment conditions play an important role in producing low flows.

The goal of this study is to trace origins of nonstationarity in low flows through developing a nonstationary low-flow frequency analysis framework with the consideration of the time-varying climate and catchment conditions and human activity. In this framework, the climate and catchment conditions are quantified using the eight indices, i.e., meteorological variables (total precipitation P, mean frequency of precipitation events λ, temperature T and potential evapotranspiration), basin storage characteristics (base-flow index, recession constant K) and aridity indexes (climate aridity index AI_EP, the recession-related aridity index AI_K). The specific objectives of this study are (1) to find the most important index to explain the nonstationarity of low-flow series, (2) to determine the best subset of TCCCs indices and/or human activity indices (i.e., population, POP; irrigation area; and gross domestic product, GDP) for the final model through the stepwise selection method to identify nonstationary mode of low-flow series and (3) to quantify the contribution of selected explanatory variables to the nonstationarity.

This paper is organized as follows. Section 2 describes the methods. The Weihe River basin and available data sets used in this study are described in Sect. 3, followed by a presentation of the results and discussion in Sect. 4. Section 5 summarizes the main conclusions.

Figure 1 The framework of nonstationary low-flow frequency analysis.

2 Methodology

The flowchart of how to organize the nonstationary low-flow frequency analysis framework is shown in Fig. 1. The whole process is divided into three steps. The first step is the preliminary analysis, including the graphical presentation of both explanatory variables and low-flow series, the statistical test for nonstationarity, and the correlations between each explanatory variable and each low-flow series. The second step is the single covariate analysis for the most important explanatory variable. The third step is the multiple covariate analysis for the optimal combination. We use a low-flow frequency analysis model and stepwise regression method to accomplish the last two steps. In the following subsections, first, the low-flow frequency analysis model is constructed based on the nonstationary probability distributions method, in which distribution parameters serving as response variables can vary as functions of explanatory variables. Second, the distribution types used to build the nonstationary model are outlined. Then, the candidate explanatory variables related to the time-varying climate and catchment conditions and human activity are clarified. Finally, estimation of model parameters and selection of models are illustrated.

2.1 Construction of the low-flow nonstationary frequency analysis model

Generally, a nonstationary frequency analysis model can be established based on the time-varying distribution parameters method (Du et al., 2015; López and Francés, 2013; Liu et al., 2015; Richard et al., 2002; Villarini and Strong, 2014). For the nonstationary probability distribution f_Y(Y_t|θ^t), let Y_t be a random variable at time $t(t=\mathrm{1},\mathrm{2},\mathrm{\dots },N)$ and vector ${\mathit{\theta }}^t=[{\mathit{\theta }}_{\mathrm{1}}^t,{\mathit{\theta }}_{\mathrm{2}}^t,\mathrm{\dots },{\mathit{\theta }}_m^t]$ be the time-varying parameters. The number of parameters m in hydrological frequency analysis is generally limited to three or less. The function relationship between the kth parameter ${\mathit{\theta }}_k^t$ and the multiple explanatory variables is expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{g}_k\left({\mathit{\theta }}_k^t\right)=h_k\left(x_{\mathrm{1}}^t,x_{\mathrm{2}}^t,\mathrm{\dots },x_n^t\right),\end{array}$$

where $x_{\mathrm{1}}^t,x_{\mathrm{2}}^t,\mathrm{\dots },x_n^t$ are explanatory variables, n is the number of explanatory variables, g_k(⋅) is the link function which ensures the compliance with restrictions on the sample space and is usually set to natural logarithm for the given negative predictions and h_k(⋅) is the function for nonstationary modeling. The generalized linear model theory (GLM; Dobson and Barnett, 2012) is used to build function relationships between distribution parameters and their explanatory variables. In GLMs, the response relationship can be generally expressed as

$$\begin{array}{}\text{(2)}& {\displaystyle }{g}_k\left({\mathit{\theta }}_k^t\right)={\mathit{\alpha }}_{\mathrm{0}k}+\sum _{i=\mathrm{1}}{\mathit{\alpha }}_{ik}{x}_i^t,\end{array}$$

where ${\mathit{\alpha }}_{ik}(i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },n,k=\mathrm{1},\mathrm{\dots },m)$ are the GLM parameters.

In order to compare the nonstationary models constructed by various combinations of explanatory variables, Eq. (2) is modified in this study using the dimensionless method for the standard GLM parameters. The value of ${\mathit{\theta }}_k^t$ could be assumed to be equal to its mean (${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k)$ when all explanatory variables are equal to their mean (${\stackrel{\mathrm{‾}}{x}}_i)$, i.e.,

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\theta }}_k^t\left(x_{\mathrm{1}}^t={\stackrel{\mathrm{‾}}{x}}_{\mathrm{1}},x_{\mathrm{2}}^t={\stackrel{\mathrm{‾}}{x}}_{\mathrm{2}},\mathrm{\dots },x_n^t={\stackrel{\mathrm{‾}}{x}}_n\right)={\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k.\end{array}$$

Equation (2) is then modified as

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{l}{g}_k\left({\displaystyle \frac{{\mathit{\theta }}_k^t}{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k}}\right)={\mathit{\beta }}_{\mathrm{0}k}+\sum _{i=\mathrm{1}}^{i=n}{\mathit{\beta }}_{ik}{z}_i^t\\ z_i^t={\displaystyle \frac{x_i^t-{\stackrel{\mathrm{‾}}{x}}_i}{s_i}},i=\mathrm{1},\mathrm{2},\mathrm{\dots },n\\ {\mathit{\beta }}_{\mathrm{0}k}=g_k\left({\displaystyle \frac{{\mathit{\theta }}_k^t}{{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k}}\left|{\mathit{\theta }}_k^t={\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k\right.\right)=g_k\left(\mathrm{1}\right),\end{array}\end{array}$$

where $z_i^t$ is the normalized explanatory variable, s_i is the standard deviation of $x_i^t$ and ${\mathit{\beta }}_{ik}(i=\mathrm{1},\mathrm{2},\mathrm{\dots },n,k=\mathrm{1},\mathrm{\dots },m)$ are the standard GLM parameters. Letting the link function g_k(⋅) be the natural logarithmic function ln (⋅) and ${\mathit{\theta }}_l^t$ be the distribution parameter in $[{\mathit{\theta }}_{\mathrm{1}}^t,{\mathit{\theta }}_{\mathrm{2}}^t,\mathrm{\dots },{\mathit{\theta }}_m^t]$ with the most significant change, the degree of nonstationarity in low-flow series can be defined as $\ln \left({\mathit{\theta }}_l^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_l\right)$. Then, the contribution $c_i^t$ of each explanatory variable $x_i^t$ to $\ln \left({\mathit{\theta }}_l^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_l\right)$ could be defined as

$$\begin{array}{}\text{(5)}& {\displaystyle }{c}_i^t={\mathit{\beta }}_{il}{\displaystyle \frac{x_i^t-{\stackrel{\mathrm{‾}}{x}}_i}{s_i}}.\end{array}$$

Table 1 The probability density functions and moments (the mean and variance) for the candidate distributions in this study.

2.2 Candidate distribution functions

We need to select the form of probability distribution f_Y(⋅) to determine what type of nonstationary frequency curves will be produced. Various probability distributions have been compared or suggested in modeling of low-flow series (Du et al., 2015; Hewa et al., 2007; Liu et al., 2015; Matalas, 1963; Smakhtin, 2001). An extensive overview of distribution functions for low flow is given in Tallaksen et al. (2004). Following these recommendations, we consider five distributions, i.e., Pearson type III (PIII), gamma (GA), Weibull (WEI), lognormal (LOGNO) and generalized extreme value (GEV) as candidates in this study (Table 1). In the case of Pearson type III distribution, considering that the parameter θ₃ of Pearson type III as lower bound should approach zero and the parameter θ₃ of GEV is quite sensitive and difficult to be estimated, we assume them to be constant in this study.

Table 2 Description of the developed nonstationary models using time, TCCCs indices and/or HA indices as explanatory variables.

2.3 Candidate explanatory variables

We look for variables $x_{\mathrm{1}}^t,x_{\mathrm{2}}^t,\mathrm{\dots },x_n^t$ that can explain parts of the variations in distribution parameters θ^t. From the perspective of low-flow generation, the dependency between low-flow regime and both climate and catchment conditions has been presented by previous studies (Botter et al., 2013; Gottschalk et al., 2013; Van Loon and Laaha, 2015). We focus on eight measuring indices: precipitation, mean frequency of precipitation events, temperature, potential evapotranspiration, climate aridity index, base-flow index, recession constant and recession-related aridity index. These indices were chosen to incorporate time-varying climate and catchment conditions in nonstationary modeling of low-flow frequency and serve as candidate explanatory variables. Climate variables (i.e., precipitation, mean frequency of precipitation events, temperature, potential evapotranspiration and climate aridity index) are related to both water supply source and water loss and are therefore selected as candidate variables. It has been shown that the base-flow index and recession constant reflect the storage and release capability of the catchments (Van Loon and Laaha, 2015). The recession-related aridity index reflects both the water supply and storage capability (Botter et al., 2013). In addition to TCCCs indices, the three indices of human activity (irrigation area, population and gross domestic product) are related to water withdrawal loss for agricultural, domestic and industrial purposes and are therefore included. The detailed reasons for selecting all indices are summarized in Table 3. The values of them at each year could be estimated from hydro-meteorological data and human activity data. Annual precipitation (P) and temperature (T) are calculated directly by meteorological data. The remaining TCCCs indices need to be estimated indirectly. Detailed estimation procedures are shown in the following subsections.

2.3.1 Annual mean frequency of precipitation events (λ)

Annual mean frequency of precipitation events is defined as an index to represent the intensity of precipitation recharge to the streamflow:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\lambda }={\displaystyle \frac{\mathrm{1}}{W}}\sum _{w=\mathrm{1}}^{w=W}{\displaystyle \frac{N_w\left(A\right)}{t_r}},\end{array}$$

where N_w(A) is the number of daily rainfall events A (with values more than the threshold 0.5 mm) in wth windows with a length t_r; W is the number of windows.

2.3.2 Annual climate aridity index (AI_EP)

The ratio of annual potential evaporation to precipitation, commonly known as the climate aridity index, has been used to assess the impacts of climate change on annual runoff (Arora, 2002; Jiang et al., 2015a). The climate aridity index largely reflects the climatic regimes in a region and determines runoff rates (Arora, 2002). Therefore, we choose the annual climate aridity index as a measure of time-varying climate and catchment conditions and estimate its value in a whole region using

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathrm{AI}}_{\mathrm{EP}}={\displaystyle \frac{\mathrm{EP}}{P}},\end{array}$$

where P is annual areal precipitation (mm) and EP is annual areal potential evapotranspiration (mm). The Hargreaves equation (Hargreaves and Samani, 1985) is applied to calculate EP using the R package “Evapotranspiration” (Guo, 2014).

2.3.3 Annual base-flow index (BFI)

The base-flow index (BFI) is defined as the ratio of base flow to total flow. This index has been applied to quantify catchment conditions (e.g., soil, geology and storage-related descriptors) to explain hydrological drought severity (Van Loon and Laaha, 2015). We also choose annual base-flow index as a measure of TCCCs. BFI is estimated using a hydrograph separation procedure in the R package “lfstat” (Koffler and Laaha, 2013).

2.3.4 Annual streamflow recession constant (K)

The recession constant is an important catchment characteristic index measuring the timescale of the hydrological response and reflecting water retention ability in the upstream catchment (Botter et al., 2013). Various estimation methods have been developed to extract recession segments and to parameterize characteristic recession behavior of a catchment (Hall, 1968; Sawaske and Freyberg, 2014; Tallaksen, 1995).

In this study, annual recession analysis (ARA) is performed to obtain the annual streamflow recession constant (K). In ARA, the linearized Dupuit–Boussinesq equation is used to parameterize characteristic recession behavior of a catchment and is written as

$$\begin{array}{}\text{(8)}& {\displaystyle }-{\displaystyle \frac{\mathrm{d}{Q}_t}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{K}}{Q}_t,\end{array}$$

where Q_t is the value at time t. Equation (8) is investigated by plotting data points $\frac{\mathrm{d}{Q}_t}{\mathrm{d}t}$ against Q_t of all extracted recession segments from hydrographs at each year. The criteria of recession segment extraction are based on the Manual on Low-flow Estimation and Prediction (WMO, 2009). Then, the annual recession rate (K⁻¹) is estimated as the slope of the fitted straight line of these data points with the least-squares method. We calculated K using the R package “lfstat” (Koffler and Laaha, 2013).

2.3.5 Annual recession-related aridity index (AI_K)

In this study, the recession-related aridity index is defined as the ratio of recession rate (K⁻¹) to mean precipitation frequency (λ), denoted as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{AI}}_K={\displaystyle \frac{K^{-\mathrm{1}}}{\mathit{\lambda }}}.\end{array}$$

This ratio plays an important role in controlling the river flow regime (Botter et al., 2013; Gottschalk et al., 2013) and serves as an indicator measuring the recession-related aridity degree of the streamflow in the river channel. For example, the faster recession process or lower precipitation frequency may lead to increased runoff loss or decreased precipitation supply. Consequently, the higher the value AI_K is, the more likely low-flow events occur, and vice versa.

2.4 Parameter estimation

The model parameters including ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k(k=\mathrm{1},\mathrm{2},\mathrm{\dots },m)$ and ${\mathit{\beta }}_{ik}(i=\mathrm{1},\mathrm{2},\mathrm{\dots },n,k=\mathrm{1},\mathrm{\dots },m)$ need to be estimated. ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_k(k=\mathrm{1},\mathrm{2},\mathrm{\dots },m)$ are estimated from outputs of stationary frequency analysis through the maximum likelihood method. We have

$$\begin{array}{}\text{(10)}& {\displaystyle }L\left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}},\mathrm{\dots },{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_m\right)=\sum _{t=\mathrm{1}}^{t=N}\ln \left[f_Y\left(y_t\left|{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{2}},\mathrm{\dots },{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_m\right.\right)\right],\end{array}$$

where y_t is observed low flow at time t and N is the number of samples. The parameters ${\mathit{\beta }}_{ik}(i=\mathrm{1},\mathrm{2},\mathrm{\dots },n,k=\mathrm{1},\mathrm{\dots },m)$ are estimated through the maximum likelihood method to produce nonstationary low-flow frequency curves:

$$\begin{array}{ll}\text{(11)}& {\displaystyle }& {\displaystyle }L\left(\begin{array}{l}{\mathit{\beta }}_{\mathrm{11}},\mathrm{\dots },{\mathit{\beta }}_{n\mathrm{1}}\\ \mathrm{\dots }\\ {\mathit{\beta }}_{\mathrm{1}m},\mathrm{\dots },{\mathit{\beta }}_{nm}\end{array}\right)=\sum _{t=\mathrm{1}}^{t=N}\ln {\displaystyle }& {\displaystyle }\left\{f_Y\left(y_t\left|{\mathit{\theta }}_{\mathrm{1}}^t\left(z_{\mathrm{1}}^t,\mathrm{\dots },z_n^t\left|{\mathit{\beta }}_{\mathrm{11}},\mathrm{\dots },{\mathit{\beta }}_{n\mathrm{1}}\right.\right),\mathrm{\dots },{\mathit{\theta }}_m^t\left(z_{\mathrm{1}}^t,\mathrm{\dots },z_n^t\left|{\mathit{\beta }}_{\mathrm{1}m},\mathrm{\dots },{\mathit{\beta }}_{nm}\right.\right)\right.\right)\right\}.\end{array}$$

The residuals (normalized randomized quintile residuals) are used to test the goodness of fit of fitted model objects (Dunn and Symth, 1996):

$$\begin{array}{}\text{(12)}& {\displaystyle }{\widehat{r}}_t={\mathrm{\Phi }}^{-\mathrm{1}}\left(F_Y\left(y_t\left|{\widehat{\mathit{\theta }}}^t\right.\right)\right),\end{array}$$

where F_Y(⋅) is the cumulative distribution of y_t and ${\mathrm{\Phi }}^{-\mathrm{1}}\left(\cdot \right)$ is the inverse function of the standard normal distribution. The distribution of the true residuals ${\widehat{r}}_t$ converges to standard normal if the fitted model is correct. A worm plot (Buuren and Fredriks, 2001) is used to check whether ${\widehat{r}}_t$ have a standard normal distribution.

Figure 2 Location, topography, hydro-meteorological stations and river systems of the Weihe River basin.

Table 3 The summary of candidate explanatory variables and reason of selection.

Figure 3 The annual minimum low flows and fitted trend lines in both Huaxian (H) and Xianyang (X) gauging stations.

2.5 Model selection

Model selection contains the selection of the type of probability distribution and the selection of the explanatory variables to explain the response variables (i.e., distribution parameters θ₁ and θ₂). In order to obtain the final optimal model, the selection of the explanatory variables for θ₁ and θ₂ is conducted by stepwise selection strategies (Stasinopoulos and Rigby, 2007; Venables, 2002): i.e., select a best subset of candidate explanatory variables for θ₁ using a forward approach (which starts with no explanatory variable in the model and tests the addition of each explanatory variable using a chosen model fit criterion); given this subset for θ₁ select another subset for θ₂ (forward). The stepwise selection strategies can get a series of stepwise models with different numbers of explanatory variables, as shown in Fig. 1. In order to detect how the number of explanatory variables influences the performance of the model for describing nonstationarity, we investigate the eight types of stepwise models as shown in Table 2: the zero-covariate model or stationary model (M0), the time covariate model (M1), the single physical covariate model M2 (single TCCCs covariate model M2a or single HA covariate model M2b), two TCCCs covariates model (M3), the optimal TCCCs covariates model (M4), the optimal HA covariates model (M5) and the final model (M6). The model fit criterion is based on the Akaike's information criterion (AIC; Akaike, 1974) as shown by the following

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{AIC}=-\mathrm{2}\mathrm{ML}+\mathrm{2}\mathrm{df},\end{array}$$

where ML is the log-likelihood in Eq. (11) and df is the number of degrees of freedom. The model with the lower AIC value was considered better.

3 Study area and data

3.1 The study area

The Weihe River, located in the southeast of the northwest Loess Plateau, is the largest tributary of the Yellow River, China. The Weihe River has a drainage area of 134 766 km², covering the coordinates of $\mathrm{33}{}^{\circ }{\mathrm{42}}^{\prime }$–$\mathrm{37}{}^{\circ }{\mathrm{20}}^{\prime }$ N, $\mathrm{104}{}^{\circ }{\mathrm{18}}^{\prime }$–$\mathrm{110}{}^{\circ }{\mathrm{37}}^{\prime }$ E (Fig. 2). This catchment generally has a semi-arid climate, with extensive continental monsoonal influence. Average annual precipitation of the whole area over the period 1954–2009 is about 540 mm and has a wide range (400–1000 mm) in various regions. Under the significant impacts of climate change and human activities in the Weihe River basin in recent decades, the hydrological regime of the river has changed over time (Du et al., 2015; Jiang et al., 2015a; Xiong et al., 2015).

In the Weihe basin, the impacts of agricultural irrigation on runoff have been found to be significant (Jiang et al., 2015a; Lin et al., 2012). Lin et al. (2012) mentioned that the annual runoff of the Weihe River was significantly affected by irrigation diversion of the Baoji Gorge irrigation district. The irrigated area of the Baoji Gorge irrigation district increased over time since the founding of P.R. China in 1949, and, due to one influential irrigation system project in that area, it became more than twice as large as the original irrigation area since 1971. Jiang et al. (2015a) demonstrated that, in the Weihe basin, irrigated area, as compared with the other indices, e.g., population, gross domestic product and cultivated land area, was a more suitable human explanatory variable for explaining the time-varying behavior of annual runoff. With the above background, it is important to consider the effects of human activities that mainly originate from irrigation diversion, especially for studying low-flow series in this basin. The estimations of annual recession rate (K⁻¹) by the daily streamflow data are expected to incorporate the information of impacts of water diversions on the low flows in the river channel.

Figure 4 Frequency distributions (using the kernel density estimations) and time series processes of TCCCs variables in both Huaxian (H) and Xianyang (X) stations.

3.2 Data

We used daily streamflow records (1954–2009) provided by the Hydrology Bureau of the Yellow River Conservancy Commission from both Huaxian station (with a drainage area of 106 500 km²) and Xianyang station (with a drainage area of 46 480 km²). Low-flow extreme events were selected from the daily streamflow series using the widely used annual minimum series method (WMO, 2009). AM_n is the annual minimum n-day flow during hydrological year beginning on 1 March. Consequently, AM₁, AM₇, AM₁₅ and AM₃₀ are selected as low-flow extreme events in this study. The original measure unit of streamflow data (m³ s⁻¹) is converted to 10⁻⁴ m³ s⁻¹ km⁻² for convenience of comparison of results between the Huaxian and Xianyang gauging stations

We downloaded daily total precipitation and daily mean air temperature records for 19 meteorological stations over the basin from the National Climate Center of the China Meteorological Administration (source: http://www.cma.gov.cn/). The areal average daily series of both variables above Huaxian and Xianyang stations are calculated using the Thiessen polygon method (Szolgayova et al., 2014; Thiessen, 1911). The annual average temperature (T) and annual total precipitation (P) over the period 1954–2009 are calculated for each catchment.

Human activity data (i.e., gross domestic product, population and irrigation area) were taken from annals of statistics provided by the Shaanxi Provincial Bureau of Statistics (http://www.shaanxitj.gov.cn/) and Gansu Provincial Bureau of Statistics (source: http://www.gstj.gov.cn/).

Figure 5 HA indices in both Huaxian (H) and Xianyang (X). (a), (b) and (c) are for gross domestic production (GDP), population (POP) and irrigated area (IAR), respectively.

4 Results and discussion

4.1 Identification of nonstationarity

The graphical representation and statistical test provide a preliminary analysis for low-flow nonstationarity. The graphical representations of time-series data help visualize the trends of related variables (i.e., low flow, TCCCs and HA variables), the density distributions of TCCCs variables, and the correlations between low-flow variables and these explanatory variables. In Fig. 3, four annual minimum streamflow series (AM₁, AM₇, AM₁₅ and AM₃₀) in both Huaxian and Xianyang gauging stations show overall decreasing trends, as indicated by the fitted (dashed) trend lines. Compared with Huaxian, Xianyang has a larger runoff modulus (the flow per square kilometer) and a larger decrease in annual minimum streamflow series. For example, the decline slope of AM₃₀ is −0.0725 (10⁻⁴ m³ s⁻¹ km⁻² yr⁻¹) in Huaxian station while Xianyang station it is −0.1338 (10⁻⁴ m³ s⁻¹ km⁻² yr⁻¹).

Figure 4 shows the kernel density estimations and time processes of TCCCs variables for both Huaxian (H) and Xianyang (X) stations. The results show that these variables have different variation patterns. For example, the mean frequency of precipitation events (λ) has a decreasing trend, while temperature (T) has an increasing trend. As presented by Fig. 5, three HA variables have a significant upward trend, especially the irrigation area which is increased greatly after about 1970, suggesting that the impact of human activities in this basin has increased over time.

Table 4 The results of trend test and change-point detection for both the four low-flow series and TCCCs variables in Huaxian and Xianyang.

Significance codes: 0 “***” 0.001 “**” 0.01 “*” 0.05 “•” 0.1 “ ” 1

Figure 6 The Pearson correlation coefficients matrix between the annual minimum flow series and candidate explanatory variables in Huaxian (H) and Xianyang (X) stations; the darker color intensity represents a higher level of correlation (blue indicates positive correlation and red indicates negative correlations).

Figure 7 Comparisons among M0, M1 and M2 based on the AIC values for the four observed low-flow series in Huaxian (H) at (a) and Xianyang (X) at (b); darker red color represents a higher goodness of fit.

The significance of trends in the four annual minimum streamflow series and TCCCs variables is tested by the Mann–Kendall trend test (Kendall, 1975; Mann, 1945; Yue et al., 2002), and the change points in these series are detected by Pettitt's test (Pettitt, 1979). The results in Table 4 show that, in both Huaxian and Xianyang stations, the decreasing trends in all the four low-flow series (AM₁, AM₇, AM₁₅ and AM₃₀) and two explanatory variables (λ and P), as well as the increasing trends in T, EP, and AI_EP are significant at the 0.05 level (Table 4), but BFI shows no significant trends. However, K and AI_K had significantly decreasing trends only in Huaxian station (p value <0.05). The results of change-point detection show that all low-flow series are located at 1968–1971 (p value <0.05) except AM₃₀ at Xianyang station whose change point is located at 1993 (p value <0.05); for the eight candidate explanatory variables, the change points of the variables related to temperature (T, EP, AI_EP) in both stations are located at 1990–1993 (p value <0.05), the change points of the variables related to precipitation (λ, P) in both stations are close at 1984–1990 (p value ≤0.186) and the change points of the variables related to streamflow recession (K, AI_K) in Huaxian station are located at 1968–1971 (p value <0.05). However, BFI in both stations and K and AI_K in Xianyang station show no significant change points.

A preliminary attribution analysis is performed using the Pearson correlation matrix to investigate the relations between the annual minimum series and eight candidate explanatory variables. Figure 6 indicates that there are significant linear correlations between the four minimum low-flow series (AM₁, AM₇, AM₁₅ and AM₃₀) and all the explanatory variables except GDP have the absolute values of Pearson correlation coefficients larger than 0.27 (p value <0.05). These potential physical causes of nonstationarity in low flows are further considered by establishing the low-flow nonstationary model with TCCCs and HA variables in the following section.

4.2 Nonstationary frequency analysis models

4.2.1 Single covariate models

Figure 7 presents the AIC values of the four types of models (M0, M1, M2a and M2b) fitted for the low-flow series (AM₁, AM₇, AM₁₅ and AM₃₀). Some interesting results are shown as follows. First, nonstationary models (M1, M2a and M2b) have lower AIC values than stationary model (M0), which suggests that nonstationary models are worth considering. Second, for Huaxian station, irrespective of the chosen explanatory variables, the distribution type plays an important role in modeling nonstationary low-flow series. For example, PIII, GA and WEI distributions in AM₁₅ and AM₃₀ cases have lower AIC values than LOGNO and GEV distributions. However, for Xianyang, choosing a suitable explanatory variable may be more important than choosing a distribution type. For example, variables t, P, T, AI_EP, POP and IAR in most cases have lower AIC values than the other explanatory variables. Finally, in Huaxian, the lowest AIC values for modeling AM₁, AM₇, AM₁₅ and AM₃₀ are found in GEV_M2b_IAR, LOGNO_M2b_IAR, PIII_M2a_AI_K and GA_M2a_AI_K, respectively, while in Xianyang the lowest AIC values for modeling AM₁, AM₇, AM₁₅ and AM₃₀ are found in GEV_M2b_IAR, GEV_M2b_IAR, PIII_M2b_IAR and GEV_M2b_IAR, respectively. These results indicated that for explaining nonstationarity of low flow in Huaxian station, IAR is the most dominant HA variable and AI_K is the most dominant TCCCs variable, while in Xianyang the most dominant HA variable is IAR and the most dominant TCCCs variables causing nonstationarity in AM₁, AM₇, AM₁₅ and AM₃₀ are K, AI_EP, AI_EP and T, respectively.

Figure 8 Performance assessments of GA_M2 for AM₃₀ in Huaxian (H) on the left panel and Xianyang (X) on the right panel. (a) and (b) are the centile curves plots of GA_M2 (red lines represent the centile curves estimated by GA_M2; the 50th centile curves are indicated by thick red; the yellow-filled areas are between the 5th and 95th centile curves; the black points indicate the observed series); (c) and (d) are the worm plots of GA_M2 for the goodness-of-fit test; a reasonable model fit should have the data points fall within the 95 % confidence intervals (between the two red dashed curves).

Figure 8 shows the diagnostic assessment of the GA_M2 model (with the optimal explanatory variable) for AM₃₀ in both Huaxian and Xianyang stations. The centile curve plots of GA_M2 (Fig. 8a and b) show the observed values of AM₃₀, the estimated median and the areas between the 5th and 95th centile. Figure 8a shows the response relationship between AM₃₀ and AI_K in Huaxian: the increase in AI_K means the smaller magnitude of low-flow events because a high value of AI_K (faster stream recession or fewer rainy days) may lead to faster water loss or less supply. In Fig. 8b, the higher values of IAR means the smaller magnitude of low-flow events, which suggests that IAR plays an important role in driving low-flow generation in Xianyang. Figure 8c and d show that the worm points are within the 95 % confidence intervals, thereby indicating a good model fit and a reasonable model construction.

Figure 9 Comparisons of performance of stationary model (M0), time covariate model (M1) and physical covariate models (M2a, M2b, M3, M4, M5 and M6 with their corresponding optimal explanatory variables) for AM₃₀ in Huaxian (H) at (a) and Xianyang (X) at (b).

4.2.2 Multiple covariate models

Figure 9 shows the AIC values of the stationary model (M0), time covariate model (M1) and physical covariate models (M2a, M2b, M3, M4, M5 and M6) for AM₃₀. As shown in Fig. 9, M4 (nonstationary GA distribution with the optimal TCCCs variables) has a good performance; after adding the HA variables, M6 with the lowest AIC value is attained; it can be found that the combination of multiple TCCCs variables plays a major role in changing the low flows of the Weihe River, but the influence of HA variables should not be ignored.

Table 5 The summary of frequency analysis using GA distribution for AM₃₀ in Huaxian and Xianyang.

Figure 10 Performance assessments of GA_M6 for AM₃₀ in Huaxian (H) on the left panel and Xianyang (X) on the right panel. (a) and (b) are the centile curves plots of GA_M6 (red lines represent the centile curves estimated by GA_M6; the 50th centile curves are indicated by thick red; the yellow-filled areas are between the 5th and 95th centile curves; the filled black points indicate the observed series); (c) and (d) are the worm plots of GA_M6 for the goodness-of-fit test; a reasonable model fit should have the data points fall within the 95 % confidence intervals (between the two red dashed curves).

Figure 11 Contribution of selected explanatory variables to $c_i^t=\ln \left({\mathit{\theta }}_{\mathrm{1}}^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}\right)$ in different periods based on GA_M6.

A summary of frequency analysis based on nonstationary GA distribution AM₃₀ is presented in Table 5. We choose to focus on M4, M5 and M6. When only using TCCCs variables to model nonstationary low-flow frequency distribution, the results of M4 show the optimal combination of explanatory variables for all low-flow series contains more than three variables. For example, for AM₃₀ of Huaxian, the optimal combination of TCCCs variables includes AI_K, BFI and AI_EP. When only HA variables are used, the results of M5 show IAR is important to the low flows in this area. And M4 has a better performance than M5. When using both TCCCs variables and HA variables, the results of M6 show the optimal combination contains multiple TCCCs variables and the irrigation area. For Huaxian, the optimal combination of all explanatory variables is AI_K, IAR, BFI and AI_EP, while for Xianyang, the optimal combination is IAR, AI_EP and BFI. We can also find that if two TCCCs variables are highly correlated, they do not seem to be selected as the explanatory variables at the same time. For example, in terms of air temperature (T), evapotranspiration and the climate aridity index (AI_EP), only one of them will appear in the optimal combination. This suggests that multicollinearity problem in the multiple variables analysis can be reduced, which will help obtain more reliable GLM parameters for contribution analysis.

The diagnostic assessment of the GA_M6 model for AM₃₀ at two stations is presented by Fig. 10. The centile curve plots of GA_M6 (Fig. 10a and b) show a more sophisticated nonstationary modeling than GA_M2 (Fig. 8). When using GA_M6 to model AM₃₀ in Huaxian (Fig. 10a), similar to GA_M2, the lower low flows are found to also correspond to a higher value of AI_K, but GA_M6 is able to identify the more complex variation patterns of low flows through the incorporation of IAR, BFI and AI_EP. Figure 10c and d show that the data points of worm plots of GA_M6 are almost within the 95 % confidence intervals, thereby indicating an acceptable model fit and a reasonable model construction.

Figure 11 presents the contribution of each selected explanatory variable to $\ln \left({\mathit{\theta }}_{\mathrm{1}}^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}\right)$ in the observation year based on GA_M6 for AM₃₀ in Huaxian and Xianyang. We can find that for Huaxian, the simulation value of $\ln \left({\mathit{\theta }}_{\mathrm{1}}^t\right)$ frequently occurs below $\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}\right)$ during the two periods of about 1970–1982 and 1993–2003, which is in accordance with the observed decrease in AM₃₀ of Huaxian station during these periods. In the former period 1970–1982, both AI_K and BFI contribute a lot of negative amount to $\ln \left({\mathit{\theta }}_{\mathrm{1}}^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}\right)$, whereas, during 1993–2003, the contribution of both AI_K and BFI decreases significantly. However, IAR has almost equal negative contribution to $\ln \left({\mathit{\theta }}_{\mathrm{1}}^t\right)-\ln \left({\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{1}}\right)$ in both periods. Unlike the former three variables, the significant negative contribution of AI_EP is only found in 1993–2003. For AM₃₀ of Xianyang, the contribution of IAR, AI_EP and BFI is similar to that at Huaxian station in two periods; however, AI_K is not included in the final model.

4.3 Discussion

The impacts of both human activities and climate change on low flows of the study area led to time-varying climate and catchment conditions. Nonstationary modeling for annual low-flow series using TCCCs variables and/or HA variables as explanatory variables is clearly different from either the stationary model (M0) or the time covariate model (M1). The result demonstrates that considering multiple drivers (e.g., the variability in catchment conditions), especially in such an artificially influenced river, is necessary for nonstationary modeling of annual low-flow series.

In this study area, nonstationary modeling considering TCCCs is supported by the following facts and findings. For human activities, an important milestone representative is the completion and operation of the irrigation system on the plateau in the Baoji Gorge irrigation district since 1971 (Sect. 3.1). Figure 5c shows the change in irrigation area in this basin. And the change-point detection test in Sect. 4.1 shows that significant change points of low-flow series occur exactly at around 1971. This result demonstrates that changes in AM₃₀ may involve a consequence of this project. In addition to human activities, climate change also makes a considerable contribution to nonstationarity of low flows, as suggested by nonstationary modeling using TCCCs variables with stepwise analysis. Actually, climate driving pattern may strengthen after nearly 1990, which is indicated by change-point detection test of both annual mean temperature (T) and annual precipitation (P) as well as the behavior of annual low-flow series after nearly 1990. There are two faster recession periods, the 1970s and the 1990s, as shown in Fig. 4. The reasons for the faster recession are likely to be related to the above-mentioned project (e.g., the increasing diversion for irrigation) and climate change (e.g., the intensified evaporation) but also could be human alterations on catchment properties, such as vegetation cover change. In conclusion, the temporal variability in irrigation area, air temperature, precipitation (the frequency and volume of rain events) and streamflow recession should be the main driving factors of generating low-flow regimes in this basin. Overall, the causes of nonstationarity in the category for two gauging stations have no clear difference but have some differences in the relative importance. As shown in Table 5, when modeling the low-flow series of Huaxian using TCCCs variables, the optimal model (M4) preferred the variables that are related to recession process; however, for Xianyang, the preferred variables are related to temperature. The reason for this may be that, as a downstream station, Huaxian station suffers more intensive human activity, so that the importance of temperature change to the low-flow change is reduced meanwhile the importance of streamflow recession (related to the capability of water storage) change is enhanced. Ignoring the negative impacts of the errors in estimating annual recession constants (K) which are caused by insufficient data points of extracted stream segments at some wet years may lead to the propagation of high errors in annual recession analysis and accordingly affect the quality of nonstationary frequency analysis when K is used as an explanatory variable. Further study will give a more reliable estimation of K through the improvement of annual recession analysis. In addition, it should be noted that the population recorded in the annals of statistics may not be equal to the actual population living in the catchment. If the population in the annals is used as the explanatory variable, this difference may lead to the uncertainty of model parameter estimations. Nonetheless, it is the best population data so far and the explanatory variable POP is excluded in the final model (M6).

5 Conclusion

There is an increasing need to develop an effective nonstationary low-flow frequency model to deal with nonstationarities caused by climate change and time-varying anthropogenic activities. In this study, time-varying climate and catchment conditions in the Weihe River basin were measured by annual time series of the eight indices, i.e., total precipitation (P), mean frequency of precipitation events (λ), temperature (T), potential evapotranspiration, climate aridity index (AI_EP), base-flow index, recession constant (K) and the recession-related aridity index (AI_K). The nonstationary distribution model was developed using these eight TCCCs indices and/or three HA indices as candidate explanatory variables for frequency analysis of time-varying annual low-flow series caused by multiple drivers. The main driving forces of the decrease in low flows in the Weihe River include reduced precipitation, warming climate, increasing irrigation area and faster streamflow recession. Therefore, a complex deterioration mechanism resulting from these factors demonstrates that, in this arid and semi-arid area, the water resources could be vulnerable to adverse environmental changes, thus portending increasing water shortages. The nonstationary low-flow model considering TCCCs can provide the knowledge of low-flow generation mechanism and give a more reliable design of low flows for infrastructure and water supply.