Nonstationary Analyses of the Maximum and Minimum Streamﬂow in Tamsui River Basin, Taiwan

: This study aims to detect non-stationarity of the maximum and minimum streamﬂow regime in Tamsui River basin, northern Taiwan. Seven streamﬂow gauge stations, with at least 27-year daily records, are used to characterize annual maximum 1- and 2-day ﬂows and annual minimum 1-, 7-, and 30-day ﬂows. The generalized additive models for location, scale, and shape (GAMLSS) are used to dynamically detect evolution of probability distributions of the maximum and minimum ﬂow indices with time. Results of time-covariate models indicate that stationarity is only noted in the 4 maximum ﬂow indices out of 35 indices. This phenomenon indicates that the minimum ﬂow indices are vulnerable to changing environments. A 16-category distributional-change scheme is employed to classify distributional changes of ﬂow indices. A probabilistic distribution with complex variations of mean and variance is prevalent in the Tamsui River basin since approximate one third of ﬂow indices (34.3%) belong to this category. To evaluate impacts of dams on streamﬂow regime, a dimensionless index called the reservoir index (RI) serves as an alternative covariate to model nonstationary probability distribution. Results of RI-covariate models indicate that 7 out of 15 ﬂow indices are independent of RI and 80% of the best-ﬁtted RI-covariate models are generally worse than the time-covariate models. This fact reveals that the dam is not the only factor in altering the streamﬂow regime in the Tamsui River, which is a signiﬁcant alteration, especially the minimum ﬂow indices. The obtained distributional changes of ﬂow indices clearly indicate changes in probability distributions with time. Non-stationarity in the Tamsui River is induced by climate change and complex anthropogenic interferences.


Introduction
Probabilistic theories have played an essential role in hydraulic facilities planning and design due to the inherent uncertainties involved in hydro-climatic series. It has become a standard practice to employ frequency analysis of extreme hydro-climatic variables in the risk evaluation of hydraulic facilities [1,2]. An important assumption in such frequency analysis-based approaches is stationarity; that is, the probability distribution used to fit the variable of interest is time invariant. However, climate change and/or anthropogenic interference-induced changing environments may lead to alterations in statistical characteristics of hydro-climatic variables, and thus, result in non-stationarity. Non-stationarity in hydro-climate series renders estimates of return period and risk used for hydraulic facilities planning and design ambiguous and questionable [3][4][5][6][7].
Focusing on Taiwan, evidence of altering rainfall regime induced by climate change was indicated by several recent studies [32][33][34][35][36][37][38]. However, little research has examined the changes in streamflow regime in Taiwan, with the exception of Yeh et al. [39] who investigated the long-term streamflow trend in northern Taiwan using the Mann-Kendall test. Nonstationary modeling of streamflow in Taiwan is rarely addressed in the literature.
The main aim of this work is to detect the presence of non-stationarity in the streamflow regime, characterized by the maximum and minimum flow indices at various time scales, of the Tamsui River located in northern Taiwan. The generalized additive models for location, scale, and shape (GAMLSS), developed by Rigby and Stasinopoulos [40], are adopted for detecting non-stationarity in streamflow characteristics by respectively incorporating time-dependent parameters and anthropogenic-impact parameters in terms of reservoir index into probability models. In addition, the 16-category distributional-change scheme, proposed by Shiau and Wu [41], is also used to classify variations of the obtained nonstationary distributions of the flow indices for distinguishing spatial variability. The obtained information may usefully guide future planning and design of water-resources engineering, adapting to the changing environment.

The Generalized Additive Models for Location, Scale, and Shape (GAMLSS)
A widely used approach to model nonstationary hydro-climatic series is the timevarying moments method, indicated by Khaleq et al. [42], which incorporates time-varying parameters into probability models with the same form of stationary condition. The GAMLSS is a popular tool to achieve this purpose in hydrology and dynamically detects evolution of probability distributions with time or other covariates [43][44][45][46].
In GAMLSS, the independent observations y i (i = 1, 2, . . . , n) are assumed to have a probability density function (PDF) f (y i |θ i ), where θ i = θ i 1 , . . . , θ i p is a parametric vector accounting for location, scale, and shape. p denotes the number of parameters and is usually less than or equal to four. The distribution parameters are related to the explanatory variables (covariates) by the monotonic link functions g k (k = 1, . . . , p), which is expressed as where g k denotes the monotonic link function for the kth parameter (here only the identity and logarithm link functions are considered); X k represents an n × m matrix of explanatory variables (covariates); β T k = (β 1k , β 2k , . . . , β mk ) is a parameter vector of length m; and h jk represents the functional dependence of the distribution parameters on explanatory variables x jk .
In this study, a semi-parametric additive formulation is adopted to model nonstationarity in flow indices in the Tamsui River basin, Taiwan; that is, functional dependence is linear or smooth and smooth dependence is based on the cubic spline function in this study. Further details of GAMLSS can be found in Rigby and Stasinopoulos [40].
A total of five widely used two-parameter distributions, including lognormal (LNO), logistic (LO), gamma (GA), Gumbel (GU), and Weibull (WEI), are employed to model the flow indices. The PDFs associated with types of link function of these distributions are summarized in Table 1. The mean and variance of these distributions in terms of parameters (θ 1 and θ 2 ) are reported in Table 2. The parameters of the probability distributions are modeled as functions of the time or reservoir index, which is described in the subsequent subsection. To avoid model over-fitting, the minimum Akaike Information Criterion (AIC) is adopted in this study to select a parsimonious model. Calculations are implemented by a free-access R-based package gamlss [47]. Table 1. Summaries of probability density function (PDF) of five distributions and the link functions used to model flow indices.

Mean Variance Note
Lognormal

Reservoir Index
To evaluate impacts of dams on the streamflow regime, a dimensionless index called the reservoir index or check dam index has been proposed in the literature [24,[48][49][50]. The dimensionless reservoir index developed by Jiang et al. [50] is adopted in this study to serve as an alternative covariate of nonstationary probabilistic models, which is defined as where N denotes the number of dams upstream of the streamflow gauge station; A i denotes the catchment area of each reservoir; A T denotes the catchment area of the gauge station; V i denotes the capacity of each reservoir; and V T denotes the total capacity of all dams upstream of the gauge station.

Overview of the Tamsui River Basin
Tamsui River, located in northern Taiwan, has a total length of 159 km and a catchment area of 2726 km 2 . Originating in Pintain Mountain, the Dahan River (main tributary) wanders northward and northwestward and receives streamflow from the Xindian River and the Jilong River before emptying into the Taiwan Strait ( Figure 1). The Tamsui River, flowing through dense population areas in northern Taiwan (Taipei City and New Taipei City), led to the construction of hydraulic facilities, initiated in the early twentieth century. Currently there are 12 dams within this basin offering domestic water-supply, flood-mitigation, and hydropower-generation purposes. Figure 1 shows the spatial locations of these 12 dams and seven selected streamflow gauge stations to reflect natural and damaltered sites. Figure 2 shows the schematic layout of the Tamsui River system including tributaries, streamflow gauge stations, and dams, which clearly demonstrates the upstream and downstream locations of dams and stations.  According to Shiau and Wu [41], the mean annual total rainfall at Taipei station, with the longest records in the Tamsui River basin, is 2169 mm for the period 1897−2017. Greater temporal variations are noted in the total annual rainfall series at Taipei station. For example, the maximum total annual rainfall is 4392 mm (1992) and the minimum is 1184 mm (2003). The mean and variance of total annual rainfall at Taipei station exhibit an increasing trend [41]. The altered rainfall regime at the Taipei station located in the Tamsui River basin is associated with complex anthropogenic interferences such as urban development and hydraulic-facility operation, leading to changes in streamflow regime.

Streamflow Data
Basic information on the seven selected streamflow gauge stations in the Tamsui River basin are summarized in Table 3. The daily streamflow of these stations, measured by the Water Resources Agency (Taiwan), are used to evaluate potential streamflow regime alterations within this basin. Due to the fact that parts of the streamflow data are unavailable (missing or unrecorded), the daily streamflow records are not continuous. At least 27-year daily streamflow records are used in this study. The longest streamflow data for 61 years are recorded at H4. It is worth noting that most stations have streamflow records for recent years (after 2014) except that H5 ends its records in 2000. Table 3. Basic information of seven selected streamflow gauge stations in Tamsui River basin, Taiwan.

Code
Station Name The streamflow regime is characterized by the annual maximum and minimum flow indices with various time scales, which include the annual maximum 1-day and 2-day flow (denoted as 1DMAX and 2DMAX) and the annual minimum 1-day, 7-day, and 30-day flow (denoted as 1DMIN, 7DMIN, and 30DMIN), since these annual flow indices evidently relate to flooding and water-supply problems in Taiwan. The annual mean and standard deviation of these flow indices for all selected stations are reported in Table 4. Among these stations, three stations (H2, H5, and H7) have upstream dams that potentially affect the streamflow regime, while the remaining four stations are free of dam impacts.

Reservoir Index
Basic information on 12 dams located within the Tamsui River basin, including dam height, capacity, catchment area, and finished year, are reported in Table 5. The damaffected stations (H2, H5, and H7) have different impacts due to the various numbers of upstream dams. Streamflow regime at H2 is influenced by D1 and D2. Streamflow regime at H5 is affected by D4, D5, D6, D7, D8, D9, and D10. Streamflow regime at H7 is influenced by D11 and D12. Figure 3 illustrates the time evolution of RIs, defined by equation (2), of these three stations within the streamflow records. It is worth noting that D1 was removed in 2007 due to dam failure. Evident high-impact RIs are observed at H2 and H5 since the maximum RI exceeds the threshold of 0.25, which is suggested by López and Francés [48]. Low-impact RI noted at H7 is attributed to smaller catchment areas.    Table 6 summaries the best-fitted probability distributions and the corresponding parameter types of the maximum and minimum flow indices for all stations in Tamsui River basin. Variations of mean and variance of the obtained distributions are also summarized in Table 6.  Stationary 1DMAX is observed at H2, H3, and H6 due to their constant distribution parameters. However, different distributional changes (i.e., different time-variation patterns of parameter) of 1DMAX are noted for the other stations. For example, H7 clearly has nonlinearly declined mean and variance, H5 has nonlinearly increasing variance, H1 and H4 have nonlinearly varying variance, and H4 and H5 have constant mean. Figure 4a-d illustrate curves of various quantiles (0.95, 0.75, 0.5, 0.25, and 0.05 from top to bottom) of 1DMAX at H1, H4, H5, and H7, respectively. Evidently different distributional changes are observed in these stations. For instance, H1 illustrates that the range of 1DMAX changes from wider to narrower and then becomes wider due to nonlinear variance. H7 clearly shows a declined range of 1DMAX caused by declined variance.

Time-Covariate Nonstationary Models of the Maximum and Minimum Flow Indices
Stationary 2DMAX is only observed at H3. The time-varying patterns of 2DMAX at the remaining stations are similar to those of 1DMAX with the exception of the nonstationary distribution parameters noted at H2 and H6. For example, a wider range of 2DMAX is observed at H1 in 1957, the narrowest 2DMAX is noted around 1975 and then gradually increases. The distributional change of 2DMAX at H1 is similar to that of 1DMAX shown in Figure 4a. A gradually reducing range of 2DMAX is noted at H7, which is also similar to the distributional change of 1DMAX shown in Figure 4d.
Non-stationarity is dominant in the minimum flow indices since none of these are categorized as stationarity. This phenomenon implies that the minimum flow is vulnerable to changing environments when compared with the results of the maximum flow indices.  Stationary 2DMAX is only observed at H3. The time-varying patterns of 2DMAX at the remaining stations are similar to those of 1DMAX with the exception of the nonstationary distribution parameters noted at H2 and H6. For example, a wider range of 2DMAX is observed at H1 in 1957, the narrowest 2DMAX is noted around 1975 and then gradually increases. The distributional change of 2DMAX at H1 is similar to that of 1DMAX shown in Figure 4a. A gradually reducing range of 2DMAX is noted at H7, which is also similar to the distributional change of 1DMAX shown in Figure 4d.
Non-stationarity is dominant in the minimum flow indices since none of these are categorized as stationarity. This phenomenon implies that the minimum flow is vulnerable to changing environments when compared with the results of the maximum flow indices. Figure 5a-d illustrate variations of quantile curves (0.95, 0.75, 0.5, 0.25, and 0.05 from top to bottom) of 1DMIN at H2, H3, H5, and H7, respectively. Different variation patterns are clearly noted in 1DMIN for various stations. For example, H1 (not shown in Figure 5), H2, and H7 exhibit an increasing trend of mean, while H2 and H3 have an increasing trend of variance. However, these trends, as well as the variation pattern of mean and variance of other stations, are nonlinear variations, except for the constant variance observed at H1. Similar distributional changes between 1DMIN and 7DMIN are observed at most stations. An exception is noted at H7 due to the nonlinear-varying variance of 1DMIN becoming a nonlinearly increasing variance of 7DMIN. Distributional changes of 30DMIN are close to those of 7DMIN at H1, H2, H4, and H6. Significantly decreasing variance of 30DMIN are observed at H3, H5, and H7. In summary, variation patterns among 1DMIN, 7DMIN, and 30DMIN are similar. The most distinct distributional changes among minimum flow indices are that the increasing means of 1DMIN and 7DMIN become a constant mean of 30DMIN, and increasing variance of 7DMIN becomes decreasing variance of 30DMIN at H7. Similar distributional changes between 1DMIN and 7DMIN are observed at most stations. An exception is noted at H7 due to the nonlinear-varying variance of 1DMIN becoming a nonlinearly increasing variance of 7DMIN. Distributional changes of 30DMIN are close to those of 7DMIN at H1, H2, H4, and H6. Significantly decreasing variance of 30DMIN are observed at H3, H5, and H7. In summary, variation patterns among 1DMIN, 7DMIN, and 30DMIN are similar. The most distinct distributional changes among minimum flow indices are that the increasing means of 1DMIN and 7DMIN become a constant mean of 30DMIN, and increasing variance of 7DMIN becomes decreasing variance of 30DMIN at H7.
To gain a consistent description of distributional changes of the maximum and minimum flow indices for various stations, the 16-category distributional-change scheme proposed by Shiau and Wu [41] is adopted in this study. This scheme categorizes distributional changes in terms of variation of mean and variance of probability. Four types of variation in mean and variance (linear or nonlinear increasing, no-change, linear or nonlinear decreasing, and complex variation not belonging to previous types) are adopted in this scheme and lead to 16 categories of distributional change, which are defined in Table 7. Categories of the maximum and minimum inflow indices of all stations are also reported in Table 6. The numbers of indices classified in each category are reported in Table 7.

RI-Covariate Nonstationary Models of Maximum and Minimum Flow Indices
Effects of dams on streamflow regime alterations are evaluated in terms of the reservoir index, which is used as the covariate in the models of the maximum and minimum flow indices for various time scales. Table 8 reports the best-fitted probability distributions of the maximum and minimum flow indices using the reservoir index as the covariate at H2, H5, and H7, since these stations are located downstream of dams. Table 8. Results of GAMLSS using RI as covariate.
The results indicate that stationary distributions (constant parameters with RI) are observed at H2 (1DMAX, 2DMAX, 7DMIN, and 30DMIN) and H7 (1DMIN, 7DMIN, and  30DMIN). Streamflow regime at H2 is influenced by D1 and D2, which are check dams used to counteract erosion. Since D1 and D2 are not used for regulating streamflow, the maximum and minimum flow indices at H2 are related less to RIs, although high RIs are noted at H2. Low RIs observed at H7 leads to the minimum flow indices at H7, independent of RIs. All the maximum and minimum flow indices at H5 depend on RI because there are seven upstream dams for the purpose of water supply and hydropower generation. Heavily regulating of the streamflow of upstream tributaries of H5 induces the highly altered streamflow regime at H5. That is, the distribution parameters of the maximum and minimum flow indices at H5 are related to RIs.

Discussion
The results of time-covariate modelling indicate that the streamflow regime of the Tamsui River basin is characterized by non-stationarity, since 31 out of 35 flow indices are best-fitted by the probabilistic models with time-dependent parameters. This phenomenon is prevalent in the minimum flow indices due to 100% of the flow indices belonging to nonstationary condition. This fact reveals that the minimum streamflow regime in the Tamsui River basin is highly vulnerable to climate change and anthropogenic activities. However, the effects of damming on streamflow characteristics of downstream stations (H2, H5, and H7) in terms of RI are insignificant since approximately half the indices (46.7%) do not relate to RIs.  Table 9 reports the AICs of the best-fitted distributions for the time-covariate and RI-covariate models. The lower-AIC model denotes the better model which is closer to the observed data. At H2, lower-AIC time-covariate models are observed in 2DMAX and all the minimum flow indices, except that an identical AIC is noted for the time-covariate and RI-covariate 1DMAXs, because both are fitted as stationary models. Lower-AICs RI-covariate models are observed in 1DMAX and 2DMAX and lower-AIC time-covariate models are observed in the minimum flow indices at H5. All the lower-AICs are noted for the time-covariate models at H7. The results indicate that the time-covariate models generally outperform the RI-covariate models.
Better time-covariate models are attributed to several factors. One of the reasons is that dams are not the most dominant factor in altering streamflow regime. Other anthropogenic factors including catchment development, population growth, and urban drainage, alter the streamflow regime. Besides, the low-RI such as 0.009 noted at H7 is also a factor leading to an insignificant covariate of RI. The purposes and operation of dams may cause different effects on streamflow changes. For instance, check dams such as D1 and D2 located upstream of H2 are used to reduce erosion, not to store and regulate streamflow, and have less effects on streamflow alterations. These factors, not involved in RI, make RI less related to streamflow alterations in Tamsui River basin.
Yeh et al. [39] indicated that the annual high flow (Q 0.9 , daily discharge exceeded by 10% of days in a year) at H1 exhibited an insignificant increasing trend and an approximate null trend at H4 based on the Mann-Kendall test. Yeh et al. [39] also indicated that annual low flow (Q 0.1 , daily discharge exceeded by 90% of days in a year) had insignificant decline and increasing trends at H1 and H4, respectively. Due to different definition and data length of flow indices, the results obtained by Yeh et al. [39] are not exactly identical to the results of this study, but parts of the results are consistent. For example, constant mean of 1DMAX at H4 and increasing mean of 1DMIN at H1 noted in this study are similar to the Mann-Kendall test results obtained by Yeh et al. [39]. Shiau and Wu [41] revealed that clearly increasing variance of the maximum 1-day rainfall is noted after 2000 at Taipei station in the Tamsui River basin. This result is similar to the increasing variance observed at 1DMAX at H1, H4, and H5, but contradicts the decreasing variance of 1DMAX at H7. This fact implies that alterations in streamflow regime are not affected only by rainfall, but also by other complex factors such as watershed development, hydraulic-facility operation, the relationship between rainfall and streamflow, and others. The obtained nonstationary probability distributions of the flow indices offer more comprehensive information than the monotonic increasing or decreasing trends. For example, 1DMIN at H2 exhibits increasing mean and increasing variance (Category I and shown in Figure 5a), which is insufficiently described by increasing trend. However, the obtained best-fitted distributions are limited to the data periods used in the study; that is, the obtained time-varying models are not properly able to predict far-future streamflow conditions, indicated by Villarini et al. [51].

Conclusions
The main aim of this study is to explore non-stationarity of the maximum and minimum flow indices in terms of distributional changes for seven streamflow gauge stations in the Tamsui River located in northern Taiwan using the GAMLSS (generalized additive models for location, scale, and shape).
Results of time-covariate models indicate that stationarity is only noted in the maximum flow index, i.e., 1DMAX at H2, H3, and H6, and 2DMAX at H3. The streamflow regime in the tributaries (H3 and H6) is less influenced than those in the main rivers. Non-stationarity is prevalent in the minimum flow indices since all minimum flow indices are classified as nonstationary distributions. Clearly different distributional changes are observed among the flow indices. Streamflow regime in the Tamsui River generally exhibits non-stationarity because that 31 out of 35 (88.6%) flow indices have time-dependent distribution parameters. Category XVI (complex variations of mean and variance) is the dominant category since 34.3% of flow indices are classified in this category.
Results of RI-covariate models at H2, H5, and H7 indicate that constant parameters are observed at H2 (1DMAX, 2DMAX, 7DMIN, and 30DMIN) and H7 (1DMIN, 7DMIN,  and 30DMIN). The reasons that these flow indices do not relate to RI are attributed to the dams located upstream of H2 used for counteracting erosion and the low RI (0.009) noted at H7. The RI-covariate models are generally worse than the time-covariate models. Only the RI-covariate 1DMAX and 2DMAX at H5 outperform the time-covariate models. This fact indicates that construction of dams is not the only factor altering streamflow regime in the Tamsui River.
Based on the streamflow records used in this study, the results reveal that the streamflow regime in the Tamsui River is significantly altered except for 1DMAX and 2DMAX at some stations. Alterations in streamflow regime are attributed to climate change as well as to anthropogenic interferences. Prevalent non-stationarity in streamflow regime (88.6 % of flow indices) leads to approximate one third of flow indices (34.3%) exhibiting complex variation patterns instead of monotonic increasing or decreasing trends. RI is insufficient to describe dam-altered streamflow regime in the Tamsui River due to time-covariate modeling of 80% flow indices outperforming RI-covariate modeling. Complex anthropogenic interferences other than dam-construction need to be considered in the Tamsui River basin to evaluate spatially diverse alterations in streamflow regime.
Identification of non-stationarity in streamflow series gives an insight into streamflow regime alterations. This information can usefully guide future design and current operation of water-resources engineering, adapting to the changing environment. This study analyzes the non-stationarity in the maximum and minimum flow indices in the Tamsui River basin located in northern Taiwan. Construction of dams is not the unique factor influencing streamflow regime. Modification of the reservoir index to include purpose or operation types and other anthropogenic effects such as population or urban drainage remain as topics for further extending of this study.