Trend estimation.
Throughout this work, we use the Theil-Sen Slope estimator38 to quantify trends and apply the non-parametric Mann-Kendall test39 to evaluate significance levels. In the damage analyses,
trends are relative to the annual mean damage of the reference period 1980–1990 in the corresponding region or
subregion.
Flood modeling and definition of subregions by discharge trends.
We subdivide the nine geographical world regions into subregions with positive and negative trends in annual discharge
maxima over the period 1971–2010. Studies on changes in global discharge patterns are rare and data coverage is
not evenly distributed around the globe. Furthermore, the susceptibility of discharge to human intervention affects
discharge records and complicates disentangling human and climatic forces in observations. We therefore derive discharge
trends from the harmonized multi-model simulations of the 12 global gridded hydrological models (GHM) participating in
phase 2a of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2a40). We here apply the naturalized experiment referred to as‘NOSOC’ in
the ISIMIP2a protocol, meaning that no human impacts, such as dams and water abstractions on river flow were considered.
This is legitimate for two reasons: 1) to ensure consistency with river routing simulations that do not account for
human regulation of rivers, and 2) based on aprevious study for some major basins that showed that the shape of the
hydrograph, for peak daily flow, is not significantly different between natural and human impact experiments41. Furthermore, this allows us to better isolate climate
induced changes in river discharge (SI Sec 3.1).
Here, the 12 GHMs were driven by four separate observational (atmospheric) weather data products for the period
1971–2010 providing daily runoff at 30arcmin resolution. For this ensemble of 46 climate data/GHM combinations
(supplementary Tab. SI1), we follow the methodology applied previously in Willner et al. (2018a,b)42, 43, and
first harmonize the output of the different GHMs with respect to their fluvial network using the fluvial routing model
CaMa-Flood (version3.6.244) yielding daily fluvial
discharge at 15arcmin resolution. Especially for peak discharges, CaMa-Flood agrees better with observed
fluvial discharges than the direct output of the hydrological models45. For the subsequent analysis, we then select the annual maximumdaily discharge
for each grid cell.
For each of the 46 simulations of daily fluvial discharge and each grid cell on 15arcmin resolution, we fit a
generalized extreme value (GEV) distribution to the historical time series ofthe annual maximum discharge using
L-moment estimators of the distribution parameters (for details see SISec 3.2 and discussion in Willner et al.
(2018a)42)allowing for a model bias correction following
the approach by46. It has been shown in several
recent publications that such a hydrological modeling chain is able to reproduce patterns in observed flood
impacts12, 13. In addition to these previous studies, we account for current flood protection
standards at the sub-national scale from the FLOPROS database47. For the final assessment, we re-aggregate the high resolution flood depth data to
a 2.5arcmin resolution by 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 waslarger than zero.
Socio-economic data sources.
We use gridded GrossDomestic Product (GDP) data reported in purchasing power parity (PPP) in 2005 USD from the
ISIMIP project48 with a spatial resolution of 5arcmin from
1971–2010 as a proxy for the distribution of assets. Gridded GDP data was obtained using a downscaling methodology
49 in combination with spatially-explicit population
distributions from the History Database of the Global Environment (HYDEv3.2)50, 51 and national GDP
estimates52. Downscaled GDP data is available in 10y
increments and linearly interpolated across decades. To estimate asset values more precisely we convert gridded GDP
data into gridded capital stock, using annual national data on capital stock (in PPP 2005 USD) and GDP from the PennWorld
Table (version 9.1, https://www.rug.nl/ggdc/productivity/pwt/).
For each country the annual ratio of national GDP and capital stock was calculated and smoothed with a 10-year running mean to generate a
conversion factor, which was then applied to translate exposed GDP into asset values.
Observed asset damages are taken from reported flood damages from the NatCatSERVICE1 database collected by MunichResince 1980, excluding flash flood events or flooding
caused by tropical cyclones. We adjusted all flood damage estimates for inflation to the reference year 2005 using
country-specific consumer price indices(CPI), i.e., expressing them in the same base year as the GDP data. To do so, we
constructed a conversion factor for each country based on all reported damages for a country-specific event in 2005 and
the regularly CPI-adjusted values reported in Munich Re’s NatCatSERVICE database in the base year 2016.
Multiplying CPI-adjusted reported flood damages by this conversion factor results in CPI-adjusted damages for
2005. Event-specific damage estimates were then aggregated to year-country and year-region level in order to
be comparable with simulated river floods for which only the annual maximum was considered. Thereby we assumed that
only
one flood event is observed at each grid cell during a year. To resolve recorded damages in the refined subregional
analysis, we make use of the geocodes given for each general flood event in the NatCatSERVICE dataset and assign the
damage to the risk area of the centroid closest to the geocode location.
Economic damage assessment.
For the estimation of direct asset damages, we apply the regional residential flood depth-damage functions developed
by Huizingaet al. (2017) 14(SI Sec. 3.3).
The quantification of flood damages consists of following steps: 1) determine exposed assets on the grid-level (150
arcmins resolution) based on the flooded fraction obtained from the inundation modelling; 2) determinethe grid level
damage by multiplying the exposed assets by the flood depth and the flood-depth damagefunction; 3) aggregate the
estimated damages spatially to regional/subregional level, and 4) analyze the aggregated damages across different GHM
simulations, assessing model medians and model spread. For steps 1to 3, the open-source probabilistic natural
catastrophe damage framework CLIMADA was used53. To
account for the inhomogeneous but a priori unknown distribution of assets within a grid cell we additionally assume
that
no assets are exposed to a two-year flood event, thus subtracting the two-year flooded fractions from the modeled
flooded fraction before multiplying with the asset value. This is equivalent to assuming that nobody would construct
valuable assets in regions flooded every two years.
The modeled damages for each GHM and grid cell are then aggregated to country andregion level. For comparability
reasons we first aggregate to 9 world regions constructed by grouping countries with geographical proximity and similar
socio-economic structure following the the income group classification of the Worldbank22(Fig. 1b). For regions and
subregions, the median across all GHMs is then compared to reported damages from MunichRe’s NatCatSERVICE
(Fig. 2).
Assessing and accounting for vulnerability.
To include time-varying vulnerability, we apply an approach proposed in previous vulnerability studies12, 13, 15.Comparing modeled and observed damages, a time trend in the
ratio of observed and modeled damages can beobserved that can most likely be explained by changes in socio-economic
vulnerability and/or adaptive capacity. These changes are not properly reflected within the modeling chain and are,
e.g., caused by the fact that the protection standards underlying the FLOPROS database are stationary in time. We apply
an 11-year smoothing on the ratio of reported and modeled damages using Singular Spectrum Analysis (SSA) (Fig.SI1 and
Fig SI2)54 assuming a maximum ratio of observed and
modeled damages of 1000 to remove outliers. Before applying the SSA, missing values were replaced by the median ratio
between 1980 and 2010.
Attributing damages to individual drivers.
Given that the overall trend in damage time series is a superposition of the trends from each individual driver, we
can separate the contributions from each driver by calculating the trend \(\alpha \) of each time
series\({D}_{Cli-1980}\), \({D}_{CliExp}\) and \({D}_{Full}\) and extract climate induced trends in hazard, as well
as trends in exposure and vulnerability, according to:
\({Cli}_{R1980}=\frac{{\alpha }_{Cli-R 1980}}{{n}_{R}}\),
\({Exp}_{R} = \frac{{\alpha }_{CliExp-R} - {\alpha }_{Cli-R1980}}{{n}_{R}}\),
\({Vul}_{R} = \frac{{\alpha }_{Full-R} - {\alpha }_{CliExp-R}}{{n}_{R}}\),
|
(1)
|
We apply a non-parametric trend analysis (Theil-Sen Slope estimator) to estimate
\(\alpha \). Trends are given relative to the annual reported average damages in the time period 1980–1990
(\({n}_{R}\)) in each region or subregion (side panels in Fig. 2).
We additionally provide climate induced trends from time series with 2010
fixed socio-economic conditions (Fig. 3):
\({Cli}_{R2010} = \frac{{\alpha }_{Cli-R2010}}{{n}_{R}}\)
|
(2)
|
Socio-economic trends are assessed from 1980–2010. As climate induced trends
are independent from observational data, we can extend it backward, making use of the full ISIMIP2a time
period and additionally assess trends from 1971–2010.
Climate oscillations and global mean temperature.
In order to avoid interferences with long term temperature increase, we use the pressure based Southern
Oscillation Index as a predictor for ENSO (https://www.ncdc.noaa.gov/teleconnections/enso/enso-tech.php).
Monthly data for AMO and NAO were extracted from the NOAA/Earth System Research Laboratory
(https://www.psl.noaa.gov/data/climateindices/list/)
and PDO from the NOAA/Climate Prediction Center (https://www.psl.noaa.gov/data/climateindices/list/).
We derived annual GMT (daily mean Near-Surface Air Temperature) as the mean of three of the four input
climate forcings provided in ISIMIP2a. We excluded the WATCH dataset because it does not capture the full
historical period.
Contribution of global mean temperature and climate oscillations to climate induced damage trends.
Following the methodology introduced by Najibi & Devineni (2018) and Armal et al. (2018), we apply generalized linear models (GLM) assuming damages to
be log-normally distributed and assuming fixed 1980 socio-economic conditions ()10, 27. In a stepwise procedure we calculated GLMs from all
possible combinations of the predictors ENSO, PDO, NAO, AMO, and GMT and a constant :
\({D}_{Cli} = { \alpha }_{ENSO}\cdot ENSO + {\alpha }_{PDO}\cdot PDO + {\alpha }_{NAO}\cdot NAO +{\alpha }_{GMT}\cdot GMT +{\beta }_{2}\)
|
(3)
|
\({D}_{Cli} = { \alpha }_{ENSO}\cdot ENSO + {\alpha }_{PDO}\cdot PDO + {\alpha }_{NAO}\cdot NAO +{\alpha }_{AMO}\cdot AMO +{\beta}_{1}\)
|
(4)
|
We then select the best model applying a leave-one-out-Cross-Validation (LooCV)
28, which allows to assess model quality outside
the fitting period calculating the “out-of-sample-error” (Supplementary Tab. SI4). The best model is
the one with the smallest “out-of-sample-error”, we additionally test different link functions
(inverse-power, identity, log). To compare the contributions of the different linear predictors across the
different link functions we compare the partial derivatives of the model with regard to the individual
predictors. Finally, we test the residuals for remaining trends applying the non-parametric trend analysis.
Code availability.
All implementations, input and output
data are available upon request and will be made available in a git repository. For damage assessment, we used
the natural catastrophe damage framework CLIMADA available at:
https://github.com/CLIMADA-project/climada_python