Climate signals in river flood damages emerge under sound regional disaggregation

Climate change affects precipitation patterns. Here, we investigate whether its signals are already detectable in reported river flood damages. We develop an empirical model to reconstruct observed damages and quantify the contributions of climate and socio-economic drivers to observed trends. We show that, on the level of nine world regions, trends in damages are dominated by increasing exposure and modulated by changes in vulnerability, while climate-induced trends are comparably small and mostly statistically insignificant, with the exception of South & Sub-Saharan Africa and Eastern Asia. However, when disaggregating the world regions into subregions based on river-basins with homogenous historical discharge trends, climate contributions to damages become statistically significant globally, in Asia and Latin America. In most regions, we find monotonous climate-induced damage trends but more years of observations would be needed to distinguish between the impacts of anthropogenic climate forcing and multidecadal oscillations.

: Overview of hydrological model runs used in this study.
Columns list observational weather data sets used to drive the different global hydrological models (GHMs) (rows): C2010 is the annual trend in damage time series accounting only for variable climate and assuming fixed 2010 socio-economic conditions (D2010) and expresses the annual additional climate contribution to damages relative to the baseline annual damage (annual mean of the period 1980-1995). We derive the climate contribution in 2010 to median damages in 2010 as difference between the start level of the regression line 1980 (1971) and its end level (2010) (Δ 2010 ), *, **, and *** denote significance at 10 %, 5%, and 1% levels, respectively. N in indicates the trend in observed damages (DObs)  normalized to the the baseline annual damage. Grey shadings indicate subregions with high explained variance (R² > 20%):    5 , and the ERA-Interim data (WATCH-WFDEI) 6 . These serve as input for 12 global grid-based hydrological models (GHMs) providing daily runoff at a 30arcmin resolution: CLM4 7 , DBH 8 , H08 9, 10, 11 , JULES-B1 12 (only for GSWP3, PGMFD, WATCH), JULES-W1 13 31 .

Return period assessment and flood inundation mapping
For the subsequent analysis, we select the annual maximum daily discharge for each grid cell. For each simulation for the available historic period (1971-2010) and each grid cell, we fit the generalized extreme value (GEV) distribution 32 to the historical time series of the annual maximum discharge using L-moment estimators 33  For model bias correction, we follow the approach by Hirabayashi et al. (2013) 35 . We map the return period of each event to the corresponding flood depth in a MATSIRO 21 model run driven by observed climate forcing 36 , in bins of 1-year (1 to 100) and 10-year (100 to 1000) return periods (linearly interpolated), providing flood depth at 15arcmin resolution. Results from this observation-driven MATSIRO output have been shown to have realistic consistency in comparison with observation-based data. For this mapping, we further respect a threshold given as current flood protection at the subnational scale. This has recently been compiled in a global database (FLOPROS database) representing the currently best global-scale knowledge in the maximum return period of flood that each country/region can prevent 37 . Here, we use the "Merged layer" of this database, which combines empirical data about existing protection infrastructure ("Design layer"), data about protection standards and requirements set by policy measures ("Policy layer"), and model output from an observed relationship between gross domestic product per capita and flood protection ("Model layer"). This threshold procedure implies that, when the protection level is exceeded, the flood happens as if there was no protection in the first place (for example, dams break); below the threshold no flooding takes place.
For the final assessment, we downscale the resulting flood depth to a 18arcsec resolution using high resolution flood difference maps and re-aggregate to a 2.5arcmin resolution retaining the maximum flood depth as well as the flooded area fraction, defined as the fraction of all underlying high resolution grid cells where the flood depth was larger than zero.

Additional information on damage functions
For the purpose of our analysis we assume that the residential damage function is representative of all other damage categories. This is motivated by the fact that (1) residential damages regularly make up the largest fraction of flood damages, (2) the variation in the damage functions for different categories is small, compared to the uncertainty of the regional distribution of specific asset classes on global level, (3) we are presently unable to differentiate the different categories in the socio-economic data within a grid cell or distribute the different categories across grid cells (except for agriculture) as globally consistent data other than GDP is not available at present.

Comparison of observed and modeled global discharge
Comparing the results of modeled discharge trends (Fig. 1a, Supplementary Fig. 3) to studies on regional and global discharge patterns 38,39,40,41 , we see a good correspondence between observed and modeled trends. In contrast to extreme precipitation, there are globally more areas with declining trends than with increasing trends for the period 1971-2010, this result is in line with empirical studies for the period 1966-2005 39 . For this period, there were more stations showing decreasing trends in North America, Australia and Africa and more stations with increasing trends in Europe, Asia and Latin America. Similarly, the modeled discharge map shows greater shares of decreasing discharge trends in North America, Africa and Australia and rather increasing trends in Europe and Asia, while in Latin America there are also slightly more areas with decreasing trends (Fig. 1a). Considering entire world regions, differences are very likely to occur as gauge stations are very unevenly distributed, while modeled results are provided on a regular 0.25° degree global grid. In Europe, the decreasing tendencies in the Mediterranean found by Blöschl 40 , all studies show the slightly decreasing trends in some areas of Scandinavia. For North America observed patterns are similar to modeled ones, with slightly increasing trends in Central Northern America (CNA) in annual maximum streamflows and rather declining trends in the coastal regions, however increasing trends in CNA are mostly insignificant in both models and observations. In Eastern North America (ENA) a significant declining tendency is visible in observations, while in modeled results we also see declining trends which are mostly not significant.
In Latin America observed maximum discharge increases in the Amazonian region and in South-eastern South America, while significant decreases are observed in North-Eastern Brazil (cf. Fig. 1a with Fig. 2 in 40 ). While a significant decline in maximum stream flows was observed in South Asia 40 , the discharge trend map also shows clearly increasing trends in some areas of India and neighbouring states, however, they are mostly insignificant. In Eastern Asia, observations show significant increases in the time period 1961-2000, which are also present in modeled maximum stream flows, when considering the period 1971-2010. In Southern Australia significant and strong decreases were observed over the same time period, these are especially reflected in South Eastern Australia in the modeled map. In areas where we identify increasing trends there seem to be very few gauge stations available. In their global analysis of peak discharges, Do et al. (2017) do not find strong effects of neither presence of dams nor changes in forest cover 39 . Considering that modeled discharge does not take into account any changes in human river engineering or urbanisation patterns, the overall agreement of signs in the trend analysis also points towards only limited effects of river engineering and urbanisation in this large-scale analysis. The differences in Eastern Europe and India could be caused by human interventions, however, also limitation in hydrological modeling or availability or data quality of observations can explain these differences.

Sensitivity of climate-induced trends in damages to the selection of the baseline exposure
In this work, we aim to quantify the contribution of climate to changes in flood damages in the periods 1980-2010 and 1971 -2010. This exercise is known as "impact attribution". The detection of climate change impacts implies an assessment of the changes in a natural or human system compared to a specified baseline that characterizes its behavior in the absence of climate change. In this definition climate change refers to any long-term trend in climate, irrespective of its cause. To meet the definition of "impact attribution"' provided in the Chapter 18 of the IPCC 42 , the impact of climate change against a non-stationary baseline with temporally varying exposure and vulnerability needs be estimated by comparing observed damages trends against damages in a simulated counterfactual world where exposure and vulnerability vary according to observations but climate is assumed to be stationary.
However, since such an attribution scenario set-up has only recently been added to the ISIMIP3a protocol 43 Table 3). This approach is as close as possible to the definition of attribution in Chapter 18 of the IPCC AR5. As our general modeling approach permits to fix exposure and vulnerability, we can compare climate induced trends as obtained from 1980 fixed exposure with the trends obtained for 2010 fixed exposure. The comparison of the climate-induced trends can be understood as a sensitivity analysis regarding the selection of the baseline scenario. Differences arise when there are assets in 2010 in areas where no assets at all existed in 1980. Assuming 2010-fixed-exposure all damages that were caused on these assets contribute to the climate-induced trend, while changes in damages on these assets in the 1980-fixed-exposure are attributed to the exposure trends. The results of the sensitivity analysis clearly reveal that the derived climate induced trends are robust with regard to the choice of the base year for the exposure.
In follow-up studies based on ISIMIP3a hydrological simulations, the availability of a counterfactual climate scenario will allow for assessments which are independent of the choice of the baseline and allow to clearly distinguish whether impacts on assets are the consequence of shifting climate conditions or of new settlements in flood affected areas.

Limitations and uncertainties
The damage modeling applied in this study shows varying quality across regions (Fig. 2). The NAM region, where the observed damage variability can be almost entirely captured by the modeling chain, highlights the general power of the modeling approach. Resolving the reasons for lower explained variances in other regions will be an important next step to improve the proposed approach as a tool for climate impact attribution. Deficits could be introduced along the entire modeling chain ranging from i) observational climate forcing data whose quality varies by region, ii) translation of daily climate data into discharge, which could arise due to inadequate process-understanding or inadequate representation of direct human influences, iii) translation of discharge into flooded areas, iv) derivation of exposed assets and translation into damages, and v) deficits of reported damages. Here, we discuss sources of uncertainties step-by-step with a focus on potential biases of the climate-induced trends in damages being central for this study: i) Climate reanalysis data constrain the performance of GHMs and may partly be responsible for the spatial variation in explanatory power 44,45 , e.g., due to spatially and temporarily heterogenous global data coverage. However, the good agreement between observed discharge trends and the median discharge trends derived from the model simulations forced by four different observational data indicates that the uncertainties do not lead to large scale systematic biases in trends (Supplementary Notes 1).
ii) There are several sources of uncertainty associated with the translation of the climate forcing into discharge. They may refer to the exact timing and scaling of peak-runoff 44,46 , and the translation into spatially explicit flood-depth on the basis of Digital Elevation Models that are known to be fraught with high levels of uncertainty 47 . Spatial variations in GHM performance may explain part of the spatial variation in the correlations between observed and simulated damages. Thus, validation studies of the GHMs included in ISIMIP2a indicate a variation in performance across hydrobelts 44,46 with a better performance in the wetter equatorial and northern hydrobelts than in drier southern hydrobelts 46 . This is partly reflected in our results, as we observe very low explanatory power in regions such as North Africa and better explanatory power in Asian and North American regions or subregions. In addition GHMs show a tendency to overestimate peak-runoff 44,48 . However, accounting for floodplain storage and backwater effects as implemented in CaMa-flood, that is also applied in this study, leads to a general reduction of peak river discharge and to a better agreement between observed and simulated annual peak discharge 48 . In general the uncertainties of hydrological modeling discussed here are not expected to introduce systematic biases in trends, as modeled and observed discharge trends show very similar patterns (Supplementary Notes 1).
iii) Regarding the exact modeling of flood extents and flooded areas, the river routing model CaMa-flood has been validated in many studies 49, 50, 51 , but not yet within the ISIMIP2a context. Previous studies on selected river basins in Nigeria and Mozambique revealed performance differences ranging from poor to good across river basins and a strong dependence on the return period of input flows 50 . The representation of flooded areas may particularly depend on the adequate representation of flood protection levels only coarsely known and implemented based on FLOPROS in this study. We quantify the sensitivity of damages using two different data versions, the merged-layer (Fig. 2-4) and the modeled-layer of FLOPROS (Supplementary Fig. 7 and 8). Overall, our main findings are robust (explained variances and significance, sign and amplitude of climate-induced trends) with regard to the choice of the FLOPROS layer. Compared to the merged layer, we find slightly higher explained variances in several less developed regions such as LAM and SEA but slightly lower explained variances in NAM when using the modeled-layer of FLOPROS (cf. Fig. 2 to Supplementary Fig. 7). These differences are most likely rooted in an overestimation of policy standards in less developed areas that were used as fill-values for missing data on implemented protection standards in the merged layer. The sensitivity does not address potential biases that could be induced by temporally changes in protection levels not captured by FLOPROS only providing a snapshot. Instead of the explicit modeling through protection levels this effect is captured by the temporal variation in the derived vulnerability factors (Supplementary Fig. 1 and 2). iv) In this study, we observe a tendency of overestimating damages in low-income areas (e.g., SSA) and of underestimating damages in high-income regions (e.g., OCE) (Supplementary Fig. 1 and 2). The major reasons here are the continental depth-damage functions 52 that assume less vulnerability in high-income countries and the protection standards in the FLOPROS database that are much higher in developed countries. However, it is important to note that these only affect the quantification of the vulnerability trends and do not affect the significance, directions, and magnitude of the climate-induced trends.
v) Finally, the reported damages from the NatCatSERVICE database cannot be assumed to be free of biases 53 . However, a likely underreporting in early years will not introduce systematic biases in the estimation of climate-induced trends, but may affect the time-variable vulnerability estimate.
It is expected that it will become possible to explicitly address iii) and v) when extending the analysis to include the last decade. Thus v) could be addressed by checking whether the correlation between observed and simulated damages increases with time indicating that discrepancies found here are reduced by improved reporting of damages. In addition, simulations of flooded areas could be evaluated based on high-quality satellite data becoming available for 1985 onwards.
Our analysis on the relative importance of different teleconnections and GMT may be limited by the fact that climate oscillations are not completely uncorrelated from changes in GMT, i.e., our predictors cannot be assumed to be completely independent. Therefore, we caution in the main text that the selection of GMT as a relevant driver does not necessarily rule out the additional relevance of climate oscillations.