Detecting Annual and Seasonal Hydrological Change Using Marginal Distributions of Daily Flows

: Changes in the hydrological regime are widely investigated using a variety of approaches. In this study, we assess changes in annual and seasonal ﬂ ow characteristics based on a probabilistic representation of the seasonal runo ﬀ regime at the daily time scale. The probabilistic seasonal runo ﬀ pa tt ern is constructed by determining quantiles from marginal distributions of daily ﬂ ows for each day within the year. By applying Fourier transformation on the statistics of the daily ﬂ ow partial series, we obtain smooth periodical functions of distribution parameters over the year and consequently of the quantiles. The main ﬁ ndings are based on the comparison of the dry, average, and wet hydrologic condition zones as de ﬁ ned by the daily ﬂ ow quantiles of selected probabilities. This analysis was conducted for ten catchments in Serbia by considering changes between two 30-year nonoverlapping periods, 1961–1990 and 1991–2020. It was found that the relative change in runo ﬀ volume is the most pronounced in the extreme dry condition zone in the winter season ( − 33% to 34%). The annual time shift is the largest in the dry and average condition zones, ranging from − 11 to 12 days. The applied methodology is not only applicable to the detection of hydrologic change, but could also be used in operational hydrology and extreme ﬂ ow studies via drought indices such as the Standardized Stream ﬂ ow Index.


Introduction
Hydrologic systems, and the water cycle in general, are subject to stresses and change caused by a range of drivers. The obvious or direct stressors on hydrologic systems include widespread land-cover change, urbanization, industrialization, and significant engineering interventions, while the indirect stressors are linked to the growing demands for drinking water, food, and energy for the population [1]. The major changes to the global hydrologic cycle over the last century are attributed to global warming and other influences of climate change in many regions of the world. The changes in frequency, intensity, and timing of precipitation directly contribute to modifications in the magnitude and timing of flow in rivers, including extreme floods and hydrologic droughts [1].
Assessment of expected change in the future hydrological regime is of interest for different projects and applications in water resources management and engineering, as well as in environmental studies. The aim and scope of the application determine the type of information needed from the assessment, as well as the baseline and future time horizons and spatial and temporal scales. The hydrologic change assessments may therefore differ significantly in aspects such as spatial extent (e.g., one location, catchment, region), assessment periods (past and/or future), timestep (e.g., daily, monthly, yearly), and hydrologic regime (mean-, high-, low-flow). In addition to the objectives of the assessment, data availability is a key factor in selecting the analysis technique.
Selection of the appropriate period for the assessment of any indicator of the hydrological regime may have a significant effect on the results because of the variability in hydrological regimes over longer periods. Climatological normal periods of 30 years often impose the use of 30-year-long periods for hydrological assessments in climate change impact studies. According to the World Meteorological Organization (WMO) [2], the 30-year normal period currently in use is acceptable for hydroclimatic applications. However, WMO [2] also warns that 30-year periods may not be acceptable for an analysis of extreme events (floods and low flows). The period of record is preferred for hydrologic and water resources engineering applications; for situational assessment and forecasting, a common period of record across the region of interest is preferred when using a record length less than 30 years is acceptable [2].
In Europe, the longest assessment period reported in the literature was related to the detection of long-term changes in the annual flow regime at the Ceatal Izmail hydrologic station on the Danube River over 1840-2015, which was the period of record for this station [3]. Renner et al. [4] used the 1930-2009 period for investigating 27 basins throughout Saxony/Germany, while the flood levels of the Danube at Novi Sad over the 1919-2007 period were analyzed in [5]. Stewart et al. [6] studied the period between 1948 and 2002 to detect changes toward earlier flow timing in a network of 302 western North American HSs, while the shortest period of twenty years (1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) was used by Laaha and Blöschl [7] for analyzing hydrographs from 325 catchments in Austria. The studies that focus on future conditions and projections of hydrologic change are mostly based on 30-year baseline catchment behavior, such as 1981-2010 [8,9], but also 1975-2005 [10], 1980-2012 [11], and 1985-2008 [12]. Having in mind long-term persistence and oscillations in the hydrologic regime, shorter assessment periods in hydrologic change studies may lead to false signals about the tendency and signs of changes [13,14]. This is particularly important when trend assessment is used as a change detection technique. According to Kundzewicz and Robson [15], a hydrologic time series of at least 50 years in length is required to distinguish between trends and variability in the hydrologic regime.
Identifying and assessing the changes in the hydrological regime usually implies looking for different trends or other types of changes in the historical hydrological data. Application of various statistical tests prevails in this group of approaches, such as tests for trend detection, homogeneity, serial correlation, etc. The use of the nonparametric Mann-Kendall test for trends has been extremely popular, especially after it was modified by Hamed and Rao [16] for autocorrelated data. Tests for abrupt changes, such as Pettitt's test [17], are also used often (e.g., [18]) and mainly for detection of human-induced changes to river regimes. The use of statistical tests for detecting changes in hydrological time series has also been criticized, mainly because their null hypotheses and assumptions may be compromised by the correlation structure of the hydrological time series [19]. While these tests are always numerically feasible, their outcomes may not be truly informative if their basic assumptions are violated. Furthermore, the outcomes of the trend tests depend heavily on the period of record and may indicate false tendencies characteristic of a part of a long-range dependence such as long-term oscillations [14,20]. Moreover, trends are rarely identified from data with a time scale finer than annual (e.g., [21]).
Seasonal runoff variation is one of the most important indicators of the hydrological regime, i.e., one of the most important hydrological signatures. Information on seasonal runoff distribution over the annual cycle is essential for water resources management and hydroecology because it provides predictable patterns of water availability [22]. While the mean seasonal pattern over the annual cycle describes the average hydrological regime at the location of interest, seasonal runoff also exhibits interannual variability that is also important for other water-and environment-related management. Therefore, detecting the change in the seasonal runoff regime is important for water management and for adaptation and mitigation strategies.
Seasonal runoff variation is typically described over the annual cycle at the monthly time scale and less often at the daily time scale. Daily streamflow data have several advantages over monthly data: they provide a more functional representation of the runoff regime [23] by offering important insights for water resources management, including linking hydrological regimes to habitats and biotic communities [24], introducing variable hydrologic drought thresholds [25], and predicting streamflow characteristics in ungauged basins [12]. Most runoff indicators of interest for studying aquatic and riparian ecological response (magnitude, frequency, duration, timing, and rate of change [24]) can be assessed from the daily representation of the annual and seasonal runoff regimes. Daily flows also facilitate estimating different measures of runoff seasonality. The annual and seasonal changes in position of the hydrograph centroid are frequently considered [6,10,26], while the floods and droughts are investigated through onset date, end date [6,27], or duration [3].
Variability in the seasonal runoff pattern from year to year is typically described by computing and drawing lines around the mean (or median) pattern that correspond to seasonal runoff of certain non-exceedance probability. These lines may be considered as boundaries between the zones of different hydrological conditions such as average, wet, or dry conditions. There are no general recommendations in the literature about specific probabilities that could serve to define the boundaries between the hydrological conditions. These boundaries are frequently associated with various threshold levels that are used in water resources management to define certain design or reference quantities. In drought-related studies, the threshold levels used to identify drought events can either represent water demand or the boundary between normal and low-flow conditions [28]. Although some applications may require fixed threshold levels, the thresholdsʹ variation throughout the year reflects seasonally different hydrological conditions, ecological requirements, or water demands. For instance, in [23] the 10th percentile of monthly flows delineates very dry and dry hydrological conditions, the 25th and 75th percentiles define the zone of average hydrological conditions, and wet conditions are those above the 75th percentile. Daily flows that exceed the 75th percentile are also classified as high flows by Pekarova et al. [3], but those below the 50th percentile are classified as low flows, with additional separation into extreme low flow and low flow categories.
The probabilistic approach to defining time-varying threshold levels is also embedded in the methodology for computing the Standardized Precipitation Index (SPI) [29], where monthly precipitation aggregated over scales of 1, 2, 3, …, 12 months for each calendar month is transformed into a standard normal quantile, SPI, based on non-exceedance probability of aggregated precipitation. Typically, a threshold of SPI = 0 is used to define meteorological droughts as periods with less accumulated precipitation than expected in a long-term period. Similarly, the Standardized Streamflow Index (SSI) [30], as a probabilistic hydrologic drought index based on monthly flows, was introduced. Computation of both the SPI and the SSI involves an intensive procedure of determining the type of probability distribution for each aggregation scale [31,32]. The standardization of precipitation and flows in this probabilistic approach explicitly introduces fixed SPI or SSI thresholds (e.g., −1, 0, 1) corresponding to certain probabilities, which translate into variable thresholds for precipitation or runoff volume in each calendar month.
Representing the annual runoff cycle at the daily time scale has its advantages and disadvantages. The advantage is that the daily resolution provides more detail and smoother transitions compared with the monthly scale [33]. On the other hand, the disadvantage is that the calculation is much more demanding and extensive because it includes analyzing and fitting 365 probability distributions [33]. Furthermore, using the daily step requires smoothing to avoid jumps in consecutive days [25,33,34].
Considering all the above, the main issues in detecting annual and seasonal hydrological change are threefold. First, the analysis of changes, and particularly trend analysis, mainly deals with the annual scale while the changes in seasonal patterns are seldom addressed. In many cases, the change is not detectable on the annual scale, but is obvious on the seasonal (monthly) scale. Second, the seasonal pattern is typically described in terms of monthly flows, which may not be useful for some applications and creates an abrupt transition in flows from one month to another. Third, the change in the seasonal runoff pattern is usually identified from the average pattern, and the year-to-year variability in the pattern is typically neglected. In cases when this is considered, it is assessed from empirical distributions of monthly flows. Moreover, the change in hydrological conditions with the assessment period has not been analyzed to our best knowledge.
To address the above identified issues, the goals of this paper are as follows: (i) to investigate seasonal change despite the lack of significant change on the annual scale in our study area; (ii) to utilize probabilistic representation of the seasonal runoff cycle at the daily time scale; and (iii) to detect changes in zones of different hydrological conditions. To do so, we use a probabilistic approach for defining the seasonal runoff pattern at the daily time scale by fitting 365 marginal distributions of daily flows (MDDFs). The parameters for the 365 marginal distributions are modeled as periodic functions to avoid abrupt changes in quantiles in consecutive days. The estimated quantiles from the 365 MDDFs result in smooth lines in the plot of the annual runoff cycle and represent the boundaries of the zones with flows for the given probability of occurrence. To detect changes in the seasonal runoff pattern, we compare wetness (hydrological) conditions as defined by probabilistic thresholds between two nonoverlapping 30-year subperiods. The two subperiods are compared in terms of the overall seasonal runoff pattern, runoff volume, and timing of the annual and seasonal runoff. The analysis is performed for the selected catchments in Serbia for which the statistical tests do not reveal any change or trend at the annual scale in the 1961-2020 period. Therefore, the purpose of this research is to investigate potential alterations in the seasonal runoff pattern despite the lack of significant changes on the annual scale.
This paper is organized as follows. Section 2 describes the study area and data sets, while Section 3 presents the methodological framework used to construct a probabilistic annual runoff cycle using MDDFs, define hydrological condition zones, and evaluate annual and seasonal hydrologic changes per zones. Section 4 shows the results of distribution fitting, zoning, and the detected changes. Section 5 discusses the results from two aspects, first by comparing them with the previous findings on hydrologic change in the study area, and then by highlighting the benefits of the MDDF approach for hydrological applications. Section 6 gives the key conclusions of the presented work.

Study Area and Data
This study was conducted on daily flows in the period 1961-2020 at ten hydrological stations that belong to the Danube River basin in the Republic of Serbia ( Figure 1) and are operated by the Republic Hydrometeorological Service of Serbia (RHMSS). Basic information about the ten stations and their drainage areas is given in Table 1.  The daily flow series were checked for data completeness and only minor gaps were found that were filled either through interpolation or through regression analysis using data from the neighboring stations as predictors. The series homogeneity and randomness, as well as the presence of trends, were tested on the annual flow series with a range of parametric and nonparametric tests. The results of all tests (p-values) are given in Table  A1 in Appendix A, showing that the null hypothesis on homogeneity, randomness, and absence of trends cannot be rejected at a significance level of 0.05 for all stations.

Methodology
Detection of changes in hydrological regime at selected stations follows the methodology organized in the workflow consisting of four steps and illustrated in Figure 2. The first three steps (a-c) result in probabilistic representation of the annual runoff cycle at daily time scale, while the fourth step (d) is dedicated to defining hydrologic condition zones, calculating indicators, and detecting the annual and seasonal change. The workflow in Figure 2 relates to the analysis for one station and is repeated for each one. The analysis is performed for the 1961-2020 period of the available data, and for two subperiods, 1961-1990 and 1991-2020, representing two standard 30-year periods used by the WMO [2]. The core of the proposed methodological approach is the use of daily flows for construction of the probabilistic annual runoff cycle. In the first step (a), daily flows are prepared for analysis by defining the data series for each day within a year, resulting in a total of 365 data series for which the corresponding statistics are computed. This step is described in Section 3.1. In the second step (b), statistics of the partial daily series are modeled as the periodic functions to achieve their smooth representation over the year. This step is described in Section 3.2. In the third step (c), the periodic daily statistics are transformed into periodic parameters of marginal distributions of daily flows (MDDFs). The quantiles of MDDFs for different probabilities are then estimated for each day, and the diagrams of these quantiles are created to represent probabilistic seasonal runoff pattern. The lines of MDDF quantiles are smooth because an implicit, anterior smoothing is achieved by using periodic probability distribution parameters, thus avoiding sudden changes in flow between two successive days. The third step is described in Section 3.3. The fourth step (d) of the methodology comprises the definition of the zones of different hydrologic conditions (Section 3.4) and detection of changes in hydrological regime from one period to another (Section 3.5). The changes in the seasonal pattern, runoff volumes, and runoff timing at annual and seasonal scales are considered.

Daily Flows as a Stochastic Process
The continuous stochastic process {xt; t ≥ 0} is mathematically formulated [35] by the time function xt = f(t; α, β, γ,…), where t is time, and α, β, γ, etc. are the parameters of a multidimensional distribution describing the temporal structure of the process. Hydrologic daily time series are a close approximation to the continuous stochastic processes, where the process xt is represented by the observed daily flows xτ,i (with i = 1, 2, …, Nnumber of years, and τ = 1, 2, …, 365-number of days). Hence, the xτ,i series is a set of N process realizations over the [0, 365 days] interval.
The observed daily flow series xτ,i at each station is arranged in a matrix, in which each row represents one date (i.e., the ordinal number of a day) within a year, and each column represents one year: The data are arranged in the matrix (1) by hydrologic years, starting on 1 October (τ = 1) and ending on 30 September (τ = 365). In the leap years, 29 February is omitted, similar to many studies (e.g., [25,36]).
For each of the 365 dates, i.e., for each row in the matrix (1), the sample statistics are estimated using the method of moments for both the original series (x) and its logarithmic transformation (y = log x). In this way, column vectors of the statistics with 365 rows are obtained for each day in the year. We will refer to these series of statistics as daily statistics. They include the mean mτ, standard deviation sτ, and coefficient of skewness Csτ (τ = 1, 2,..., 365).

Periodicity Analysis
The series of daily means , = 1,2, … ,365, represents the mean seasonal runoff pattern. Due to variability in daily flows, the annual cycle of may not resemble a periodic-like function that would be obtained in a similar process with monthly data. The same is valid for the other daily statistics as well. Therefore, the seasonal runoff pattern can be described by introducing the periodic functions of the mean, standard deviation, or other parameters. These periodic functions are estimated from the column vectors of the statistics of the observed daily flow series. For any parameter υ, its periodic component υτ,per may be expressed using the Fourier series [35,37]: where υ̅ τ is the mean of υτ, h is the number of significant harmonics, Aj and Bj are the Fourier coefficients, fj = ј/ω is the frequency of the j-th harmonic, ω is the base period of 365 days, and τ = 1, 2,..., 365. The Fourier coefficients are calculated as follows: while the amplitude of the j-th harmonic is computed from the Fourier coefficients: Modeling the periodicity of the daily flow series includes finding significant harmonics for the periodic component shown in Equation (2) [35]. The number of significant harmonics adopted in this study is three, as identified in the previous study for the same stations [38] based on the periodicity analysis for the mean via the periodograms, as suggested by Yevjevich [35].

Marginal Distributions
The row vectors of the process xt in matrix (1), i.e., the vectors , , . . . , , = 1,2, . . .365, represent N realizations of the process for day in the annual cycle. The process in matrix (1) can therefore be seen as a 365-dimensional random variable. For each day, based on vector (6), a one-dimensional or marginal distribution is identified representing the distribution law of the process at the point in time. A family of 365 distributions defines a multidimensional distribution of the observed stochastic process. The series of daily flows is nonstationary because of the intra-annual periodicity of its parameters. Therefore, the parameters of the 365 marginal distributions also change periodically during the year because they are estimated from the periodic daily statistics.
There are several approaches for identifying marginal distribution functions [35]. The simplest approach that is suitable for extensive multisite computations at regional level is to apply a unique theoretical distribution type for all marginal distributions, while allowing its parameters to change periodically during the year. Radić and Mihailović [38] have shown that the log-Pearson type 3 (LPT3) distribution is the most appropriate choice among the two-and three-parameter distributions in the studied region of Serbia not only according to the statistical tests, but also because it satisfies two conditions related to the nature of the modeled physical process: (i) the upper bound of the marginal distribution should be higher than the observed maximum flows; and (ii) the lower bound of the marginal distribution should not be less than zero.
The flow quantiles Qpτ of a non-exceedance probability p = p(x) are estimated from the 365 fitted LPT3 marginal distributions. The probability density function of LPT3 distribution is given with where , , and are the shape, scale, and location parameters of LPT3 distribution, respectively, and ( ) denotes a gamma function. If the natural logarithm of data is used ( = ln ), then = and = 1. In this study, we use the common (decadic) logarithm with = 10 and = 0.434. The LPT3 parameters are estimated using the method of moments with log-transformed data as where , , and are the mean, standard deviation, and coefficient of skewness of the log-transformed observed data = , respectively. The bounds of LP3 distribution depend on the sign of the shape parameter :

Hydrological Condition Zones
The family of 365 MDDFs allows us to understand how the daily flows at some locations can change within the annual cycle. It also provides the basis for constructing the intervals in which the daily flows can occur with certain probability, and thus it could be used to define an interval, or a zone, of typical, "normal" conditions at different times in the year. The zone is bounded by the specific lower and upper quantiles of MDDFs. The zones outside the "normal" condition zone are characterized by deviations of daily flows from the typical regime for a given day [39].
Hydrological condition zones are defined here for specific probabilities of occurrence of daily flows in the annual cycle. The zone thresholds, i.e., boundaries, are quantiles of specific probabilities (Figure 3a), selected to represent dry, average, and wet conditions, based on the previous research [3,7,23,36,40]. The average condition zone is the interval between the 0.3 and 0.7 quantiles. This zone is divided into two subzones by the median flow, p(x) = 0.5, an important hydrological indicator recommended by WMO as a measure of central tendency when undertaking assessments that focus on characterizing typical conditions [2]. The 0.05 quantile delineates the dry condition zone from the lowest zone representing extremely dry conditions, while the 0.99 quantile delineates the wet condition zone from the upper zone representing extremely wet conditions.

Estimating Change in Hydrological Regime
The changes in seasonal runoff patterns are first examined visually by comparing the average zone and the median threshold (Figure 3a), as well as the periodical function of the mean daily flows. The changes in the annual and seasonal runoff volumes and timing between the two periods are computed by examining the area below the upper threshold of the hydrologic condition zone. The runoff volume is obtained by integrating this area (i.e., the flows representing specific MDDF quantiles), while the timing is described by the time coordinate of the centroid of this area [6,26]. The following two indicators of change are computed at annual and seasonal scales: 1. Relative change in runoff volume, ∆Vp (%): where Vp is runoff volume below the p-th MDDF quantile representing the upper threshold of the zone for the given period.
where Cp is the time coordinate (ordinal number of day within year) of the centroid of the area under the p-th MDDF quantile for the given period.
The seasons within the hydrological year used for evaluation of the changes are defined as follows: October, November, and December represent autumn, January, February, and March-winter, April, May, and June-spring, and July, August, and September-summer.

Daily Flow Statistics and Their Periodicity
The daily statistics computed for each day within a year include the mean (mτ), standard deviation (sτ), and coefficient of skewness (csτ). Figure 4 shows the seasonal variation in the daily statistics for HS3 in all three considered periods, together with their corresponding periodic functions (smooth lines). The mean and standard deviation both exhibit distinct seasonal patterns which agree in phase. The highs occur in the spring and the lows in the summer (Figure 4a,b). This pattern is characteristic for all but two considered stations (HS1 and HS6). The skewness does not exhibit a distinct seasonal pattern like the mean and standard deviation, but a complex one with more highs and lows. It is also worth noting that there is a significant positive skewness over the entire annual cycle at all stations, as shown for HS3 in Figure 4c.   Figure 5 shows the daily statistics of the log-transformed data at HS3: mean (myτ), standard deviation (syτ), and skewness (csyτ), as well as their corresponding periodic functions. Expectedly, logarithmic transformation suppresses the variability in the daily statistics compared with the ones in the original space. The distinct seasonal pattern is still visible for the mean, but not for the standard deviation, for which the amplitude of the intra-annual oscillation is significantly smaller. The logarithmic transformation reduced the skewness and even produced negative values in both the daily skewness series and the estimated periodic function during certain parts of the year. This was noticed at all stations in some part of the year for at least one period considered. This further affects the parameters and the shape of the LPT3 marginal distributions. The statistics and their periodic functions for all stations are given in the Supplementary Materials (Figures S1 and  S2).

Periodic Parameters of the Marginal LPT3 Distributions
The parameters of marginal LPT3 distributions are periodic functions, resulting from the periodic statistics myτ,per, syτ,per, and csyτ,per. Figure 6 shows the seasonal variation in the three LPT3 distribution parameters at HS3. The LPT3 parameters highly depend on the value and sign of the skew, csy. Therefore, a comparison of periodical functions of statistics ( Figure 5) with LPT3 distribution parameters ( Figure 6) expectedly shows that the differences between studied periods are transferred from the statistics (especially the skew) to the distribution parameters in accordance with the expressions in Equation (8). The results for all stations are given in the Supplementary Materials ( Figure S3).

Marginal Distributions of Daily Flows
The flow quantiles Qpτ of a non-exceedance probability p = p(x) are estimated from the 365 fitted LPT3 distributions, one per each day within a year, with its periodic parameters. The smooth lines shown in Figure 7 represent the quantiles of MDDFs for the three periods considered at HS3. The results for all stations are given in the Supplementary Materials ( Figure S4).  Figure 5c), marginal LPT3 distributions are bounded from above (Equation (9), also, e.g., [42]). The values of the upper bound were computed for all such cases and it was found that they do not cause underestimation of daily flow upper-tail quantiles for any non-exceedance probability of interest in hydrologic applications.
Visual inspection of the MDDFs in Figure 7 reveals that the differences between the three periods are the smallest for the median daily flows and in the average zone, and grow larger for the upper and lower quantiles. For this station, HS3, the maximum quantiles for 1991-2020 are higher and occur later compared with those for 1961-1990. For the joint 1961-2020 period, the quantile lines lie expectedly between the lines of the two subperiods.

The Zones of Hydrological Conditions
The zones of hydrological conditions for all stations are shown in Figure 8 in accordance with the adopted classification.
The patterns of thresholds (quantile lines) in Figure 8 obviously change with the period. The changes occurring in the wet condition zone at all stations are more visible due to the scale, but the changes in the dry condition zone are also significant, as can be seen in Figure 7b. At all stations, the pattern remains the same for quantiles 0.3, 0.5, and 0.7, which define the normal hydrological conditions. The seasonal pattern is not changed at any station from a simple unimodal shape to a mixed bimodal shape, or vice versa. The bimodal shape of the seasonal pattern is the most visible at HS1, where the most interesting change is for the 0.99 quantiles. Here, the primary peak in 1961-1990 was in the autumn while the secondary one was in the spring; in 1991-2020, the primary peak was in spring, while the secondary was in winter.

Annual and Seasonal Hydrologic Condition Changes
The graphical presentation of hydrologic condition zones defined in Figure 3a is helpful for detecting changes in daily flow patterns, while the quantification of change in respect to runoff timing and volume (the areas defined in Figure 3b) is required for practical purposes. Change in runoff timing is assessed from the shifts of the centroids of the areas below the upper thresholds of the hydrologic condition zones, as shown in Figure  9a for the annual scale and in Figure 9b for the seasonal scale at HS3. Relative changes in runoff volume and time shifts of the centroids, calculated from Equations (10) and (11), are shown in Table 2 and Table 3, respectively, for all stations. Table 2. Relative change in runoff volume (%) in the recent 1991-2020 period compared with the 1961-1990 period. The first column indicates the probabilities p(x) of the upper thresholds of the hydrologic condition zones. Changes are shown for annual (Ann) and seasonal scales (Autumn, Winter, Spring, and Summer), and colored in blue (increase) and red (decrease).

p(x) Season
Legend: Δt (days) A general impression from Table 2 is that there is less runoff volume in the recent 30year period compared with the previous one. HS6 is the only station at which runoff volume has increased in the recent period. HS2 and HS9 have distinct increases in extremely dry conditions in the winter and spring seasons, respectively. Additionally, HS2 exhibits a significant runoff volume increase in the wet condition zone in summer, as does HS10. The most unfavorable changes are found in HS10, with significant volume reductions in dry condition zones both annually and seasonally, and in all hydrologic condition zones in the spring. This station is located approximately 90 km downstream of the outflow of the Pirot hydropower plant (HPP) that started operation in 1990, at the very end of the first 30-year study period, as a peaking HPP [43]. It should be noted that we did not detect a change in the mean annual flow series of HS10 in the data preparation phase. This is most likely due to the small effect of this HPP and its reservoir on the regime of the downstream HS10 for several reasons: the small drainage area of the reservoir compared with the drainage area of HS10, the annual water balancing mode of the reservoir, and subdaily disturbance of the regime at HS10 by HPP outflows.
The shift in runoff timing shown in Table 3 is less pronounced compared with the relative volume change. The changes in more than 7 or 10 days might be significant for some applications in the water sector such as irrigation and hydropower operation planning. In that respect, the most significant changes are found at the annual scale in the dry and average condition zones in HS1 and HS2 of the Sava River basin, where the runoff centroids occur up to 11 and 8 days earlier, respectively. A later occurrence of dry conditions is found at the annual scale in HS7, HS8, and HS10, which are all situated in the upper and middle parts of the Južna Morava River basin. Wet hydrologic conditions occur 7 days later in HS2 and HS7 in summer, and 9 and 8 days earlier in HS7 and HS8, respectively, in the spring season.
The data from Tables 2 and 3 are shown as boxplots in Figure 10 and Figure 11, respectively. The results are organized in such a way that annual and seasonal changes in the hydrologic condition zones are easily grasped for the whole study area.  A general decreasing tendency in annual runoff volume in 1991-2020 compared with 1961-1990 at all stations is clearly visible from the boxplot in Figure 10. The most pronounced decrease is for the wet condition zone. Runoff volume has also decreased in all seasons and for all hydrological condition zones, except the autumn runoff in the average zone (for the 0.5 and 0.7 quantiles), and the spring runoff in the extremely dry zone. There is less water available in wet zones (0.7 to 0.99 quantiles) during the high-flow season in winter and spring. Similarly, in dry zones (0.3 and 0.05 quantiles) there is less water available in the summer and autumn seasons in the recent period 1991-2020.
In terms of runoff timing at the annual level (Figure 11), the centroids for dry and average conditions occur earlier in 1991-2020 compared with 1961-1990, while the centroids for wet conditions occur somewhat later. The seasonal changes in timing are not pronounced for all zones, but the centroids for wet conditions appear earlier in spring at all stations, and later in winter at most stations.
The annual changes in both runoff volume and runoff timing at all stations are mapped for spatial analysis in Figure 12 and Figure 13, respectively.  The spatial distribution of the relative volume change between the two periods shown in Figure 12 reveals that the catchments in western Serbia have experienced an increase in runoff volume in the recent period for the dry condition zone. It should be noted that these catchments are small to medium areas in size. The decrease in annual runoff volume is more pronounced at the transition from the average to wet condition zone in all catchments. The most western station, HS1, behaves differently compared with the other stations. It has the largest catchment area and represents the only station with runoff generated almost completely out of the country. This difference in behavior is pronounced in respect to runoff timing, as shown in Figure 13. The runoff timing for extreme dry conditions is later in the recent period at the annual level for all other stations. For the average wetness conditions, runoff generally occurs earlier, but this change is insignificant, while it generally occurs a bit later for wet conditions. The most southern station, HS7, shows earlier runoff timing for wet conditions. Expectedly, the timing shifts at the annual level are more pronounced compared with the shifts within seasons.

Long-Term Changes in Hydrological Regime
The majority of previous studies considering the same area used trend detection of annual/seasonal/extreme flows to study hydrologic change in the past or in the future, or the "climate-hydrology-assessment chain" approach [44], in which a selected climate scenario plays an important role. Comparison of the approach presented here with other studies can therefore only be made indirectly, i.e., based on the detected changes. The methodology presented here can be used to analyze the seasonal runoff pattern, but it is less convenient for the analyses of phenomena at short, event-based time scales (such as flow flashiness) or at time scales longer than one year (such as interannual variability), as underlined by Blum et al. [36]. The probabilistic seasonal runoff pattern obtained from the MDDFs is used here to define zones of different hydrologic conditions and to detect a change in this pattern by looking into the changes in runoff timing and volume between the two 30-year periods. Our study is therefore aimed at detecting a structural change in the hydrological regime instead of the simple detection of trends or projections via climate change scenarios. The proposed approach is also not aimed at detecting the changes in intensity and frequency of flood events, which are often considered important consequences of climate change. However, our approach can be useful for a general understanding of the hydrological regime at a given location and can provide indications of the changes in the high-water regime, which is related to the floods.
In general, our results are consistent with the findings of Blöschl et al. [45], who reported decreasing flows in the Balkan region, and with Lobanova et al. [46], who used hydrologic simulations with a daily timestep under future climate conditions and concluded that shifts in runoff seasonality, particularly earlier occurrence of high flows in the snowmelt-driven catchments, can be expected. Also, according to [46], winter and early spring runoff in Central and Eastern Europe (in the high-flow season) are expected to increase. This is fairly consistent with our results (Figure 9), which show an increase in runoff in the wet season over 1991-2020 that is not significant (yet), but later timing of winter runoff and earlier timing of spring runoff ( Figure 10) are detected, indicating the shifts in the wet period projected for the future in [46].
In a comprehensive analysis of trends in mean annual and seasonal flows in Serbia over 1961-2010 [47], no trend in mean annual flows was detected at all stations used in our study except HS8 and HS10, for which moderate decreasing trends were found. This is not in agreement with our findings from the data preparation phase (Table A1 in Appendix A) on the absence of trends over the 1961-2020 period at all stations considered. On the other hand, we found the greatest negative changes in annual runoff volume at HS8 and HS10 for the 0.5 quantile zone (−10% and −15%, respectively; see Table 2). In the same zone, HS6 exhibits the highest positive change in annual runoff volume of 10%, while there was no trend detected at this station in [47]. These discrepancies may stem from the difference in the periods for analysis (1961-2010 in [47] vs. 1961-2020 in our study), because the additional 10 years in our study cover a wetter period in Serbia, including an extremely wet year, 2014 (also known for extreme floods). On the seasonal level, a low-significance positive trend is found in [47] in autumn and winter runoff at HS6, whereas we found changes in seasonal runoff volume of 11% and 15% for autumn and winter, respectively. Additionally, a moderate negative trend in spring runoff at HS10 was found in [47], while our results show a change in spring runoff volume of -25%. However, we detected several considerable seasonal changes exceeding 10%, both positive and negative, in the average condition zone, while these were not detected as trends in [47].
A variety of studies about different aspects of hydrologic change for the Sava River (e.g., [14,[48][49][50]) are in line with our results for HS1, showing that the mean flow is slightly declining. Hydrologic simulations in [49] for the near future (2011-2041), which partially overlaps with the second period in our study, indicate changes consistent with our results about increases in winter runoff volume in all hydrologic condition zones. An increase in winter runoff at HS1 is also found as a weak (statistically insignificant) positive trend over the long period of 1928-2017 [48]. A substantial decrease in river flows expected in the spring and summer seasons [49] is also visible in our results for all zones except for wet condition zones. The analysis of long-term trends of climate drivers and assessment of runoff in [50] also shows declining runoff in the Sava River basin as a consequence of increased air temperature and evapotranspiration. General findings in [13] for HS1 and HS4 indicate decreasing trends in summer and autumn runoff, and mostly increasing trends in winter and spring runoff, which are partially in accordance with our results. Our study indicates that HS1 and HS4 exhibit decreases in spring runoff volume and an increase in autumn runoff volume for the 0.5 quantile zone.
The hydrologic projections for the Kolubara (HS2) and Toplica (HS9) catchments indicate a decrease in snow storage and a substantial decrease in runoff in 2001-2030 (near future) compared with 1961-1990 (baseline period) [51]. Our results do not show a substantial but rather a slight decrease in runoff volume in the average and wet zones, which is closer to the findings of Idrizović et al. [52] for HS9. In [51], the largest decrease in runoff at HS9 is associated with the spring season, which is in accordance with our results for the wet zone only. The flood flow series simulated in [51] indicate little change for HS9 and a small increase for HS2 in median annual maximum floods. Our results for the wet zone at the annual scale for HS9 are consistent with these results, and also with the results in [52]. The analyses in [51] additionally show that the future droughts are expected to be more frequent, to last longer, and to start earlier in the summer for both the HS2 and HS9 catchments. Judging by the combination of change in runoff volume and time shifts in the dry and extremely dry zones (Tables 2 and 3, 0.05 and 0.3 quantile), runoff reduction in these zones is already evident, but the time shift toward an earlier date can be seen in the extremely dry zone at HS9 only.

Probabilistic Annual Runoff Cycle as an Indicator of Hydrological Conditions
The probabilistic annual runoff cycle obtained from the MDDFs is used in this study to assess the changes over time in hydrological regime at selected rivers and locations in Serbia. It can also be used for other purposes, including operational monitoring of the hydrological regime aimed at early warnings and water management, and especially in operational reservoir management and management of droughts. The diagram with the MDDF of the annual runoff cycle (Figure 7) shows the zones of different hydrological conditions. Plotting actual hydrological data over this diagram clearly indicates current conditions in the catchment, but it can also reveal departures from normal conditions and anomalies. It can therefore be a support for responding to unfavorable situations and their mitigation.
The zones of hydrological conditions in the probabilistic annual runoff cycle are delimited by different probabilistic thresholds, which vary over the annual cycle because they represent the quantiles of MDDFs. This type of zoning based on probabilistic thresholds is very similar to the application of the SPI and SSI as drought indicators [29,30]. The main rationale behind the SPI and SSI as nondimensional (standardized) probabilistic indicators is to enable comparative analyses of precipitation or streamflow at different locations and at different times. Application of these two indices in monitoring hydrological conditions therefore requires the transformation of current precipitation or streamflow observations into standardized values. Contrary to this, the probabilistic annual runoff cycle allows direct evaluation of current streamflow observations for a given location.
A common characteristic of the methodologies for assessing the probabilistic annual runoff cycle and the SPI/SSI is related to fitting probability distributions to the observations. The SPI and SSI are usually calculated from monthly data over different accumulation periods from 1 to 12 months, whereas precipitation or streamflow accumulations for each month in a year are fitted by separate parametric or nonparametric distributions. Stagge et al. [33] attempted to fit a set of 365 parametric probability distributions for daily SPI estimation, but they retained the monthly averaging period. When compared with the SPI/SSI applications, using the daily time scale for the probabilistic annual runoff cycle provides information on runoff conditions with finer resolution and helps to avoid discontinuities when using the monthly scale and transiting from one month to another.
Regardless of temporal resolution, finding appropriate distribution type(s) for precipitation or streamflow is a crucial issue because it can considerably affect the computed SPI/SSI series and the derived drought characteristics [30,53]. Inappropriate distribution can lead to inconsistencies in the SPI or SSI such as the average not being equal to 0 and the variance not being equal to 1 [30]. When assessing the probabilistic annual runoff cycle for a given location, finding the most appropriate distribution is still a great challenge, but the targets of this process are less demanding regarding spatial comparability and in terms of constraints related to statistical properties.
Furthermore, the probabilistic annual runoff cycle is obtained here using the periodic functions for statistical properties of daily streamflow. This approach leads to smooth periodic-like lines of streamflow quantiles within a year and additionally ensures continuity of the thresholds that delineate certain hydrological condition zones.
As the SPI and SSI are devised as indices that should facilitate comparability of conditions in space and time, relevant thresholds that define specific conditions (e.g., droughts) are obtained via the equiprobability transformation from a selected distribution into the standardized normal quantiles. Typical thresholds are SPI = 0 or SSI = 0, and they are constant over the annual cycle. Conversely, the thresholds in the probabilistic annual runoff cycle are expressed directly in the physical units of streamflow at a given location and are therefore variable throughout a year. This provides a more intuitive presentation of the hydrological conditions in the annual cycle, and the zones delimited by these thresholds readily show deficits and surpluses in available water quantities.

Conclusions
Hydrologic change detection involves different identification approaches, temporal scales, and spatial scales. This study presents an analysis of changes in the annual and seasonal runoff of ten catchments in Serbia in the 1961-2020 period at the daily time scale. The analysis is based on the probabilistic representation of the seasonal runoff pattern that is obtained from marginal distributions of daily flows (MDDFs) with periodic distribution parameters inherited from the periodic statistical parameters of daily flow series. The hydrologic (wetness) conditions are defined as the zones between the MDDF quantiles of specific probabilities, selected to depict dry, average, and wet normal conditions in the catchment.
Based on a comparison of these zones between the two 30-year periods at the annual and seasonal scales, the following may be concluded: 1. The seasonal runoff pattern changed from one period to another in terms of temporal shift and the occurrence of more extreme flows. However, the general pattern of seasonal runoff remained the same. The prevailing pattern is simple and unimodal, while the less present mixed regime is bimodal. 2. In most of the catchments, runoff volume has decreased in the recent 1991-2020 period at both the annual and seasonal scales. The critical season is summer for dry and average conditions, with volume reduction in all catchments. 3. The most pronounced shift in runoff timing is found on the annual scale. Dry and average conditions occur earlier at this scale. The change in runoff timing is found to be insignificant for all seasons and zones, except for wet conditions, which occur earlier in spring.
The proposed zoning of hydrologic conditions based on probabilistic variable daily flow thresholds enables more precise analysis of runoff volume and timing between the selected periods. It also enables more detailed analysis of zones that are critical for some applications in a certain period of the year. For instance, such an analysis may be of interest in drought-related studies in the seasons when water deficit conflicts with water demand.
The results of this study are generally in line with the results of the previous hydrologic change analyses in the study area, in terms of both past changes and those predicted for the near-future climate, but our results may be considered more detailed due to the applied approach based on MDDFs with periodic distribution parameters. It is demonstrated that the proposed approach allows for detecting a structural change in the hydrological regime, and detecting changes in the regulated flow regime when statistical tests cannot detect the change. The daily time scale of the seasonal runoff pattern representation makes this approach more useful for monitoring actual catchment wetness and runoff conditions and therefore for operational water management and early warning systems. In situational assessment and forecasting, the wetness condition zones may be used for flow forecasting by focusing on the zone to which current daily flows belong. It is worth noting that the boundaries of the wetness zones could be changed easily depending on the purpose required. Furthermore, the zone boundaries, i.e., the probabilistic thresholds, can also be converted to hydrologic drought indicators such as the SSI.
Supplementary Materials: The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/w15162919/s1. Figure   variances (Fisher's F-test and Levene's test), the Mann-Whitney test for homogeneity, the Wald-Wolfowitz run test for independence, and the Mann-Kendall test for trends (with the null hypothesis that there is no trend in the data).