Trend Detection in Annual Streamflow Extremes in Brazil

Abstract

Changes in streamflow extremes can affect the economy and are likely to impact the most vulnerable in society. Estimating these changes is crucial to develop rational adaptation strategies and to protect society. Streamflow data from 1106 gauges were used to provide a comprehensive analysis of change in eight different extreme indices. The modified trend-free prewhitening and the false discovery rate were used to account for serial correlation and multiplicity in regional analysis, issues shown here to distort the results if not properly addressed. The estimated proportion of gauges with significant trends in low and high flows was about 23% and 15%, respectively. Half of these significant gauges had more than 60 years of data and were associated with changes greater than 5% per decade. A clear spatial pattern was identified, where most increasing trends in both low and high flows were observed in Southern Brazil, and decreasing trends in the remaining regions, except for the Amazon, where a pattern is not clear, and the proportion of significant gauges is low. Results based only on gauges unaffected by reservoirs suggest that reservoirs alone do not explain the increasing trends of low flows in the southern regions nor the decreasing trends in high flows in the remaining hydrographic regions.

1. Introduction

Understanding how climate and human activities impact flow and storage of water on the planet is the subject of current research efforts by the international scientific community. It was the main subject of “Panta Rhei—Everything Flows” (2013–2022), an initiative of the International Association of Hydrological Sciences (IAHS) dedicated to the relationship between hydrologic behavior and changes in society. Similarly in the United States [1], research priorities of the U.S. Geological Survey for water shifted to focus on how climate and anthropogenic activities affect and impact quantity, quality, and reliability of water resources systems.

Water resources systems are normally planned, designed, and operated based on the premise that streamflow series are stationary [2]. Although very convenient, the assumption of stationarity has been the topic of intense discussion in the literature. Milly et al. [2] challenged the stationarity assumption and initiated a long debate over whether the stationary assumption is adequate to model streamflow, or if a nonstationary framework should be implemented [3,4,5].

There are many natural and anthropogenic factors that can change the behavior of streamflow over time. Land use and land cover changes due to expansion of urban or agricultural areas, construction of water infrastructure for water supply, energy production and flood control, and increase in surface imperviousness and groundwater extraction for agricultural and industrial uses are all examples of human activities in the basin that alter streamflow behavior. Additionally, increases in greenhouse gases and aerosols in the atmosphere, largely due to human activities, are altering the frequency and intensity of heavy storms and the pattern of monsoon precipitation, which may impact the frequency and magnitude of the occurrence of floods and droughts [6]. Furthermore, natural and large-scale climate phenomena, such as El Niño–southern oscillation (ENSO), Pacific decadal oscillation (PDO), and Atlantic multidecadal oscillation (AMO), affect the precipitation pattern in many regions of the world, which can change the behavior of flow extrema in time scale of a few years [7,8,9].

The nonstationary behavior of flow extremes is probably the most critical aspect of streamflow as it can have a large impact on society [5,10,11,12,13]. Changes in low-flow characteristics may weaken the reliability of water supply systems, reduce crop production, aggravate water quality problems, constrain the transportation of goods through a river system, and diminish hydroelectric production. Conversely, a change in frequency and magnitude of flood events may cause an increase in property damages, transportation disruption, and human and livestock losses [5,11,12].

Understanding the behavior of changes in streamflow characteristics is important to estimate mean annual risks for the future, crucial for decision makers to develop adaptation and mitigation strategies [5,12,13,14]. Hence, detecting changes in streamflow series is an important step in water resources planning and management [11,12,15,16].

Changes in streamflow is often assessed by applying null hypothesis significance tests (NHSTs) to the observed series of the flow index of interest [4,17]. NHSTs are based on the null distribution (H₀) of the test statistic, which is derived based on the assumption that the series is stationary. One uses the null distribution and the sampling estimate of the test statistic to decide whether the series is stationary (accept H₀) or nonstationary (reject H₀). There is no definite answer to this question because there is always a probability of mistakenly rejecting H₀ when the series is in fact stationary (Type 1 Error), which may result in overpreparedness or erroneously accepting H₀ when the trend is present (type 2 error), which may induce underpreparedness [11]. Traditionally, NHST procedures employ a significant level (e.g., α = 0.05) to control type 1 error, neglecting type 2 errors.

For the last two decades, trend detection studies of streamflow series have been carried out in many parts of the world at different spatial scales, including local and regional scales within a single country, such as in the United States [18,19,20,21,22], Australia [23,24,25], China [26,27,28,29], France [30,31], Great Britain and the United Kingdom [32,33], Germany [34], Spain [35,36], Turkey [37], Ireland [38], Poland [39], and Canada [40,41,42,43]. There are also studies carried out at the continental scale, such as those executed for Europe [44,45,46,47,48], and even at the global scale [7,15,49,50,51].

However, there are many challenges that need to be addressed in typical trend detection methods to strength the approach. For instance, the presence of autocorrelation, common for annual low flows, increases the chance of finding a significant trend when a trend is not present, making it difficult to associate with the nominal significance level of the test [7,52,53]. The effects of serial correlation in trend detection analysis have long been mentioned in the statistical literature [54,55,56], even though a definitive solution is yet to be found. To the best of the authors’ knowledge, Lettenmayer [57] was the first to point out in the hydrologic literature that the presence of positive serial correlation results in a probability of type 1 error that is larger than the level of significance of the test. These results were later shown by von Storch [58] and discussed in several studies [41,59,60,61,62,63]. Two main strategies have been developed to deal with this issue: (a) change the original time series by removing the serial correlation and (b) adjust the null distribution of the test statistic to account for the serial correlation. The former strategy is adopted in popular methods, such as the prewhitening (PW), introduced by von Storch [58], and the trend-free prewhitening (TFPW), developed by Yue et al. [64] to resolve the low statistical power provided by PW. The second strategy can be implemented in different ways, for instance, using a correction of the variance of the test statistics [65], by deriving the entire null distribution of the test statistics [62] or by using sampling techniques, such as block bootstrap, suggested by Onoz and Bayazit [63]. A modified version of the TFPW, introduced by Onoz and Bayazit [63], is used in this study as it has been shown to provide a better balance between the type 1 and type 2 errors than the TFPW, the former being directly related to the statistical power. A more in-depth discussion of this method is provided later in the paper.

Another aspect of NHSTs that has been neglected for a long time in trend detection studies is related to the simultaneous evaluation of multiple null hypothesis tests, known as multiplicity of tests [66,67]. Like autocorrelation, multiplicity in statistical tests tends to overestimate the number of gauges with a significant trend in a region. One common strategy to deal with multiplicity in statistics in general is to use either the Walker or the Bonferroni test [66]. Although both tests correctly control the type 1 error in a region, they have small statistical power, which means that the probability of identifying a trend when a trend is present in the series is low. Over the years, many different methods have been proposed in the statistical literature to deal with multiplicity [68,69,70]. The concept of the false discovery rate (FDR), introduced by Benjamini and Hochberg, [70] based on the expected value of the proportion of false positives among the rejected null hypotheses became very popular in many fields of science, especially in genetics and epidemiology. The FDR controls the type 1 error in a framework where multiple null hypotheses are evaluated simultaneously but has greater statistical power than the Walker and Bonferroni tests. The multiplicity issue has only recently been recognized as an important problem in trend analysis, with an increasing number of papers [29,31,41,71,72,73,74,75,76] adopting methods developed in other fields of research, such as the false discovery rate developed by Benjamini and Hochberg [70].

In Brazil, studies that investigate changes in streamflow extremes based on observed flow data at the national scale are limited. Unlike studies that have focused on precipitation extremes, studies on streamflow have relied on relatively small datasets [77,78,79,80,81,82,83,84,85]. Apart from [83], we could not find annual maxima flood trend studies at the national scale. The issue of autocorrelation is rarely considered, and there are practically no studies that demonstrate the impact of temporal dependence on trend detection results carried out for Brazil. Although all trend studies have evaluated the proportion of significant results at a given region, the issue of multiplicity of tests is ignored.

This paper aims to fill the gaps discussed in the previous paragraph by providing a comprehensive trend detection analysis for eight streamflow extreme indices in the 12 hydrographic regions of Brazil using daily data from 1106 gauges. The effects of serial correlation in the data and the multiplicity of test in regional analysis are discussed and properly addressed in the methodology, avoiding the excessive number of gauges with false positive trends.

2. Materials and Methods

2.1. Study Area

The objective of this study is to provide a comprehensive picture of changes in streamflow extremes for the whole country so water managers can easily identify regions or subregions where action is needed. Brazil has the highest total renewable freshwater supply on the planet [86]; however, water resources are unequally distributed in space and time. Thus, water availability must be viewed as an indicator of the average situation. In fact, 80% of the water in the country is in the Amazon River Basin, where only 5% of the population lives. Approximately 93% of the population depends on just 30% of the total water available [86]. These uneven distributions of supply and demand are the source of many water-related problems in the country. Changes in the amount of water available, especially in the pattern of low flows, may impact water quality, reliability of water supply systems, and hydropower generation. Additionally, the occurrence of floods is also a major concern in the country, where approximately 7.7 million people have been impacted by flooding from 2013 to 2016, highlighting how changes in the patterns of high flows are also an import issue.

Because of the continental size of the country, Brazil is divided into 12 hydrographic regions for water resources planning and management, as illustrated in Figure 1: Amazon (AMZ), Western Northeast Atlantic (ANC), Parnaíba (PNB), Eastern Northeast Atlantic (ANO), Tocantins–Araguaia (TOA), São Francisco (SFO), East Atlantic (ALE), Paraguai (PRG), Paraná (PRN), Southeast Atlantic (ASD), South Atlantic (ATS), and Uruguai (URU). These hydrographic regions are used in this study as independent spatial units, acting as the basis for the analysis of trend detection in annual streamflow extremes.

2.2. Streamflow Data

This study employed daily streamflow data from 1106 selected gauges located across Brazil. Only gauges with at least 30 years of record length and with at least 5 years of data in the period 2000–2015 were selected for the study. The latter criterion was used to avoid working with gauges that had been inactive for many years.

Figure 2 presents the location of these gauges, where colors indicate record length. Gauges with a long record length are located primarily in the eastern part of Brazil. The longest record length was 92 years. A quarter of the 1106 gauges used in the study has at least 60 years of data, and half have at least 43 years. Catchment sizes range from 12 to 4,670,000 km² with a median of 3220 km². The majority (80%) of the gauges have drainage area less than 20,000 km², and only 10% have drainage areas larger than 50,000 km². Figure 3 presents the proportion and number of gauges for each hydrographic region stratified by classes of record length. The PRN hydrographic region has almost 30% (325) of all gauges used in the study, but the ASD has the most gauges (92) with a record length greater than 60. The PNB hydrographic region has the smallest number of gauges with only 16.

2.3. Streamflow Indices

Table 1. Streamflow indices.

Index	Definition
QX1d	Annual maximum daily flow
QX5d	Annual maximum 5-day consecutive average flow
QX30d	Annual maximum 30-day consecutive average flow
Qmean	Mean annual streamflow
Q7	Annual minimum 7-day consecutive average flow
Q30	Annual minimum 30-day consecutive average flow
Q7Wtri	Annual minimum 7-day consecutive average flow (wettest trimester)
Q7Wsem	Annual minimum 7-day consecutive average flow (wettest semester)

presents eight indices that were selected to describe different aspects of measured daily streamflow. The use of such indices is important to provide a clear picture of possible changes in flow patterns and how these changes may impact different parts of society. For instance, changes in low flows can have an impact on the reliability of water supply and hydroelectric generation, whereas changes in high flows may dramatically impact flood risk.

2.4. Presence of Reservoirs

Although attribution is not the objective of this study, we decided to perform the analysis for two datasets, one with the whole set of selected gauges and another that contains only gauges that are not affected by reservoirs.

To accomplish this task, we used a database that contains information about natural and artificial water bodies in Brazil (MDA). There are 240,899 water bodies currently registered in the MDA, summing up to a total surface area of 173,749.56 km². About a quarter of them (66,372) are composed of natural water bodies with a total area of 128,165.89 km². The remaining are artificial water bodies that account for a total area of 45,583.76 km².

The criterion adopted to decide whether reservoirs located upstream potentially affect the flows of a gauge was the degree of regulation (DoR). DoR is defined as the ratio between the sum of all reservoirs’ storage capacity located upstream and the long-term mean annual flow [87,88,89] given by:

$$Do{R}_j=\frac{1}{Qml{t}_j}{\displaystyle \sum }_{i=1}^{n_j}{V}_{ij}$$

(1)

where $V_{ij}$ is the storage capacity of reservoir i located upstream at gauge j, $n_j$ is the number of reservoirs upstream, and $Qml{t}_j$ is the mean annual streamflow, in volume, at gauge j.

There is no predefined DoR threshold value above which one can say flows are affected by reservoirs upstream. Here, we follow [89,90] and declare that gauges with $DoR\le 0.02$ are unaffected by the presence of reservoirs. Out of the 1106 gauges selected for the study, 670 were considered unaffected by reservoirs, about 60% of the total number of gauges. The choice of the threshold of 0.02 for DoR is subjective, and the most appropriate value for this analysis may depend on the flow index of interest. This threshold value should be better evaluated in future assessments. Figure 4 provides a general picture of the values of DoR at the national scale. In this case, DoR was computed for each river stretch of the national river network.

2.5. Trend Detection Framework

This section provides the methodological steps used in this study to identify trends and estimate their magnitudes. Figure 5 gives a general overview of the strategy used. This section continues with a description of the statistical test used to assess trend significance and the estimator of the magnitude of a trend, followed by a discussion on the methods used to deal with both serial correlation and multiplicity of tests.

2.5.1. Mann–Kendall Trend Test

The Mann–Kendall (MK) trend test [91,92] was adopted in this study to identify statistically significant monotonic increasing or decreasing trends. The MK test is a nonparametric approach that does not require any assumption about the probabilistic model. The test is also applicable to data containing outliers or nonlinear trends [46]. The MK test is widely used in trend detection studies because it is distribution-free and robust against outliers and has a higher power for non-normally distributed data [63,93]. In addition, it has been used in most previous streamflow trend analyses [7,48,52,94]. Considering {x₁, x₂, x₃, . . ., x_n} a time series of length n, then the MK test statistic S is given by:

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sign\left(x_j-x_i\right)$$

(2)

where

$$sign\left(x_j-x_i\right)=\left\{\begin{array}{c}1 ,se x_j>x_i\\ 0, se x_j=x_i\\ -1, se x_j<x_i\end{array}\right.$$

(3)

The null hypothesis H₀ for the test is ‘‘there is no trend in the time series”. If H₀ is true, then S is normally distributed with zero mean and variance ${\sigma }_0^2$,

$${\sigma }_0^2=\frac{n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{j=1}^m{t}_{j}j\left(j-1\right)(2j+5) }{18}$$

(4)

where n is the number of data points, m is the number of tied groups, and t_i is the number of ties of size i. A tied group is a set of sample data having the same values. In cases where the sample size n is greater than 10, S is approximately normally distributed, which justifies the use of the standard normal variate, Z, as the test statistic:

$$Z=\left\{\begin{array}{c}\frac{S-1}{{\sigma }_0},\hspace{1em}\mathrm{if} S>0\\ 0\hspace{1em},\hspace{1em}\mathrm{if} S=0\\ \frac{S+1}{{\sigma }_0},\hspace{1em}\mathrm{if} S<0\end{array}\right.$$

(5)

where the unity in the numerator is used for continuity issues in small samples since S is discrete. A positive (negative) value of Z indicates an increasing (decreasing) trend. If $z<z_{\alpha /2}$ or $z>z_{1-\alpha /2}$, the null hypothesis that trend is absent can be rejected at α% significance level. In this study, the α significance level was 5%.

2.5.2. Theil–Sen Slope Estimator

Here, we use the Theil–Sen approach (TS) [95] to estimate the actual magnitude of the trend, $\beta$, which is important to understand the practical implications of detected trends and to evaluate the costs of the associated impacts. The nonparametric TS estimator of the magnitude of trend is based on the slope of all $n\left(n-1\right)$ distinct pairs of observations in a sample of size n, denoted here as $d_{ij}$:

$$d_{ij}=\frac{\left(x_j-x_i\right)}{\left(j-i\right)}$$

(6)

where x_i and x_j are the data values at times $1\le i<j\le n$. The decadal relative TS slope estimator of β is given by:

$$\widehat{\beta }=\frac{\left(med\left\{d_{ij}\right\}\right)}{\overline{x}}\times 10$$

(7)

where $med\left\{d_i\right\}$ is the median of all values of $d_{ij}$ and $\overline{x}$ is the sample mean.

2.5.3. Adjustment for Autocorrelated Data

The MK test assumes that the input data are serially independent. The presence of a positive serial correlation in the data structure makes the test more liberal, meaning it rejects the null hypothesis more often than its significance level suggests [58,60,96,97].

The presence of serial correlation is assessed by a statistical test based on the following approximation of the 95% confidence interval for the population lag-1 serial correlation ${\varrho }_1$ [98,99,100]:

$$\frac{-1-1.96\sqrt{n-2}}{n-1}\le {\varrho }_1\le \frac{-1+1.96\sqrt{n-2}}{n-1}$$

(8)

The true and unknown lag-1 serial correlation ${\varrho }_1$ is usually estimated using the estimator r₁, given by:

$$r_1=\frac{{{\displaystyle \sum }}_{t=1}^{n-1}\left(x_t-\overline{x}\right)\left(x_{t-1}-\overline{x}\right)}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2}$$

(9)

where $\overline{x}$ is the sample mean. If the sample lag-1 serial correlation, r₁, falls within the interval given by Equation (8), the sample data are assumed to be serially independent. Otherwise, the sample data are considered serially correlated.

As the estimator of the autocorrelation r₁ has a downward bias (negatively) [92], we also use a bias-corrected estimator based on [101]:

$$r_1^{\ast }=\left(\frac{n{r}_1+2}{n-4}\right)$$

(10)

In order to eliminate the influence of serial correlation on the MK test, von Storch [58] proposed the removal of the serial correlation component from the time series based on the assumption that the series comes from a lag-1 autoregressive process, AR(1). The MK test is then applied to the modified series to assess the significance of a possible trend. The PW strategy is very popular in trend detection studies, which accounts for serial correlation.

Yue and Wang [97], as well as Yue at al. [102], showed that the PW procedure can seriously distort the possibility of the MK test to detect a trend. When a trend is present in the series, the estimation of the correlation component of the AR (1) model is inflated. When PW removes the inflated correlation component of the series, it ends up removing part of the trend, causing the distortion in the results of the MK test, reducing the power of the test.

In order to mitigate these limitations, Yue et al. [64] proposed a modified PW procedure called TFPW. In this procedure, the magnitude (slope) of the trend is first estimated, and the original series is detrended. Then the lag-1 serial correlation coefficient of the detrended series is estimated, and the original series is prewhitened using this estimate. Finally, the identified trend is added to the prewhitened series. The Mann–Kendall test is then applied to the resulting series to assess the significance of the trend. The authors argued that the removal of the trend as a first step may allow for a more accurate estimate of the population’s lag-1 autocorrelation coefficient and a subsequent better estimation of the significance of the trend.

The TFPW procedure was considered the solution to the low level of power attributed to the PW strategy. However, later studies [62,65] demonstrated that the sampling uncertainty in the estimation of the trend led to the violation of the type 1 error. In some cases, the TFPW’s inability to control the type 1 error was similar to the MK’s in the original problem.

Önöz et al. [63] proposed a modified TFPW (MTFPW) procedure to better control type 1 error. Both procedures use the estimated serial correlation from the detrended series, $X_t^d$, thus avoiding overestimation due to the presence of a trend. However, unlike TFPW, which applies the PW procedure to $X_t^d$ to obtain a detrended and prewhitened series, $X_t^{dPW}$, the MTFPW applies the PW procedure to the original series to obtain a trended prewhitened series to which the MK test is applied to assess the significance of the trend. The outline of the MTFPW procedure is as follows:

• Estimate the magnitude of the trend, $\beta$, using Equation (7).
• Obtain a detrended series, $X_t^d$, by removing the estimated trend from of the original series, $X_t^d=X_t-\widehat{\beta }{\overline{X}}_{t}t/10$, where t is the time interval.
• Estimate an unbiased sample autocorrelation (r₁) of the detrended series, $X_t^d$.
• If r₁ is not statistically different from zero, then the MK test is applied to the original series, $X_t$. Otherwise, the PW procedure is applied to the original series, $X_t$, to obtain a prewhitened series, $Y_t^{MTFPW}=X_t-r_1{X}_{t-1}$.
• Apply the MK test to $Y_t^{MTFPW}$ to check the significance of the trend.

The MTFPW provides a good tradeoff between type 1 and type 2 errors. The TFPW has a much larger statistical power than the PW method, but it fails to control the type 1 error. The MTFPW, based on a relatively small change in the TFPW method, provides a better control of the type 1 error without losing much statistical power.

2.5.4. The Multiplicity Problem

Trend detection studies usually consist of applying simultaneous hypothesis tests at many gauges in a region. A significant level ${\alpha }_L$ for the local test (one specific gauge) is specified a priori based on the risk of making a type 1 error. Therefore, if the null hypothesis is true, the probability of making a type 1 error for the specific gauge is simply ${\alpha }_L$.

When the goal is to detect the proportion of gauges with trends in a region with many gauges, it is important to control the type 1 error for the whole region, not just for each individual gauge. Ignoring the control of the type 1 error for the region results in a violation of the nominal significance level ${\alpha }_R$ at the regional level. This problem is called multiplicity of hypothesis tests in statistics.

Here, we use the concept of false discovery rate (FDR), first introduced by Benjamini and Hochberg [70] and defined as the expected value of the proportion of false detections among all rejections:

$$\mathrm{FDR}=\left\{\begin{array}{c}E\left[\frac{N_{10}}{R}\right], if R>0\\ 0 , if R=0\end{array}\right.$$

(11)

where $N_{10}$ is the number of gauges where the null hypothesis was mistakenly rejected. The FDR is widely used in the genetics and epidemiology fields, but only recently, it has gained popularity in hydrologic problems [31,66,67,72,103,104].

The Benjamini–Hochberg (BH) procedure that uses the FDR concept is based on the ordered p-values $P_{\left(1\right)}\le P_{\left(2\right)}\le \dots \le P_{\left(m\right)}$, associated with the hypotheses $H_{\left(1\right)}\le H_{\left(2\right)}\le \dots \le H_{\left(m\right)}$ and the critical values $d_1\le d_2\le \dots \le d_m$, considered to be equal to $d_i=i{\alpha }_R/m$. The BH procedure is sequential, meaning the decision to reject or accept the null hypothesis is made gauge by gauge starting from the one with the largest p-value (less evidence to reject the null). If $P_{\left(m\right)}<{\alpha }_R$, then all gauges are declared nonstationary; otherwise, the condition $P_{\left(j\right)}<d_j$ is tested until it is true. For the first time $P_{\left(j\right)}<d_j$, all hypotheses $H_{\left(i\right)<j}$ are rejected.

3. Results and Discussion

3.1. Impacts of Serial Correlation and Multiplicity of Tests

This section illustrates the impacts the presence of serial correlation and the multiplicity of hypothesis tests have on trend detection from the MK test. To assess these impacts, we used three different approaches to detect trends in the eight streamflow indices (Table 1): MK, MK-MTFPW, and MK-MTFPW-FDR. The MK approach consists of applying the MK test to the series with no concern for whether observations are correlated or independent or to control the type 1 error at the regional level. The approach denoted as MK-MTFPW, as its name suggests, deals with the presence of autocorrelation using the MTFPW method but ignores the multiplicity issue. Finally, the MK-MTFPW-FDR approach deals with both issues using the FDR to control the type 1 error at the regional level. Both local and regional significance levels used in this study were equal to 5%.

Figure 6 shows the proportion of gauges with a significant trend for each one of the eight streamflow indices using the three different approaches. One can see, as expected, that ignoring the presence of a serial correlation and the control of the type 1 error at the regional level results in a larger proportion of gauges with significant trends, possibly with a large number of false positives, when compared with the other two approaches.

When we compare the results obtained by the MK (red bar) with those obtained by the MK-MTFPW (black bar), the differences are larger for low-flow indices compared with high-flow indices. This makes sense as a serial correlation is much stronger among indices of low flow. These results show the importance of using a procedure to consider the serial correlation when it is present in the data. For instance, the proportion of gauges with significant trends in Q7 series is about 55% if the MK procedure is used but only approximately 35% when the serial correlation is appropriately considered in the analysis. Ignoring the serial correlation can heavily distort the results and may provide a completely different picture of reality.

Ignoring the multiplicity problem when simultaneously applying hypothesis tests for many gauges in a region can also result in large numbers of false positives at the regional level. One can see important differences in the proportion of gauges with significant trends when comparing the MK-MTFPW (black bar) and the MK-MTFPW-FDR (blue bar). However, unlike the differences between the MK and the MK-MTFPW, due to serial the correlation alone, here, these differences are equally proportional regardless of the streamflow indices.

Although results presented in Figure 6 provide an enlightening understanding of the effects of serial correlation and multiplicity in the trend detection, they fail to provide a complete picture of the situation. Neither is the signal of the significant trends explicitly given nor is the spatial distribution of the significant gauges presented on a map.

Both the signal of the significant trends and the distribution of gauges on a map are presented in Figure 7. Because of the lack of space, results are presented only for Q7. Red (blue) circles represent gauges with decreasing (increasing) trends, whereas black dots represent gauges with nonsignificant results. The three panels in Figure 5 show the results obtained by each of the three approaches used in the analysis.

Regardless of the approach used, one can notice a clear spatial pattern with the majority of decreasing significant trends happening in the northern and eastern parts of the country, except for the Amazon region, where no clear pattern is identified, and a greater part of gauges with increasing significant trends in the southern part of the country. As we move from the MK results to the MK-MTFPW and MK-MTFPW-FDR results, the density of red or blue circles diminishes, but the spatial consistency remains. In fact, it was observed that the use of the FDR procedure, which controls the expected number of false positives, avoided some regional inconsistencies obtained with the MK and MK-MTFPW, where some gauges located close to each other presented significant trends with different signs. The use of the MK-MTFPW-FDR approach resulted in more clearly identified clusters of gauges with significant trends.

3.2. Trend Analysis in Hydrographic Regions

The previous section showed the impacts that both serial correlation and multiplicity have on the results of trend detection. From this point onwards, all the results presented here were obtained by the MK-MTFPW-FDR approach, which accounts for both issues. Moreover, the control of the type 1 error at the regional level using the FDR procedure was applied separately for each hydrographic region.

Table 2. Percentage of MK trend test results for the 1106 stations across Brazil (NS = nonsignificant and S = significant) for the eight streamflow indices stratified in three classes of record length and three classes of relative change per decade ($\widehat{\beta })$ .

Indice	Res	$\widehat{\mathit{\beta }}$			$\widehat{\mathit{\beta }}$			$-5\% > \widehat{\mathit{\beta } }>-10\%$			$5\% <\widehat{ \mathit{\beta } }< 10\%$			$0\% > \widehat{\mathit{\beta } }>-5\%$			$0\% < \widehat{\mathit{\beta }}< 5\%$			Tot
		30–44 years	45–59 years	>60 years	30–44 years	45–59 years	>60 years	30–44 years	45–59 years	>60 years	30–44 years	45–59 years	>60 years	30–44 years	45–59 years	>60 years	30–44 years	45–59 years	>60 years
QX1d	NS	6%	1%	0%	2%	0%	0%	86%	2%	1%	5%	3%	2%	12%	6%	9%	11%	7%	10%	86%
		7%			3%			11%			9%			27%			28%
	S	4%	1%	1%	1%	0%	1%	0%	0%	1%	1%	1%	2%	0%	0%	0%	0%	0%	0%	14%
		6%			2%			2%			3%			0%			1%
QX5d	NS	7%	1%	0%	2%	0%	0%	9%	3%	1%	5%	3%	2%	10%	6%	9%	11%	6%	10%	87%
		8%			2%			13%			10%			26%			27%
	S	4%	2%	1%	0%	1%	0%	0%	0%	1%	1%	0%	2%	0%	0%	0%	0%	0%	1%	13%
		6%			1%			2%			3%			1%			1%
QX30d	NS	6%	1%	0%	2%	0%	0%	9%	3%	2%	5%	3%	2%	10%	6%	11%	11%	6%	8%	84%
		8%			2%			13%			9%			27%			25%
	S	5%	2%	1%	0%	0%	0%	1%	0%	2%	1%	0%	3%	0%	0%	0%	0%	0%	0%	16%
		7%			1%			3%			4%			1%			0%
Qmean	NS	8%	2%	0%	2%	0%	0%	7%	3%	2%	6%	2%	1%	10%	8%	10%	9%	5%	7%	81%
		11%			2%			11%			8%			28%			21%
	S	6%	2%	0%	0%	0%	0%	1%	1%	2%	1%	1%	4%	0%	0%	0%	0%	0%	0%	19%
		8%			1%			3%			6%			0%			0%
Q7	NS	10%	3%	1%	4%	0%	0%	7%	3%	2%	3%	3%	1%	8%	6%	8%	7%	3%	7%	77%
		14%			5%			12%			7%			22%			16%
	S	7%	2%	1%	2%	1%	2%	1%	1%	1%	0%	0%	3%	0%	0%	0%	0%	0%	0%	23%
		11%			4%			3%			3%			1%			0%
Q30	NS	8%	3%	0%	4%	1%	0%	7%	4%	2%	4%	2%	1%	9%	6%	9%	6%	3%	6%	77%
		12%			5%			13%			7%			24%			16%
	S	8%	2%	1%	2%	1%	1%	2%	1%	1%	0%	1%	3%	0%	0%	0%	0%	0%	0%	23%
		11%			4%			4%			4%			0%			0%
Q7WTri	NS	8%	3%	0%	2%	0%	0%	5%	4%	2%	3%	2%	1%	8%	7%	9%	9%	5%	8%	78%
		11%			2%			12%			7%			24%			21%
	S	7%	5%	2%	1%	0%	0%	0%	1%	2%	0%	0%	2%	0%	0%	0%	0%	0%	0%	22%
		14%			1%			3%			2%			1%			0%
Q7WSem	NS	7%	2%	0%	3%	0%	0%	8%	4%	2%	5%	2%	2%	9%	5%	8%	7%	4%	7%	76%
		9%			3%			13%			9%			23%			19%
	S	7%	4%	1%	1%	1%	1%	1%	2%	2%	0%	0%	2%	0%	0%	1%	0%	0%	0%	24%
		13%			2%			5%			3%			1%			0%

Figure 8 displays the results across the 12 hydrographic regions using stacked bars to show both the proportion and the actual number of gauges with significant trends for each hydrographic regions, and colors to differentiate decreasing (red) from increasing (blue) trends. Only the results for Q7, Q_mean, QX1d, and QX30d are presented, but they provide a representation of the situation.

presents the results in terms of the proportion of gauges with significant trends at the national level for all streamflow indices. The results are stratified based on three classes of record length: short (30–44 years), medium (45–59 years), and long (>60 years), and three classes of relative change per decade: small ($\left|\widehat{\beta }\right|<5\%$), medium ($5<\left|\widehat{\beta }\right|<5\%$), and large ($\left|\widehat{\beta }\right|>10\%$). Partitioning the results into classes of record length and relative change per decade helps the interpretation of the results.

Gauges with longer record length are associated with more precise results in terms of detection and estimation of trend magnitude since they usually result in tests with more statistical power. Therefore, gauges with large relative changes per decade and long, or at least medium, record lengths stand out from the rest as results are more precise and potential impacts tend to be more relevant due to the magnitude of changes. Water managers and other professionals in charge should study gauges in these classes more carefully and perhaps initiate an assessment of possible impacts in society. Gauges with short record length should be treated carefully as their trend detection results are more uncertain due to low statistical power and the likely overestimation of trend magnitude. Further, gauges with small relative change per decade, at least in a preliminary analysis, are of less concern. For instance, even a small relative change in Q7 in an already water-stressed region may be enough to threaten the water supply system of a community. These discussions, which are out of the scope of this paper, should be kept in mind when interpreting the results. The reader interested in issues related to differences between statistical and practical significance is referred to [105,106,107] for a general discussion and to [64,93,96,108,109] for a hydrological perspective. We just want to make sure that we provide results with enough details so that water managers can identify entire regions or a set of gauges that demand a deeper examination to assess a need for action.

Based on the results presented in Table 2, one can see that the proportion of gauges with significant trends at the national scale is larger for low flows (22–24%) than for high flows (14–16%). When we focus only on large and medium relative changes per decade ($\left|\widehat{\beta }\right|>5\%$), the proportion of significant gauges for low flows (20–23%) is even larger than that for high flows (12–15%). Moreover, at the national scale, low-flow indices seem to be of a greater concern, at least in terms of the magnitude of change, than high flows. This should not be interpreted as changes in high-flow indices cannot be a problem in some specific basins.

A deeper look into the results for the low-flow indices indicates that approximately two-thirds (63–70%) of the gauges with significant trends have a large relative change per decade ($\left|\widehat{\beta }\right|>10\%$). However, only approximately one-fourth (26–30%) of the gauges with significant trends have long record lengths (n > 60). As the record length is strongly associated with the uncertainties in the results, its value is important to understand the degree of confidence we have in these results. Results for low-flow indices presented in Table 2 show that about 30% of the gauges with significant trends have at least a medium relative change per decade ($\left|\widehat{\beta }\right|>5\%$) and long record length (n > 60). This proportion of gauges represents about 70 potentially problematic basins with the highest degree of confidence in the results due to the long length of record. If we include gauges with medium record length, the proportion increases from 30% to 50% of the gauges with significant trends, which represents about 130 basins. This statement does not mean that results obtained with short record lengths should not be used, but caution is required.

In addition, the results show that there is a tendency to have a larger proportion of gauges with a decreasing trend than with increasing trends for all indices, especially if the relative change per decade is large. However, this tendency is clearest for low flows. In general, the dominant trend in Brazil points toward drier conditions, although there are some regions that are becoming more humid, as we will show later.

The results from Figure 8 clearly show a large predominance of decreasing significant trends in SFR, PNB, TOA, and ALE hydrographic regions and to a lesser degree in ANC and ANO, regardless of the streamflow indices. All these regions are in northeastern and southeastern Brazil or encompasses part of the Central–Western and Northern parts. Conversely, we find the large majority of positive increasing trends in URU, ATS, and PRG hydrographic regions, all located in the south or the southern part of the Central–Western region. However, the degree of regional predominance is not as strong as the one for negative trends. This general spatial pattern agrees with results previously published for the southern [77,81,82,83,84,110] and for the northeastern and northern parts of Brazil [78,80,82,85].

Although most significant gauges in the AMZ region have increasing trends, except for Q7, for which there is a balance between increasing and decreasing trends, the proportion of gauges with significant trends is small. Finally, one observes a mix of increasing and decreasing significant trends for gauges located in the PRN hydrographic region, which encompasses the south and the southern part of Southeastern Brazil. A further discussion about the spatial distribution of these results is provided in the sequence.

The spatial distribution of the results for the streamflow indices discussed in the previous paragraphs is displayed in Figure 9 with the inclusion of information on the classes of record length and relative change per decade (see legend).

One can see that the majority of the significant trends have medium or large relative change per decade ($\left|\widehat{\beta }\right|>5\%$) with a few yellow or green gauges. Except for the Amazon hydrographic region, the country can be divided into two great regions by drawing an NW–SE line that cuts the PRN region in half. Above this line, one predominantly observes gauges with decreasing trends, including those located in the northern half of the PRN region. In the region below this line, which includes the PRG region, one predominantly observes gauges with increasing trends. This general pattern is true not only for the four indices presented in Figure 9 but also for the remaining four indices.

It is worth mentioning that the presence of a large number of gauges with long records and with large or medium relative increasing changes per decade trends for Q7 (dark and light blue circles) in Southern Brazil. The same pattern is observed for high-flow indices but with a smaller proportion of gauges. It is also possible to see a great number of gauges with long records and with large or medium relative decreasing change for Q7 in gauges located mostly in the ASD, northeastern part of PRN, and southern part of both SFR and ALE hydrographic regions. These regions and subregions that contain gauges with long record length, resulting in more precise trend results, and with large or medium relative changes per decade, putting them in a subset of potentially problematic regions, should have the attention of water resources managers, especially in cases of decreasing trends for Q7 and increasing trends for both QX1d and Qx30d.

3.3. Analysis with Gauges Unaffected by Reservoirs

The results presented so far were based on a dataset that contains gauges that are potentially affected by the presence of reservoirs upstream, a fact that may explain, at least partly, the findings of significant trends presented earlier. Although attribution is not the objective of this paper, we decided to provide a view of trend patterns for those gauges that are unaffected by reservoirs. To perform this second part of the analysis, we created a second dataset that contains only gauges that are unaffected by reservoirs.

The results of this analysis are presented in Figure 10, which displays the spatial distribution of gauges with significant trends for Q7 and QX1d. Results for two datasets are presented side by side to ease the comparison. Panels (a) and (c) provide the results for gauges that are unaffected by reservoirs, while panels (b) and (d) provide the results for all 1106 gauges.

In general, the spatial pattern obtained in the previous analysis for Q7 and QX1d using all 1106 gauges is still present in the analysis that uses only gauges unaffected by reservoirs. It is possible to say, for instance, that the large number of gauges with increasing trends for Q7 in Southern Brazil is not the result of reservoirs alone as we observe a clear pattern of increasing trends at unaffected gauges. Something similar can be said about the results for QX1d for the same region. In this case, the presence of reservoirs is expected to reduce the annual maximum flows, but we observe increasing trends at gauges both affected and unaffected by the presence of reservoirs.

Figure 10c shows that almost all gauges unaffected by reservoirs located in the SFR, ALE, and ANC regions show a decreasing significant trend for QX1d. This result supports the general conclusion that the decreasing trends for high flows observed in the first part of the analysis cannot be explained solely by the presence of reservoirs.

4. Conclusions

This paper used daily streamflow data from 1106 gauges across the Brazilian territory to assess the significance and estimate the magnitude of changes in eight extreme streamflow indices that cover different and important aspects of flows for water resources planning and management. Unlike most studies carried out in the past using Brazilian hydrological data, the effects of both serial correlation and the multiplicity of hypothesis tests in the regional analysis were assessed and properly considered so as to control the expected number of false positives.

Results on the proportion of gauges with significant trends obtained with three different approaches illustrated the distorting effects serial correlation and multiplicity may have in the analysis, showing how important it is to properly consider these issues in trend detection studies at local and reginal scales.

Trend detection analysis carried out for the 12 hydrographic regions of Brazil, stratified by record length and estimated relative change per decade, provides a comprehensive picture of the existence of trends in flows classified by different levels of confidence in the results, with results obtained for gauges with long record length being the most reliable. Results presented this way are believed to be helpful for water managers, so they can better judge the balance between uncertainty in the estimate and the magnitude of change, both important for making decisions regarding investments in further analysis and development of adaptation strategies.

The estimated proportion of gauges with significant trends in low and high flows was about 23% and 15%, respectively. About two-thirds of all gauges with significant trends for low flow, regardless of their record length, observed relative changes of more than 10% per decade. Results also showed that 50% of gauges with medium or long record lengths observed significant relative changes of at least 5% per decade.

The study identified a very clear spatial pattern. Regardless of the flow index, the majority of gauges with increasing trends were located in the southern part of Brazil, more specifically in the south of the PRN, PRG, and URU hydrographic regions, while the large majority of gauges with decreasing trends were found to the north of that area. A final analysis based on the concept of degree of regulation (DoR) shows that the spatial pattern of trends is also present when only gauges that are unaffected by reservoirs are used, suggesting that the presence of reservoirs alone cannot explain the increasing trends in high flows observed in the southern gauges or the increasing trends in the low flows observed in the gauges located in the northern hydrographic regions.

Further research is needed to investigate the role of potential natural and human-induced factors in the estimated changes in streamflow extremes. Climate variability and change, as well as human activities, such as agricultural expansion, urbanization, and construction of reservoirs, are all factors that must be examined in an attribution analysis that must be carried out in a more refined spatial scale.