2.1 Observed and simulated streamflow datasets

The GSIM archive is used as daily observational discharge for this analysis. Daily streamflow simulations available through the ISIMIP are used, with historical simulations (forced with observational climate in ISIMIP2a and bias-corrected climate model outputs in ISIMIP2b) spanning from 1971 to 2005 (Gosling et al., 2019) and future simulations (ISIMIP2b) covering the 2006–2099 period (Frieler et al., 2017). Six GHMs are considered: H08 (Hanasaki et al., 2008a, b), LPJmL (Schaphoff et al., 2013), MPI-HM (Stacke and Hagemann, 2012), ORCHIDEE (Guimberteau et al., 2014, 2018), PCR-GLOBWB (Wada et al., 2014; Sutanudjaja et al., 2018) and WaterGAP2 (Müller Schmied et al., 2014, 2016). These models were selected as they have provided discharge data within phases 2a and 2b of ISIMIP at the time this study began (June 2018). A summary of the similarities and differences across participating GHMs is provided in Sect. 1.2 in the Supplement.

To assess the model structural uncertainty across GHMs, trends in streamflow extremes simulated under observational atmospheric forcing, available through the Global Soil Wetness Project Phase 3 (GSWP3) reanalysis (Kim, 2017), were compared to observed trends. The influence of the high uncertainty in climate models (Kumar et al., 2013; Kiktev et al., 2003) on streamflow simulations was assessed by comparing observed trends and trends simulated when using atmospheric forcing from four GCMs for the historical period (“hindcast” simulations; hereafter referred to GCMHIND atmospheric forcing). These GCMs were bias-corrected but their simulations have different sub-monthly, inter-annual and decadal variability, and thus the hindcast simulations reflect both GHM and GCM uncertainty. To quantify the implication of model uncertainty for future projections of flood hazard, trends simulated under projected climate change by the end of this century (using the same four GCMs) were also assessed for two greenhouse gas concentration scenarios, RCP2.6 (hereafter referred to GCMRCP2.6 atmospheric forcing) and RCP6.0 (hereafter referred to GCMRCP6.0 atmospheric forcing). As a result, four simulation settings were used in this study, denoted by the atmospheric forcing; an overview is given in Table 1. These settings comprise two historical runs (GSWP3 and GCMHIND runs) and two future runs (GCMRCP2.6 and GCMRCP6.0), collectively amounting to a total of 69 simulations (see Table S3 with full list of simulations).

Table 1Summary of streamflow observation and simulation datasets used in this study. GSIM was used as the observed streamflow database. Streamflow simulations were obtained from six GHMs (H08, LJPmL, MPI-HM, ORCHIDEE, PCR-GLOBWB and WaterGAP2). One observational atmospheric forcing dataset (GSWP3) and outputs of four GCMs were used as input for streamflow simulations.

Download Print Version | Download XLSX

For GSWP3 simulations, a preliminary analysis (see Sect. 4 in the Supplement) shows that both “naturalised runs” (i.e. human water management not taken into account) and “human impact runs” (i.e. human water management inputs were used) exhibit similar characteristic of trends in MAX7 index. Some potential reasons for negligible impacts of human water management are the spatial distribution of stream gauges (may be biased toward regions with insignificant changes in water management during the 1971–2005 period), or the inclusion of small catchments (more that 3000 catchments with reported area less than 9000 km²); thus, floods are more sensitive to changes in climate forcing relative to the accumulated basin-wide influence of human impacts. Naturalised runs were therefore chosen, since this setting is available for more GHMs (six) when compared to the human impact setting (four). Although significant efforts were made by ISIMIP to keep the setting across simulations as consistent as possible, there were some differences in model versions and input data (e.g. WaterGAP2.2 (ISIMIP2a) was used in ISIMIP2a while WaterGAP2.2c was used in ISIMIP2b; ORCHIDEE (Guimberteau et al., 2014) was used in ISIMIP2a while ORCHIDEE-MICT (Guimberteau et al., 2018), with improvements on high latitude processes, was used in ISIMIP2b). Although the influence of versioning is minor for WaterGAP2, the potential effects of technical discrepancies cannot be checked in the context of this study, as not all required simulations are readily available (see our discussion in Sect. 3.3 in the Supplement). In addition, owing to technical requirements across GHMs, different models do not have the same set of coastal cells, which may lead to some minor effect to the statistics when averaged across all simulation grid cells.

2.2 Catchment selection and simulated streamflow extraction for observation–model comparison

To enable an observation–model comparison, simulated discharge needs to be extracted from gridded model output. Large-scale hydrological models, however, generally do not simulate discharge accurately over small-to-medium size catchments due to the coarse resolution of river network datasets in their routing schemes (Hunger and Döll, 2008). To address this limitation, previous GHM evaluations usually selected large catchments (a threshold of 9000 km² was adopted, approximating the size of a 1^∘ longitude–latitude grid cell), and routed discharge (unit: m³ s⁻¹) at the outlet of the catchment was used as simulated streamflow for a specific catchment (Zhao et al., 2017; Veldkamp et al., 2018; Zaherpour et al., 2018, 2019; Liu et al., 2017). For evaluation studies that used relatively small catchments (e.g. area less than 9000 km²), the un-routed runoff simulation (unit: mm d⁻¹) was extracted while observed discharge was converted to runoff using catchment area prior to comparison (Gudmundsson et al., 2012b; Beck et al., 2017). To increase the sample size for the model–observation comparison (the first objective), the present study used both (i) daily un-routed runoff for small catchments and (ii) daily routed discharge simulations for large ones, and thus two extraction procedures were adopted. A summary of these extraction procedures is provided below while detailed technical descriptions are provided in Sect. 2 in the Supplement.

• For catchments with an area from 0 to 9000 km²: un-routed runoff (mm d⁻¹) was extracted and then converted into discharge (m³ s⁻¹) by multiplying averaged runoff with catchment area reported in the station metadata. Specifically, catchment boundaries were superimposed on the GHM grid to obtain the weighted-area tables, which were then used to derive averaged runoff from the un-routed runoff simulation. To avoid double-counting runoff from the same grid points, runoff for catchments that share similar weighted-area tables (i.e. similar simulated streamflow would be extracted – see Sect. 2 in the Supplement for a detailed description) was averaged (using catchment areas as weights) and a single “averaged time series” was used in place of the runoff from the component catchments.

• For catchments with an area greater than 9000 km²: the “discharge output” approach (Zhao et al., 2017) was adopted to extract routed discharge (m³ s⁻¹) from the GHM cell corresponding to the outlet of each catchment.

To ensure sufficient data are available for historical trend analysis, only GSIM stations with at least 30 years of data available during the 1971–2005 period were considered (each year having at least 335 d of available records, implying that annual maximum of a specific year is identified only when more than 90 % of the daily record is available). These relatively strict selection criteria also enable a comparison between this study and preceding observation-based investigations (Gudmundsson et al., 2019; Hodgkins et al., 2017). As catchment boundary shapefiles (Do et al., 2018a) were used to extract simulated streamflow for small catchments, stations were further filtered using two criteria: (i) availability of reported catchment area and (ii) catchment boundary being accompanied by a “high” or “medium” quality flag (i.e. the discrepancy between reported and estimated catchment area is less than 10 %).

Figure 2Locations of 3666 streamflow observations (blue dots: 3024 non-averaged time series; yellow dots: 624 averaged time series, where geographical coordinates were averaged from all component gauging coordinates) selected from GSIM archive for the model–observation comparison. Grey dots indicate GSIM time series that were removed due to insufficient data availability or quality.

A total of 4595 stations satisfied the quality selection criteria, of which large catchments (i.e. area greater than 9000 km²) where no suitable grid cell could be identified were further removed (11 catchments). For cases of two or more small catchments (i.e. area less than or equal to 9000 km²) with similar weighted-area tables, the “averaged time series” (using catchment areas as weights) was calculated. A total number of 1542 time series fell in this category and were aggregated into 624 “averaged time series”. Figure 2 shows the spatial distribution of the final dataset for model–observation comparison, containing data for 3666 locations (3042 non-averaged time series and 624 averaged time series). The majority of available catchments are located in North America and Europe, with some regions over Asia, Oceania and South America also covered.

2.3 Detecting trends in annual maximum streamflow

For each streamflow dataset, daily discharge was smoothed to 7 d averages to reduce variability in simulated streamflow, which can arise from the coarse routing parameters of GHMs (Dankers et al., 2014). The annual maximum time series of 7 d averaged discharge (labelled as the MAX7 index in the GSIM archive) was then derived to represent peak flow events. For gridded datasets, the “centre averaged approach” (e.g. averaged streamflow of 7 January is the mean value of 4–10 January) was used (the common setting of the CDO software, freely available at https://code.mpimet.mpg.de/projects/cdo, last access: 1 March 2020), and the MAX7 time series was therefore derived for each GSIM station using this same approach. As a result, the derived value of the MAX7 index is slightly different to the value available in the online version of GSIM (Gudmundsson et al., 2018b), which applied a “backward-moving average” technique (e.g. averaged streamflow of 7 January is the mean value of 1–7 January). Our preliminary analysis (not shown), however, indicated that this difference did not lead to substantial changes in the key findings (i.e. similar spatial composition between increasing and decreasing trends).

The magnitude of trends in the MAX7 index at a specific catchment or grid cell was quantified using the normalised Theil–Sen slope (Gudmundsson et al., 2019; Stahl et al., 2010), and the results are expressed in percentage change per decade. The significance of the local trend was assessed using a Mann–Kendall test at the 10 % two-sided significance level (Wilks, 2011). The null hypothesis (no trend) is rejected if the two-sided p value of the test statistic (Kendall's τ) is lower than 0.1, while the direction of the trend (i.e. increasing or decreasing) was determined using the sign of τ.

2.4 Statistical techniques

To explore GHMs' capacity to simulate observed trends and the implication of model uncertainty for projected trends, trends in streamflow extremes obtained from GSIM (observed trends) and ISIMIP simulations (simulated trends) are analysed. The observed trends were available for 3666 observation locations. Simulated trends were available for all 59 033 GHM grid cells (estimated from routed discharge of each grid cell; Antarctica and Greenland were removed). To enable a model–observation comparison, we also extract a subset of simulated trends over the 3666 observation locations (described in Sect. 2.2).

2.4.1 A hypothesis-test approach for comparison of trend characteristics

A range of hypothesis tests (summarised in Table 2; GSWP3 simulations were used to assess GHM uncertainty while GCMHIND simulations were used to assess the combined GCM–GHM uncertainty) was applied to address the first two objectives, which require comparing trend characteristics exhibited from different streamflow datasets. Four characteristics of trends were assessed.

• Trend mean: the mean (percentage change per decade) of trends in streamflow extremes across all gauge- or cell-based time series over a spatial domain. A hypothesis test was adopted to assess whether the trend means exhibited from two specific streamflow datasets (e.g. model vs. observed) are significantly different from each other.

• Trend standard deviation: the standard deviation (percentage change per decade) of trends in streamflow extremes across all gauge- or cell-based time series over a spatial domain. A hypothesis test was adopted to assess whether the trend standard deviations exhibited from two specific streamflow datasets are significantly different from each other.

• Percentage of significant trends (%): the percentage of trends in a domain that are statistically significant, with gauge- or cell-based significance calculated using the Mann–Kendall test at the 10 % significance level. To assess whether the percentage of significant (increasing or decreasing) trends exhibited from a specific streamflow dataset is produced by random chance, a field significance test (Do et al., 2017) was adopted (described in Table 2).

• Trend spatial pattern: the spatial distribution of trends in streamflow extremes over a spatial domain. Pearson's correlation (r statistic) (Galton, 1886; Kiktev et al., 2003) between trends of MAX7 index obtained from two datasets was used as a measure of similarity in the trend spatial structure. The hypothesis test (pattern similarity test) was adopted to assess whether (i) the correlation between simulated trends introduced by GHMs and observed trends is significantly higher than zero, and (ii) the correlation between trends simulated under hindcast atmospheric forcing and observed trends is significantly lower than that between trends simulated under observational atmospheric forcing and observed trends.

Table 2Hypothesis tests conducted to address the first two objectives.

Download Print Version | Download XLSX

2.4.2 Estimating uncertainty of trend characteristics across ensemble members

The third and final objective, which focused on the implications of GCM–GHM uncertainty on projected changes in flood hazard, was addressed by quantifying the spread of trend characteristics (i.e. trend mean, trend standard deviation and percentage of significant trends) exhibited from routed discharge projections under two representative concentration pathways.

The spatial uncertainty of projected trends (GCMRCP2.6 and GCMRCP6.0) was also quantified by calculating intra- and inter-model correlation of the trend patterns across all ensemble members available under the two projections. Intra-model correlation represents spatial uncertainty introduced by the GCM and was calculated from simulated trends introduced by the same GHM (using different simulated atmospheric forcing). Inter-model correlation represents the combined GCM–GHM spatial uncertainty and was calculated for each pair of simulated trends that were (i) introduced by the different GHMs and (ii) forced with different projected atmospheric forcing.

To assess the robustness of GHMs in projecting changes in flood hazard, each grid cell available in the discharge simulation grid was then categorised into one of the five “flood-risk” (here “flood-risk” level is defined as the number of ensemble members projecting significant increasing trends) groups based on the number of GCMRCP2.6 and GCMRCP6.0 simulation members projecting a significant increasing trend (Group 1: no members, Group 2: from 1 to 5 members, Group 3: from 6 to 10 members, Group 4: from 11 to 15 members and Group 5: from 16 to 18 members).

Finally, to assess whether locations projected with an increasing trend by the majority simulations are adequately monitored, each GSIM gauge was sorted into one of these five groups based on the gauge's geographical coordinates. The allocation of gauges to these groups was then analysed to determine whether the most comprehensive global database of daily streamflow records to-date was evenly distributed across the five “flood risk regions”. An inadequate coverage of stream-gauge networks over high-risk regions indicates potentially high vulnerability to future changes in flood hazards, as insufficient data are available to inform decision makers.

Figure 3Normalised Theil–Sen slope for historical trends in flood magnitude (MAX7 index) exhibited over 3666 locations across three streamflow datasets (a: GSIM; b: GSWP3; c: GCMHIND). Multi-model average is shown for simulated trends. Trend is expressed in percentage change per decade. Scatter plots between trends obtained from GSIM and GSWP3/GCMHIND simulated streamflow are provided in (d) and (e).

3.1 Capacity of GHMs to reproduce observed trends in flood hazards

Visual inspection of the normalised Theil–Sen slope across the GSIM time series (Fig. 3a; regional maps provided in Fig. S4) shows a spatial pattern that is consistent with recent findings on trends in observed flood magnitude (Mangini et al., 2018; Do et al., 2017; Mallakpour and Villarini, 2015; Gudmundsson et al., 2019; Burn and Whitfield, 2018; Ishak et al., 2013). Specifically, decreasing trends tend to dominate Asia (most stations located in Japan and India), Australia, the Mediterranean, the western and north-eastern US, and northern Brazil, while increasing trends appear mostly over central North America, southern Brazil and the northern part of western Europe (including the UK). Note that the observation locations are not evenly distributed (86 % in North America and Europe), and thus the confidence of this assessment varies substantially across continents.

Table 3Characteristics of trends in the MAX7 index over the 1971–2005 period across 3666 locations for GSIM observed trends and GSWP3 simulated trends (six GHMs available). Trend mean and trend standard deviation are expressed in percentage change per decade. Correlation was obtained from GSIM observed trends and GSWP3 simulated trends for each GHM. Boldface texts represent values that reject the null hypotheses outlined in Table 2 (hypothesis 1 to 4).

Download Print Version | Download XLSX

The multi-model average of GSWP3 simulated trends (trends simulated under observational atmospheric forcing; Fig. 3b and d) has generally good capacity to reproduce spatial patterns of observed trends. The multi-model average of GCMHIND simulated trends (trends simulated under hindcast atmospheric forcing; Fig. 3c and e), however, could not reproduce some spatial agglomerations of trends in streamflow maxima (e.g. the decreasing trends in south-eastern Australia, increasing trends over north-eastern Europe). This feature indicates the inconsistent climate variability between GCMs and the real world, suggesting GCM climate forcing cannot account for observed trends at sub-continental scale. In addition, GCMs uncertainty can potentially contribute to this inconsistency. Interestingly, the multi-model average of both GSWP3 and GCMHIND simulations generally exhibits a lower magnitude of changes (i.e. closer to “zero change”) compared to the observed trends. This feature is more prominent in GCMHIND (21 simulations available) compared to GSWP3 (six simulations available) and can be explained by two possibilities. The first possible explanation is the nature of averaging, which tends to smooth out variability in trend magnitude across ensemble members, leading to a relatively “close to zero” change across the globe (given that each GCM has stochastic decadal climate variability, so that averaging results forced by GCMs tends to cancel trends). An alternative explanation is that individual simulations also exhibit a lower magnitude of change relative to observation. As Fig. 3 is not sufficient to evaluate the latter possibility, a more detailed comparative analysis between observed trends and individual simulated trends using both historical climate forcings (via GSWP3) and GCM hindcasts was conducted. Specifically, four characteristics of trends in extreme flows (i.e. trend mean, trend standard deviation, percentage of significant trends and trend spatial structure) were assessed for individual simulations and the results are reported in following sections. At the global scale, GSIM observed trends exhibit a mean and standard deviation of −2.4 % and 9.9 % change per decade over the 1971–2005 historical period. Furthermore, there are 7.5 % (12.1 %) stations showing significant increasing (decreasing) trends (detected by the Mann–Kendall test at the 10 % significance level). These numbers, however, are not statistically significant at the global scale.

Table 3 shows the results of the global model–observation comparison using GSWP3 simulated trends across the six GHMs. Compared to observed trends, most simulated trends have a significantly higher global trend mean at the observed locations and lower trend standard deviation. The percentage of locations showing significant trends varies substantially across simulations, but the values were not statistically significant. All GHMs demonstrate low-to-moderate capacity in simulating the spatial pattern of trends (spatial correlation coefficients range from 0.35 to 0.50, indicating that GSWP3 simulated trends account for between 12 % and 25 % of the cross-location variability in the observed trend signal). There is, however, a notable difference in terms of the overall sign of trends simulated by each GHM. This feature indicates that using different GHMs can lead to different interpretations about the overall change in flood hazard at the global scale, despite having a common boundary forcing. Therefore, the “closer to zero” trends of ensemble averages (illustrated in Fig. 3) likely reflect the implication of averaging rather than a systematic bias of GHMs toward a low magnitude of change. As an implication, ensemble averages, though useful, should not be used as the sole reason to infer changes in floods, as it may undermine the actual magnitude of simulated trends. As a result, the following analyses will report the full range (and mean) of each trend characteristic estimated across all ensemble members to communicate the uncertainty underlying the results.

Table 4Characteristics of trends in the MAX7 index over the 1971–2005 period across 3666 locations for GCMHIND simulated trends. Trend mean and trend standard deviation are expressed in percentage change per decade. Intra-model averages of trend characteristics are shown for each GHM. Values in the parentheses show the number of simulations rejecting the null hypothesis (from 1 to 4) outlined in Table 2 (out of four GCMs). Multi-model minimum, maximum, and average values together with those exhibited from GSIM are also provided.

Download Print Version | Download XLSX

Table 4 provides the results of the model–observation comparison using GCMHIND simulated trends (intra-model averages are shown while results of individual simulations are reported in Sect. 4 in the Supplement). Similar to GSWP3 trends, intra-model averages (i.e. calculated from simulations of one GHM) of GCMHIND trends tend to have a higher global mean and lower trend standard deviation than observed. The composition between the percentages of locations showing significant trends varies substantially across simulations and statistical significance was found only for decreasing trends in 3 out of 21 simulations (2 LPJmL simulations and 1 MPI-HM simulation). The multi-model ranges encapsulate the observed trend mean and percentage of significant trends, while the observed trend standard deviation is clearly above the range exhibited from all GCMHIND simulations. The significantly lower simulated trend standard deviation can be partially attributable to the coarse resolution of GHMs' atmospheric and land surface inputs, which may not sufficiently reflect the variation of hydrological processes across small-to-medium catchments.

Among 21 GCMHIND simulations, the “zero similarity” hypothesis (hypothesis 5) was rejected over 13 simulations, indicating that GCM–GHM ensemble members possess some capacity to simulate the spatial structure of observed trends in streamflow extremes. The correlation between GCMHIND simulated trends and GSIM observed trends, however, is significantly lower than that exhibited from GSWP3 simulated trends across all GHMs (reported at Table 3). The results of the similarity assessment are illustrated for a single GHM (H08, as the results were similar for other GHMs) in Fig. 4, where the correlation between observed trends and GSWP3 simulated trends is significantly different from zero. In contrast, the correlation between observed trends and each of the simulated trends under hindcast atmospheric forcing (GCMHIND simulations) is much lower, with two of the four not being statistically higher than zero. These results confirm the substantial influence of atmospheric forcing on the simulated trend pattern relative to the GHM's structure.

Figure 4Model–observation correlation between observed trends and simulated trends across all simulations (GSWP3 and four GCMHIND simulations) of a single model (H08; similar results for other GHMs). Coloured dots indicate actual correlation between a specific simulated trend pattern and observed trend pattern across 3666 locations. Colour lines represent the PDFs of correlation between simulated trend pattern and observed trend pattern obtained through a bootstrap resampling procedure (B=2000).

Download

Table 5Characteristics of trends exhibited from the GSIM/GSWP3/GCMHIND streamflow dataset at the continental scale (each observation location of 3666 sites was sorted into one of the six continents). For simulated trends, only the multi-model average is shown for each region. Trend mean and trend standard deviation are expressed in percentage change per decade. Values in the parentheses show the number of simulations rejecting the null hypothesis described in Table 2 (up to 6 for GSWP3 simulations and 21 for GCMHIND simulations). For GSIM, field significance of increasing and decreasing trends was highlighted by boldface texts. For GSIM, field significance of increasing and decreasing trends was highlighted by boldface text.

Download Print Version | Download XLSX

To further quantify changes at the regional scale, a model–observation comparison (identical to that at the global scale) was conducted over six continents, and the results are summarised in Table 5 (multi-model averages are shown). The trend mean exhibited from GSIM ranges from −10.7 % (Oceania) to 2.4 % change per decade (Europe), while trend standard deviation ranges from 8.3 % (Europe) to 15.8 % change per decade (Oceania). The percentage of significant increasing (decreasing) trends exhibited from GSIM ranges from 3.2 % to 22.6 % (from 6.3 % to 29.1 %), and the composition of significant trends across the six continents is consistent with a previous investigation (Do et al., 2017). The observed percentage of significant trends is found to be above random chance for Europe (increasing flood magnitude) and Australia (decreasing flood magnitude), and this feature is captured quite well by GSWP3 simulated trends, with at least half of the simulations confirming field significance detected from GSIM. Trend characteristics simulated by GHMs at continental scale confirm some important findings from global-scale assessments, suggesting substantial uncertainty of trends in streamflow extremes introduced by GHMs at the continental scale:

• both GSWP3 and GCMHIND simulations generally exhibit a higher trend mean and lower trend standard deviation compared to the observed trend at the continental scale (see also Sect. 3.1 in the Supplement);

• GCMHIND simulations generally exhibit lower capacity to reproduce trend characteristics relative to GSWP3 simulations due to the combined GCM–GHM uncertainty.

For GSWP3 simulations, the spatial correlation is weakest in Asia, as no simulation rejects the null hypothesis of “zero similarity”, while the spatial correlation is strongest in Oceania (mainly southern Australia; correlation of 0.63). Oceania, however, exhibits the highest model–observation discrepancy in trend mean and trend standard deviation, indicating the capacity of a given GHM in terms of the trend spatial structure is not necessarily consistent with its performance in terms of the mean and spread of trends.

GCMHIND trends also suggest the opposite composition between percentages of significant trends compared to GSIM trends (e.g. simulated trends suggest more locations showing significant increasing trends while observed trends suggest the opposite). Among six continents, GCMHIND trends exhibited the lowest correlation (−0.14) in Oceania, whereas GSWP3 suggested the strongest correlation in this continent. This assessment further indicates the substantial impact of atmospheric forcing relative to GHM model structure on the simulated trends in high flow events. It is informative to note that this result is expected, as GCMs (despite having been bias-corrected) generally have low capacity in reproducing the timing of wet or dry periods or the spatial distribution of climate extremes (Kiktev et al., 2007), and GHMs are likely to inherit these limitations when using GCMs' outputs as climate forcing data.

3.2 Determining the representativeness of observation locations in the GHM simulations

To assess the representativeness of observation locations in GHM grid cells, trend characteristics obtained from all simulated grid cells were compared to those estimated from the observation locations (3666 sites globally). For GSWP3 simulations, the results suggest a significant difference between trend characteristics from all model grid cells compared to those obtained from the observation locations (Table 6; multi-model averages shown). This feature is consistent at both global and continental scales, including North America and Europe – the continents with the best stream-gauge density. Specifically, the trend mean tends to get closer to zero, while the trend standard deviation obtained from all grid cells tends to be higher than that over observation locations. The difference between the percentages of significant increasing and decreasing trends across all grid cells also gets smaller. For instance, the percentage of observation locations showing significant increasing (decreasing) trends over Oceania is 3.7 % (22.1 %) for GSWP3 multi-model averages (reported in Table 5), while the corresponding values are 10.7 % (15.1 %) when all grid cells are considered (reported in Table 6). Additionally, field significance for increasing (decreasing) trends is detected in two (four) out of six simulations over Oceania, while the same feature could not be detected over the observation locations. These findings confirm that trends exhibited from observation locations are not a representative sample of trends obtained from all simulation grid cells, which has also been suggested through Fig. 2. As a result, a common model–observation picture of changes in global flood hazard remains elusive. To enable a holistic perspective of changes in extreme flows, it is therefore crucial to improve not only models' capacity, but also data accessibility and expand streamflow observational networks to ensure unbiased samples are available for large-scale investigations.

Table 6Characteristics of simulated trends across all grid cells at both continental and global scales (multi-model averages are shown). For each simulation, cell-based trend mean and trend standard deviation were compared to those of gauge-based trends (reported in Table 4). Values in parentheses represent the number of simulations that reject the null hypothesis described in Table 2 (up to 6 simulations for GSWP3 and 21 simulations for GCMHIND). GSIM results are also provided for reference. For GSIM, field significance of increasing and decreasing trends was highlighted by boldface text.

Download Print Version | Download XLSX

The findings using GCMHIND simulations are similar in terms of the trend mean (closer to zero) and trend standard deviation (higher) across all grid cells relative to the observation locations. Across all land areas, the composition of the percentages of land mass showing significant trends exhibited by GCMHIND simulations contradicts that obtained from the GSWP3 simulations for many continents. For example, GSWP3 simulations suggest more land areas showing significant decreasing trends than increasing trends over Asia and Oceania while GCMHIND simulations indicate an overall increasing change in extreme flows over the same continents. This feature further confirms the importance of uncertainty in atmospheric forcing in driving the spatial structure of the simulated trends, which will be explored further in the next section.

3.3 The implication of simulation uncertainty on the projection of trends in flood hazard

This section focuses on the uncertainty in simulated trends under projected climate forcing at the global scale. For MPI-HM (no simulation for HadGEM2-ES forcing), streamflow was only simulated across the main stream network (approximately 45 % of the global land grid cells), and thus three simulations of this GHM were removed from the analysis. As a result, only 18 ensemble members were used to explore the uncertainty in projected trends (GCMRCP2.6 and GCMRCP6.0 – trends estimated for the 2006–2099 period and all cells were considered).

Table 7The uncertainty in the characteristics of projected trends (GCMRCP2.6 and GCMRCP6.0) across 18 members at the global scale (five GHMs). Trend mean and trend standard deviation have unit of %-change per decade. At-site significance of trend was identified using Mann–Kendall test at 10 % level and the percentage of grid cells showing significant increasing and decreasing trends was reported (no field significance test was conducted). Intra-model average value of each metric across is shown for each GHM (numbers of simulations are provided in the first column).

Download Print Version | Download XLSX

Table 7 shows a relatively low spread of the global trend mean (ranging from −1.3 % to 0.8 % change per decade; multi-model average of 0.0 % change per decade for both GCMRCP2.6 and GCMRCP6.0) and trend standard deviation (ranging from 1.8 % to 4.1 % change per decade) across ensemble members. LPJmL and ORCHIDEE generally suggest a decreasing trend at the global scale, evident through the negative global mean and more grid cells showing significant decreasing trends. The standard deviation of trends in future simulations is substantially lower than the historical run (reported in Table 6). This feature is potentially due to the capacity of longer time series in capturing the inter-decadal variability of the streamflow regimes, with both dry and wet periods being considered (Hall et al., 2014). Projected trends under the RCP2.6 scenario generally have lower mean and lower standard deviation closer to zero compared to those introduced by the RCP6.0 scenario, reflecting the nature of an ambitious “low-end warming” scenario, when anthropogenic climate change reaches its peak in the middle of the 21st century followed by a generally stable condition.

Interestingly, although most models suggest relatively moderate changes in the global trend mean, the composition between percentages of grid cells showing significant trends varies substantially, ranging from 7.5 % (7.1 %) to 30.1 % (35.0 %) for significant increasing (decreasing) trends at the 10 % level, with RCP6.0 generally exhibiting higher values. This finding indicates that inferences of changes focusing on global averages may mask significant regional trends, as there was a substantially high percentage of locations exhibiting significant increasing and decreasing trends exhibited in individual models.

Uncertainty in the spatial structure of trends in streamflow extremes is further investigated using both intra-model (to reflect GCM uncertainty) and inter-model correlations (to reflect the combined GCM–GHM uncertainty). A more robust spatial pattern of projected trends under RCP6.0 was found, indicated through generally higher intra- and inter-model correlation compared to those exhibited from trends simulated under RCP2.6 across all GHMs. This feature potentially reflects the less contrasted regional climate change of RCP2.6 relative to RCP6.0. The inter-model correlation is consistently lower than intra-model correlation due to the combined uncertainty of both GHMs and GCMs.

Figure 5Number of simulations showing statistically significant trends at the 10 % level at each grid cell. Panels (a) and (b) show results for the assessment of increasing trends, while (c) and (d) show results for significant decreasing trends. (a, c) Results of GCMRCP2.6 simulations; (b, d) results of GCMRCP6.0 simulations.

To quantify the robustness in terms of regions with significant trends in streamflow extremes, the number of simulations showing significant increasing and decreasing trends was counted for each grid cell (values ranging from 0 to 18). As shown in Fig. 5a and c, the projections under RCP2.6 do not suggest many regions with an increasing trend for most ensemble members, but consistently suggest decreasing trends over the majority of Africa, Australia and western North America. Although both scenarios suggested a similar spatial pattern, projections under the RCP6.0 scenario (Fig. 5b and d) show a substantially higher robustness in terms of regions with significant changes over time in streamflow extremes. For instance, significant increasing trends are projected consistently over southern and south-eastern Asia, eastern Africa, and Siberia, while high agreement of decreasing trends is found over southern Australia, north-eastern Europe, the Mediterranean and north-western North America. These findings share some similarity with a previous investigation that used the ISIMIP Fast Track simulations (published before the ISIMIP2a and 2b simulations used here) to identify regions projected with an increasing magnitude of 30-year return level of river flow (Dankers et al., 2014). Specifically, both studies suggest overall (1) an increasing trend over Siberia and South-East Asia and (2) a decreasing trend over north-eastern Europe and north-western North America. The present study, however, additionally highlights a dominant decreasing trend over Australia, which was not shown previously. The different numbers of ensemble members (45 in Dankers et al., 2014, and 18 in the present study) and greenhouse gas concentration scenario (RCP8.5 in Dankers et al., 2014, and RCP2.6 and RCP6.0 in the present study) between two studies indicate that the choice of GCM–GHM ensemble and greenhouse gas concentration scenarios could lead to substantially different projections of changes in flood hazard at the regional scale.

These results suggest the key role of GCM uncertainty in projections of changes in flood hazards, emphasising the importance of a flexible adaptation strategy at the regional scale that can take this uncertainty into account (Dankers et al., 2014) such as increasing flexibility in reservoir operations, focusing on improved infrastructure resilience and preparing for uncertain changes in flood hazards. Such a strategy is achievable only through a reliable and robust understanding of the change in flood hazards. The assessment of the representativeness of streamflow observations (Sect. 3.2), however, demonstrated that the observation locations selected for this assessment are not a representative sample of the entire land mass. As a result, inference of changes in flood hazard may be biased toward well-observed regions. To further highlight the potential impact of limitations in observed streamflow datasets, the proportion of available stream gauges located in regions with different levels of projected “flood risk” was assessed. We first categorised each simulation grid cell into one of the five “flood-risk” groups. Note that in this analysis, “risk” is defined as the number of simulations projecting a significant increasing trend, rather than the prominent definition of risk as the combination of hazard, exposure and vulnerability (Kron, 2005). In this analysis, the RCP6.0 scenario was chosen as it yielded a higher global “risk” of flood hazard relative to the RCP2.6 scenario.

Figure 6 presents the percentage of all simulated grid cells (a) categorised in each of the five groups, and of GSIM stations located in each group (b). As can be seen, 11.7 % of grid cells fell into the “high-risk” groups (8.9 % from Group 4 with 11–15 ensemble members, and 1.8 % in Group 5 with 16–18 ensemble members), while 68.9 % of grid cells fell into the “low-risk” groups (22.0 % for Group 1 with no ensemble members, and 46.9 % for Group 2 with 1–5 ensemble members). Of all GSIM stations, only 0.9 % are located in high-risk grid cells (no station located in Group 5 grid cells) compared to 89.5 % of stations located in low-risk grid cells (35.4 % for Group 1 and 54.1 % for Group 2). The uneven distribution of stream gauges indicates potential difficulties in using observational records to provide an assessment of global or regional changes in flood hazard, which in part arises from data caveats associated with the spatio-temporal coverage and quality of observed gauge records across the globe. This finding further suggests the urgent demand for ongoing efforts to make streamflow observation more accessible. In addition, new innovations in remote sensing (Gouweleeuw et al., 2018) or development of runoff reanalysis (Ghiggi et al., 2019) should also be supported to complement the understanding of changes in floods for locations that were not observed by stream gauges.

Figure 6Percentage of grid cells (“Landmass”) grouped by the number of simulations projecting a significant increasing trend under the RCP6.0 scenario, and the percentage of streamflow stations (“GSIM”) assigned into each group. The range of possible simulations is from 0 to 18 and binned into five groups (Group 1: no members, Group 2: from 1 to 5 members, Group 3: from 6 to 10 members, Group 4: from 11 to 15 members and Group 5: from 16 to 18 members). To identify which group a specific station belongs to, the geographical coordinates of that station was superimposed on top of the global “flood-risk” map.

Download