Data
The data sources used in this study are described in Table 1. We used catchment boundaries and streamflow time series for 8255 stations, retrieved from the Global Streamflow Indices and Metadata Archive (GSIM62,63). The GSIM dataset provides streamflow time series at monthly, seasonal, and yearly resolutions. In this study, we used time-series indices at monthly resolution, and we analysed only stations with monthly maximum and mean daily streamflow observations for at least 20 years and high catchment-delineation quality (attribute in the GSIM dataset). These observations cover different periods between 1900–2016 (Supplementary Fig. 14). For the selected catchments, we derived time series of hydro-meteorological and biophysical variables. These include catchment averages of precipitation (P) from MSWEP64, median soil moisture (SM) near the surface (0–5 cm depth) and the root zone (0–250 cm depth) from GLEAM65, the total water storage (TWS) from GRACE66 and the normalized difference vegetation index (NDVI) from NOAA67. We removed urban and cropland areas68 (retrieved from https://ghsl.jrc.ec.europa.eu/download.php?ds=bu; https://dataverse.harvard.edu/dataverse/geospatial_human_dimensions_data) from the NDVI maps to focus on ecological drought impact signals in the data. Additionally, we retrieved surface water extent time series from Donchyts et al., 201669, to analyse both lakes with an area between 1 and 2 km2 and water supply and irrigation reservoirs (no area restriction) located within the investigated catchments. The restriction in the minimum and maximum lake area is the result of a trade-off between a manageable number of lakes to be analysed and a widespread presence of these surface water resources in the investigated catchments. Finally, we only analysed the time series of the reservoirs and/or lakes with the largest volumes within each catchment. All data were extracted from satellite products or in-situ observations, the only exception being soil moisture of the root zone for which only modelled data are available on a global scale, but satellite soil moisture data are assimilated in the GLEAM model. For each station, we used the full extent of available time series data (time periods varying between 1900–2016) to compute streamflow statistics for determining benchmarks (i.e., threshold analysis). For the analysis of the drought propagation and the detection of drought and flood events, we selected the period 1980–2016, for which most of the hydro-meteorological and biophysical data are available. Adopting a uniform analysis period, along with the selection of catchments with high-delineation quality, resulted in limited data coverage in regions like Asia, Australia, and northern and central Africa.
Table 1
Spatial resolution and temporal coverage data used in the analysis
Dataset | Temporal resolution | Spatial resolution | Temporal coverage |
In-situ river streamflow data – GSIM62,63 | Daily statistics per month (MAX, MIN, MEAN) | nodes (catchment outlets/ hydrometric stations) | Varying (1900–2016) |
Precipitation – MSWEP64 | monthly | 11 km | 1979–2022 |
Soil Moisture - GLEAM (Global Land Evaporation Amsterdam Model) v3.6a65 | monthly | 25 km | 1980–2020 |
Surface water extent (Global Surface Water)69 | monthly | polygons/ centroids | 1984–2020 |
No noise (smoothed) Normalized Difference Vegetation Index – NOAA67 | 7-day composite | 4km | 1981–2022 |
Total Water Storage (TWS) anomaly is computed as standardized deviation of the GRACE satellite Liquid Water Equivalent (GRACE)66 | monthly | 50 km | 2002–2020 |
Drought and flood definition and propagation
We identified extreme high and low anomalies in the hydro-meteorological and biophysical time series using a threshold level method70,71 (TLM). This approach allows us to quantify the severity and duration of floods and droughts without fitting a distribution to the data, because fitting a distribution is a process that can bias the estimated anomalies72–74. Further, the TLM approach allows to compare between hydrometeorological variables, which is required when studying drought propagation70,75.
For the detection of high streamflow anomalies, we applied the fixed threshold level method to the monthly maximum daily streamflow values. Using maximum daily flows allows capturing within-month streamflow variability and thus the highest flows that could lead to flooding. For the detection of low anomalies, we applied the variable threshold level method to all hydrometeorological and biophysical variables analysed, except for streamflow. Threshold analysis can be difficult to apply to streamflow time series in arid and cold basins due to many zero values. To overcome this limitation, we applied a combined variable TLM and CDPM (consecutive dry period) method to the monthly mean daily streamflow values, following Van Huijgevoort et al. (2012)76. The method allows detection of low anomalies for ephemeral rivers (more common present in arid regions as well as in rivers located in cold areas, where streamflow equals zero when the temperature drops below zero for a longer period76).
The combined variable TLM and CDPM method is only applied to stations whose time series, grouped by month, have more than 5 percent of monthly mean daily streamflow data equal to zero. In this case, the consecutive number of dry months is calculated by considering the zero values in the streamflow timeseries. For non-zero flows, we compared the threshold to the respective streamflow percentile calculated using only the non-zero flows. The consecutive number of dry month time series is then transformed into percentile values and rescaled according to the fraction of zero values. Percentiles of positive flow values are then combined with the rescaled percentile of zero flow. An overview of the methods used to detect high and low anomalies for each variable is provided in Table 2.
Finally, the threshold method is sensitive to the subjective selection of the threshold. In our analysis, we used an exceedance threshold value of 0.85 for detecting extremely low conditions in monthly mean daily streamflow and other hydro-meteorological and biophysical variables. A threshold level of 0.15 was used to detect extreme high conditions in monthly maximum daily streamflow. We chose these values after assessing the sensitivity of the results to threshold choice (see section Sensitivity Analysis in Methods).
After computing anomalies for all variables in all catchments, we detected drought and flood events. We defined drought (D) as a negative anomaly below its threshold in at least one hydrological variable (i.e., soil moisture surface, soil moisture root, monthly mean daily streamflow, total water storage). Other hydro-meteorological and biophysical variables such as precipitation, surface water extent, and NDVI were used to indicate the severity of the drought and its propagation through the hydrological cycle. We defined flood (F) events as positive anomalies of monthly maximum daily streamflow above its threshold.
For each drought event, we estimated the start and end date based on the start date of the first negative anomaly detected in any variable and the end date of the last anomaly accounting for the propagation of the drought in the hydrological cycle. Consequently, drought duration was calculated as the number of consecutive months in which low anomalies and their propagation in the hydrological system occur. The severity of the drought was calculated as the maximum deficit (in percentile), detected during the duration of the drought for each of the variables under negative anomaly. For each flood event, we determined its start and end date based only on the occurrence of the streamflow peak. Accordingly, flood duration corresponds to the duration of the streamflow peak (in months) and flood severity to its surplus in relation to the threshold (in percentile). To investigate whether and how flood timing shift varies by event magnitude, we categorize the detected drought events into “moderate” and “severe” events. Moderate events as are defined as anomalies between the 15th and the 5th percentile and severe events as above the 5th percentile.
Consecutive and Compound drought-flood detection
We defined compound drought-flood (D&F) when drought conditions occur at the same time as high streamflow, and consecutive drought-to-flood (DtoF) when drought conditions are immediately followed by high streamflow in the month following the end of the drought. For each station, both D&F and DtoF were compared with all detected floods to determine the percentage of floods compound or occurring shortly after drought conditions (i.e., ratio between D&F or DtoF events and the total F events, Fig. 1). Since the complex visualization of the results due to the high concentration of stations in certain regions, we computed the average ratio over a 10x10 km grid. Results were also aggregated by climate region according to the categorization provided by GSIM for each catchment, which was based on the Köppen–Geiger climate classification system77.
Seasonality analysis
The mean date of occurrence of detected drought and flood events was computed through circular statistics34,37,78–80. First, we converted the start month \({M}_{i}\)of occurrence of drought (flood) event \(i\) at station \(j\) into an angular value:
\({\theta }_{i,j}= {M}_{i,j} \frac{2\pi }{12}\) ; \(0{\le \theta }_{i,j}\le 2\pi\) and 12 the number of months in a year (1)
Then, we computed the average date of occurrence \({\stackrel{-}{M}}_{j}\) of drought (flood) at a station \(j\) as:
$${\stackrel{-}{M}}_{j}=\left\{\begin{array}{c}{\text{t}\text{a}\text{n}}^{-1}\left(\frac{{\stackrel{-}{y}}_{j}}{{\stackrel{-}{x}}_{j}}\right) \bullet \frac{12}{2\pi } {\stackrel{-}{x}}_{j}>0, {\stackrel{-}{y}}_{j}\ge 0\\ \left({\text{t}\text{a}\text{n}}^{-1}\left(\frac{{\stackrel{-}{y}}_{j}}{{\stackrel{-}{x}}_{j}}\right)+ \pi \right) \bullet \frac{12}{2\pi } {\stackrel{-}{x}}_{j}\le 0 \\ \left({\text{t}\text{a}\text{n}}^{-1}\left(\frac{{\stackrel{-}{y}}_{j}}{{\stackrel{-}{x}}_{j}}\right)+ 2\pi \right) \bullet \frac{12}{2\pi } {\stackrel{-}{x}}_{j}>0, {\stackrel{-}{y}}_{j}<0\end{array}\right.$$
2
With
$${\stackrel{-}{x}}_{j}= \frac{1}{n} \sum _{i=1}^{n}\text{cos}\left({\theta }_{i}\right)$$
3
$${\stackrel{-}{y}}_{j}= \frac{1}{n} \sum _{i=1}^{n}\text{sin}\left({\theta }_{i}\right)$$
4
where \({\stackrel{-}{x}}_{j}\) and \({\stackrel{-}{y}}_{j}\) are the cosine and sine components of the average date, and \(n\) is the number of years of data used at that station. These components were used to calculate the concentration\({R}_{j}\)
$${R}_{j}= \sqrt{ {{\stackrel{-}{x}}_{j}}^{2}+ {{\stackrel{-}{y}}_{j}}^{2}}$$
5
which expresses the consistency in timing around the average date of occurrence. \({R}_{j}\) values equal to one indicate that the data are strongly directional and hence, all drought (flood) events in the analysed station occur in the same month of the year. \({R}_{j}\) values equal to zero indicate that the data are not seasonal, and hence drought (flood) timing is widely dispersed throughout the year.
The approach allowed us to show for each station globally the average start month of the flood for consecutive (DtoF) and compound (D&F) drought-flood events (Fig. 2a and Fig. 3a). We then computed the percentage of occurrence of drought (flood) events for each month of the year with respect to all the months where we identified a drought (flood) in the analysed station. The results were aggregated according to Köppen–Geiger climate classification and latitude values to consider catchments with similar precipitation and temperature patterns throughout the year (Fig. 2c and Fig. 3b).
Flood timing shifts during drought-flood events
We analysed the shifts in the timing of flood events (i.e., flooding occurring earlier or later in the year than expected), for single floods and for floods that preceded or compounded with drought conditions (i.e., drought-flood events). These shifts were computed as the circular difference between the date of occurrence of individual flood events for a particular catchment and the mean date of occurrence of all single floods of that catchment.
To investigate the difference between the timing of drought-flood events and that of single floods, we quantify their difference in two ways. First, we calculate the absolute difference in the mean timing of the drought-flood events and the single floods. This yields a shift in flood timing expressed in months, which indicates how much drought-flood events are shifted compared to single floods (Figs. 5 and 6). To quantify how anomalous these timings are, we also rank ordered all flood timings per catchment from the smallest to the largest (absolute) difference between the flood timing of the individual event and the mean date of flooding of all single flood events. We subsequently calculate the mean rank order of the drought-flood events per catchment, and derive a non-exceedance probability of size of the timing difference between drought-flood shift compared to single floods. In the case of one single drought-flood event in a catchment, this non-exceedance probability equates to the rank order of the single event. However, since averaging out values would cause the percentile mean to converge to the 50th percentile, we employed Monte Carlo simulation to quantify how anomalous the timing of drought-flood events is by comparing the likelihood of non-exceedance of their rank order in their timing. This simulation involved generating 20,000 random samples of n uniformly distributed numbers between 0 and 100 to estimate their mean percentile, giving insight on the frequency at which a particular mean occurs in relation to all other calculated means (Fig. 4).
Further, we tested the sensitivity of our results (earlier and later floods) to the flood timing concentration, accounting only for catchments with unimodal (one distinct peak over a year) flood patterns (Supplementary Fig. 15). The analysis yielded similar spatial distribution of earlier and later flood timing in both consecutive and compounded drought-flood events, compared to the distribution obtained by also including catchments with bimodal or uniform flood patterns (Figs. 5 and 6).
Sensitivity analysis according to different drought thresholds
Since the results of the threshold method can be sensitive to the choice of the threshold level, we applied different threshold values to test the robustness of our results. Specifically, the 15th, 10th, 5th and 1st percentiles were used to detect drought events, and the 85th, 90th, 95th and 99th percentiles to identify flood events. The probability density function plots of flood timing shift for DtoF and D&F events (Supplementary Fig. 16) reveal a consistent pattern across various drought and flood thresholds.
In DtoF events, the kernel density estimates of the flood timing shift exhibit right-skewed distributions, suggesting a prevalence of later floods. Curves related to drought thresholds between the 15th and 5th percentile present a similar concentration of data points toward later floods. The same applies for curves related to flood threshold between the 85th and 95th percentile. The only exceptions are the curves representing drought anomalies below the 1st percentile and flood anomalies above the 99th percentile, which exhibit a slightly flatter shape with more pronounced left tails. This pattern may be attributed to a lower number of detected events compared to those identified using a higher and lower threshold respectively for drought and flood detection.
The analysis of the probability density curves for flood timing shift in D&F events shows consistently a concentration of values towards later floods (peaks located at the positive x-values), for different drought and flood thresholds. The curve representing flood anomalies above the 99th percentile displays a slightly flatter shape compared to the other curves. The distributions of the curves exhibit a distinct right-skew pattern as drought severity increase, except for the curve representing drought anomalies below the 1st percentile. This curve exhibits a higher concentration of values towards earlier floods compared to the other curves.