The effect of surge on riverine flood hazard and impact in deltas globally

Current global riverine flood risk studies assume a constant mean sea level boundary. In reality high sea levels can propagate up a river, impede high river discharge, thus leading to elevated water levels. Riverine flood risk in deltas may therefore be underestimated. This paper presents the first global scale assessment of the joint influence of riverine and coastal drivers of flooding in deltas. We show that if storm surge is ignored, flood depths are significantly underestimated for 9.3% of the expected annual population exposed to riverine flooding. The assessment is based on extreme water levels at 3433 river mouth locations as modeled by a state-of-the-art global river routing model, forced with a multi-model runoff ensemble and bounded by dynamic sea level conditions derived from a global tide and surge reanalysis. We first classified the drivers of riverine flooding at each location into four classes: surge-dominant, discharge-dominant, compound-dominant or insignificant. We then developed a model experiment to quantify the effect of surge on flood hazard and impacts. Drivers of riverine flooding are compound-dominant at 19.7% of the locations analyzed, discharge-dominant at 69.2%, and surge-dominant at 7.8%. Compared to locations with either surge- or discharge-dominant flood drivers, locations with compound-dominant flood drivers generally have larger surge extremes and are located in basins with faster discharge response and/or flat topography. Globally, surge exacerbates 1-in-10 years flood levels at 64.0% of the locations analyzed, with a mean increase of 11 cm. While this increase is generally larger at locations with compound- or surge-dominant flood drivers, flood levels also increase at locations with discharge-dominant flood drivers. This study underlines the importance of including dynamic downstream sea level boundaries in (global) riverine flood risk studies.


Introduction
Currently, global flood risk studies either examine riverine or coastal floods (Jongman et al 2012, Hallegatte et al 2013, Hirabayashi et al 2013, Ward et al 2013, Hinkel et al 2014, Vitousek et al 2017, Vousdoukas et al 2018, Dottori et al 2018. As such, these studies have not accounted for compound events, in which the combination of multiple drivers and/or hazards can interact to modulate risk (Zscheischler et al 2018). Compound flood events can occur from the interplay between riverine and coastal flood drivers, for instance when: high sea levels propagate up a river leading to elevated water levels; and/or the drainage of high river discharge is impeded by elevated sea levels. Current riverine flood hazard models ignore these interactions and potential dependencies between riverine and coastal flood drivers, which may result in an under-or overestimation of flood risk (Wahl et al 2015. A first step towards accounting for compound events in global flood risk assessments is to understand where, and under which conditions, compound events modulate flood hazard. Several studies have addressed this by examining statistical dependence between different riverine and coastal flood drivers. They find dependence between: storm surge and precipitation in Australia (Zheng et al 2013, Wu et al 2017, the United States (Wahl et al 2015, Moftakhari et al 2017, Europe (Petroliagkis 2018, Bevacqua et al 2019, and the Netherlands (van den Hurk et al 2015, Ridder et al 2018); and storm surge and discharge in various parts of the United Kingdom (Svensson and Jones 2002, Lamb et al 2010, Hendry et al 2019, the Netherlands (Kew et al 2013, Klerk et al 2015, Khanal et al 2019, Texas (USA) (Couasnon et al 2018) and Italy (Bevacqua et al 2017). At the global scale significant dependence between storm surge and discharge based on observations was found at more than half of the locations studied  and based on simulations at 26% of the locations studied (Couasnon et al 2020).
A limitation of dependence-based analyses of compound events is the need for event selection based on the flood drivers (e.g. surge or discharge) rather than flood levels. This introduces bias in the joint probability estimate, as events are either conditioned on one driver or on the other (Hawkes 2008, Zheng et al 2014. Furthermore, extreme water levels might be driven by events that are not extreme themselves (Serafin et al 2019). Van den Hurk et al (2015) were the first to carry out an impact-based analysis of compound events (i.e. based on the impact of compound flood drivers rather than their dependence) for a case study of a near-flood event in the Netherlands. An ensemble of surge and precipitation timeseries were simulated with a regional climate model and used to force a hydrodynamic model of the inland water system. The simulated time-series were shuffled to remove dependence between surge and discharge. By comparing simulated water levels from original and 'shuffled' time-series, the effect of surgeprecipitation dependence on extreme inland water levels was examined. However, analysis of compound events based on simulated flood levels rather than flood drivers requires models that realistically simulate interactions between multiple drivers.
At the global scale, the first river routing model to account for surge-discharge interactions was presented by Ikeuchi et al (2017). They included dynamic downstream sea level conditions in the global river routing model CaMa-Flood (Yamazaki et al 2011) by coupling it to the Global Tide and Surge Model (GTSM; Muis et al 2016). They show a significant difference in the annual maxima of riverine water levels between simulations using dynamic sea level boundary conditions and those using static mean sea levels. However, they did not assess the drivers of extreme water levels nor the effect of surge on flood levels specifically, leaving the question unanswered as to where, and to what extent, compound surge affects flooding.
To date, no global analysis of surge-discharge interactions based on simulated water levels exists. To fill this gap, we developed a global compound flood model framework with the aim to identify dominant flood drivers in deltas globally and assess the effect of surge on riverine flood hazard and impact. This is an important step towards including compound flood events in global flood risk modelling.

Methods
We developed a model framework consisting of a global river routing model forced by a multi-model ensemble of global hydrological models and bounded downstream by a global tide and surge model (section 2.1). We analyzed simulated water levels from the model framework to classify the dominant driver of riverine flooding in deltas globally (section 2.2); to assess the effect of surge on flood hazard (section 2.3); and flood impact in terms of population exposed (section 2.4). . These runoff and dynamic sea level (surge and tide) data were used to force the global hydrodynamic river routing model CaMa-Flood (Yamazaki et al 2011) to simulate riverine water levels and flood depths which are input for the analysis. Each model component is further discussed in this section.

Model framework
We used CaMa-Flood version 3.6.4, which has a 1D routing scheme derived from HydroSHEDS (Lehner et al 2008). CaMa-Flood has an explicit representation of floodplains (Yamazaki et al 2011), which is crucial to correctly simulate discharge extremes (Zhao et al 2017). It solves the local inertial equation ( We ran CaMa-flood for the period 1980-2014, with a spin-up period of two years using repeated forcing from the first year. Daily values of instantaneous discharge and water levels at 00:00 GTM are stored and used as input to the analysis in this study. CaMa-Flood is forced by runoff, which is introduced at the head of each river channel section. The runoff data are obtained from the state-of-the-art global multi-model 15 ′ resolution ensemble dataset of E2O tier 2, which serves as a state of the art in current global hydrological modelling. The multi-model range in runoff stems from a combination of different total evaporation values and different storage dynamics in the models due to the different concepts and parameterization of runoff generation, see table 1, representing the uncertainty in land surface and hydrological processes (Schellekens et al 2017). For more details about the individual E2O models we refer the reader to Dutra et al (2017) and Schellekens et al (2017). From the available models we selected five that assume natural conditions, i.e. without anthropogenic water extractions. The runoff data were preprocessed to be on an identical grid from 90 North to 60 South, re-defined as a positive flux, and negative runoff values were set to zero in the JULES and ORCHIDEE data after discussions with the data owners (personal communication, 2018). Runoff data are then interpolated to the CaMa-Flood unit catchments using mass-conservative area-weighted averaging. We validated simulated discharge from CaMa-Flood forced by the E2O runoff ensemble against observations from the Global Runoff Data Centre with a focus on the magnitude and timing of discharge extremes. Although we find a large spread between individual models, the ensemblemean performance statistics generally shows low model bias and small time lags compared to observations (see supplementary information (available online at stacks.iop.org/ERL/15/104007/mmedia)).
We introduced dynamic sea level boundary conditions from GTSM at the downstream end of each river in CaMa-Flood. GTSM is the first global hydrodynamic model to simulate surge levels, i.e. the response of the sea surface to changes in atmospheric pressure and wind speed (Pugh and Woodworth, Pugh et al 2014), with sufficiently high temporal and spatial resolution for this application (i.e. near-shore resolution of 2.5 km). It has good performance compared to tide gauge data and other models , Wahl et al 2017, Cid et al 2018 and the timing and magnitude of storm surge peaks display sufficient performance for global scale compound flood analysis (Couasnon et al 2020). FES2012 simulates tides based on 32 tidal constituents and assimilation of satellite altimetry data (Carrere et al 2012) and is proven to have good near-shore performance (Stammer et al 2014). Mean sea level, tide, and surge are linearly superimposed to yield time-series of total still water levels at a 30-minute temporal resolution, thereby ignoring non-linear surge-tide interactions.  (2017). CaMa-Flood and GTSM do not have a perfectly joined interface: the most downstream river point in CaMa-Flood (hereafter referred to as river mouth) is often located inside the estuary, whereas GTSM output locations are slightly offshore. We therefore assumed a simplified estuary to schematize the missing link between the CaMa-Flood river mouth and GTSM. As the exact shape and bathymetry of estuaries globally is unknown, we extrapolated the channel width and depth from the CaMa-Flood river mouth, keeping the depth constant (Savenije 2005) and with a set length of 10 km. This estuary channel length is based on extensive validation by Ikeuchi et al (2017). River mouths in CaMa-Flood were coupled to the nearest GTSM output location within a maximum distance of 75 km. This distance threshold was selected as a trade-off between including as many river mouths as possible and excluding unrealistic links with GTSM output locations. Due to the relatively coarse resolution of the hydrological models, we focused on catchments with a minimum catchment size of 1000 km 2 . Using these criteria, a downstream boundary was set for 3433 river mouths based on 2352 GTSM output locations.

Flood drivers
We classified the dominant drivers of flooding at each river mouth location, represented by annual maximum riverine water levels (h AM ) extracted from the simulated time series, into four classes: surge-dominant, discharge-dominant, compounddominant or insignificant. The classification is based on the rank correlations between both h AM and discharge and h AM and skew surge (i.e. vertical difference between maximum still water level and high tide in a tidal cycle, stored as the maximum value at a daily time step). We used skew surge as it is the quantity of total water levels that might lead to flooding (

Flood levels
We developed three experiments, see table 2, to assess the difference in extreme riverine water levels with (scenario A) and without (scenario B) surge components included in the dynamic downstream sea level boundary. Surge levels were separated into a daily and seasonal component to assess their relative effects on flood levels. The seasonal component is associated with seasonal gyre circulation driven by synoptic pressure and wind differences at time scales longer than one month (e.g. Yang et al 1998, Palma et al 2004 and computed as monthly mean surge levels. The daily component is associated with surge due to short term meteorological variations in wind speed and sea level pressure and is computed as the difference between the total variation and seasonal component. In order to derive return periods of extreme water level beyond the length of our simulate time series, we fitted the 2-parameter Gumbel distribution using the L-moments method (Hosking and Wallis 2005). We find a difference in flood level to be significant if the sign of the difference is the same for all ensemble members. Additionally, confidence intervals (5th-95th percentiles) are obtained from bootstrapping with a sample size of 1000, where the Gumbel parameters are bias-corrected for the mean of bootstrap parameter samples.

Population exposed
We analyzed the population exposed to flooding by overlaying downscaled inundation and population maps at 18" resolution, assuming no flood protection. The downscaled inundation maps are calculated based on the HydroSheds elevation (Lehner et al 2008). Cells are flooded when the flood depth of a unit-catchment is larger than the relative height above the outlet elevation of that unit-catchment. We used the 2010 WorldPop 30" resolution gridded population dataset (Tatem 2017) and resampled it to the resolution of the inundation depth maps using bi-linear interpolation, i.e. linear interpolation in x and y direction, of population density. We assume that if flood depth is larger than zero the total population in that grid cell is exposed. Flood depths are underestimated if surge is ignored in basins where we find a positive difference in simulated flood depths between a scenario with compared to a scenario without surge levels, see experiment 1 in table 2. We find flood depths to be significantly underestimated if the difference is positive for all ensemble members. Finally, we construct a risk curve based on the population exposed and flood exceedance probability at return periods ranging from 1 to 100 years. Expected annual population exposed is then calculated as the area under the risk curve using the trapezoidal rule (e.g. Ward et al 2011). Results of the ensemblemean expected annual population exposed are presented.

Flood drivers
Globally, flood drivers are classified as compounddominant at 19.7% of the 3433 river mouth locations (figure 3(a)), with an average correlation of 0.57 between h AM and skew surge and 0.63 between h AM and discharge at these locations. Flooding is discharge-dominant at 69.2% of locations, with an average correlation of 0.84 between h AM and discharge at these locations, and surge-dominant at 7.8% of locations, with an average correlation of 0.60 between h AM and skew surge at these locations. The remaining 3.3% of locations are classified as insignificant as neither driver displays significant correlation in a majority of the ensemble members. Generally, compound flood drivers are found around large parts the USA, north-west Europe, the east coast of China, the east coast of Thailand and Malaysia, and around the Australian coastline. These regions are largely similar to those identified with high compound flood potential based on statistical dependence between simulated (Couasnon et al 2020) and observed  surge and discharge. Notable differences occur along the east coast of the USA and the coast of the Baltic sea, likely due to the different selection criteria for compound events between the studies. For the UK we find a similar spatial pattern of locations with compound drivers compared to locations with a frequent joint occurrence of high skew surges and high river discharge (Hendry et al 2019), which are found more often along the west and south coasts relative to the east coast of the UK. Next, we examined relationships between characteristics of river mouth locations and flood driver classification and examined whether these are   (2019) focus on high surge and high discharge, we also sample events with high surge and moderate discharge. Under these conditions, surge is more likely to propagate up rivers with flat topography.

Flood levels
Our results show that 1-in-10 years (T10) flood levels are generally exacerbated due to surge, with an overall ensemble-mean difference in riverine water level at the river mouth (∆h) of 11 cm between Scenario A and B of experiment 1 in table 2. ∆h is significant and positive at 64.0% of the 3433 river mouth locations studied, and significant and negative at 12.2%, while at 23.9% the ensemble members do not agree on the sign of ∆h, and are therefore classified as insignificant ( figure 5(a)). Moreover, ∆h is larger than the 5%-95% bootstrap confidence intervals for all ensemble members at 17.3% of the locations. ∆h is largest at locations with surge-dominant (28 cm) or compound-dominant flood drivers (30 cm), while ∆h is small (3 cm) at locations with dischargedominant flood drivers. Generally speaking, regions with the largest positive ∆h are the coasts of Alaska (US), North-West Europe, the Chinese coast at the Yellow Sea, and the coast on the Gulf of Carpentaria (Australia), which are all characterized by large surge extremes. Generally, at higher return periods the number of locations with significant ∆h decreases, while the overall ∆h at locations with a significant difference increases, see table 3. To better understand ∆h, we separate it into a difference in riverine water level due to a seasonal (∆h seasonal , see experiment 2 in table 2) and daily component (∆h daily , see experiment 3 in table 2). The daily component is mainly associated with surge due to short term meteorological variations in wind speed and sea level pressure, while the seasonal component is associated with seasonal gyre circulation (e.g. Yang et al 1998, Palma et al 2004. Large positive values of ∆h are mainly caused by ∆h daily , which for T10 is positive at 73.1% of the locations with a mean increase of 14 cm ( figure 5(b)). Negative ∆h is mainly caused by ∆h seasonal , negative at 50.3% of the locations, with a mean decrease of 3 cm ( figure 5(c)). In some areas, ∆h seasonal and  Ikeuchi et al (2017). who reported on the effect of total sea level variations on riverine water levels, we find similar areas with large ∆h. Notable differences between the two studies include the Gulf of Carpentaria and North Sea coast, where we find larger ∆h, which can be attributed to relatively high surge levels.

Population exposed
If surge is ignored, flood depths (and thus flood risk) are significantly underestimated for 30.7 million out of 332.0 million of the total expected annual population exposed (ensemble-mean), i.e. 9.3%. In absolute numbers, most people for whom flood depths are underestimated live along the densely populated coasts of east and south Asia. In relative numbers, flood depths are underestimated for a large percentage of the total expected annual population exposed (ensemble-mean) in small coastal basins with compound-or surge-dominant drivers, but also larger basins along the Hudson Bay coastline (Canada), the Neva (Russia), and the Elbe and Weser (Germany), see figure 6.

Limitations of the datasets and methods
The magnitude and timing of annual maxima surge and discharge estimates from GTSM and CaMa-Flood are not perfectly resolved, see section 2.2. To account for some of these uncertainties, we used the   E2O tier 2 multi-model ensemble. We only used a single surge model as there is less uncertainty in the timing of surge compared to discharge simulations (Couasnon et al 2020) and to date there is only one global hydrodynamic surge model with sufficient temporal and spatial resolution for this application.
Some processes that could affect the classification of flood drivers are currently missing in the model framework. GTSM does not account for non-linear surge-tide interactions, or inter-annual variability in mean sea levels due to steric effects or waves, which can be important drivers of coastal flooding at regional scales (e.g. Arns et al 2017, Vitousek et al 2017. CaMa-Flood does not take the operation of reservoirs into account, while these will significantly change the magnitude and timing of discharge peaks (Mateo et al 2014, Fleischmann et al 2019. Local variations in bathymetry that are not addressed in the CaMa-Flood and/or GTSM models may cause bias in the absolute water levels locally. Near-shore and estuarine areas are still very difficult to resolve accurately in global bathymetry datasets (Weatherall et al 2015) and therefore provide large uncertainty for global compound flood risk analysis. The model framework does not account for the influence of discharge on local sea levels as these are derived independently. A two-way coupling between GTSM and CaMa-Flood would be required to assess the complete interactions.
Furthermore, we did not account for uncertainties in the meteorological forcing. While the MSWEP V1.2 precipitation dataset is known to have a good performance compared to many other stateof-the-art global precipitation datasets, it has some caveats, including spurious drizzle and attenuated peaks (Beck et al 2017). GTSM is known to underestimate surge in areas with tropical cyclones due to the coarse spatial resolution of ERA-Interim (Muis et al 2016, Dullaart et al 2020. This might lead to an underestimation of the contribution of surge to riverine flooding in areas with high cyclone activity. The classification of flood drivers is less sensitive to this underestimation as it is based on the relative rank of the flood drivers. Only if different annual maximum flood events are selected as a result of this underestimation, it may affect the classification. In general, recent updates of meteorological forcing datasets providing higher resolution data and/or longer timeseries, including MSWEP v2 (Beck et al 2018) and ERA5 (the successor to ERA-Interim), could further improve the robustness of our results.
We estimated riverine flood extent and subsequent flood impact based on downscaled flood depths from CaMa-Flood. To focus on the effect of surge on flood impact we assumed no flood protection as accurate global data on protection standards are sparse (Scussolini et al 2016) while simulated flood impacts very sensitive to flood protection (Ward et al 2013). Furthermore, direct coastal or pluvial flooding could further influence the simulated water levels. To resolve complex hydrodynamic interactions between different flood drivers in coastal areas, higher resolution 2D flood models are required. A nested modelling approach (e.g. Hoch et al 2019) could be a possible avenue to explore in order to take a more integrate approach to improve flood modelling in coastal areas without compromising too much on computationally efficiency.

Conclusions and future work
In this study we present the first mapping of the dominant drivers of riverine flooding in deltas globally and assessed the effect of surge on riverine flood hazard and impact. The research highlights the importance of including dynamic sea level boundary conditions in riverine flood risk models. Drivers of riverine flooding are compound-dominant at 19.7% of the locations analyzed, discharge-dominant at 69.2% and surge-dominant at 7.8%. Compared to locations with either surge-or discharge-dominant flood drivers, locations with compound-dominant flood drivers generally have larger surge extremes and are in basins with faster discharge response and/or flat topography. Globally, surge exacerbates T10 flood levels at 64.0% of the locations analyzed, with a mean increase of 11 cm. While this increase is the largest at locations with compound-or surge-dominant flood drivers, surge also affects flood levels at locations with discharge-dominant flood drivers. A small decrease in T10 flood levels is observed at 12.2% of locations analyzed due to negative surge levels associated with dominant seasonal gyre circulations. Finally, we show that if surge is ignored, flood depths are underestimated for 30.7 million out of a total of 332.0 million (9.3%) expected annual population exposed (ensemble-mean).
In general, large scale flood risk studies would improve from a more holistic representation of flooding in our models, including direct coastal flooding from storm surges and waves as well as pluvial and fluvial flooding. This may require more detailed 2D hydrodynamic modelling in coastal areas to resolve complex hydrodynamic interactions between these different drivers. While we focused on classifying the drivers of riverine flooding per location, investigating the drivers and meteorological conditions of individual flood events would further enhance our understanding of compound events.

Data Availability
The data that support the findings of this study are openly available. The source code for the simulation, pre-and postprocessing and the analysis is available on GitHub at https://github.com/Dirk Eilander/compound_hotspots (DOI: 10.5281/zenodo.3665811). The dataset of simulated water levels and discharge at 3433 river mouth locations globally, including several components of nearshore still water levels is available on Zenodo (DOI: 10.5281/zenodo.3665734).
Research (NWO) in the form of a VIDI grant (Grant No. 016.161.324) and the TOUGOU program by MEXT Japan (JPMXD0717935457), and JSPS KAKENHI (Grant No. JP16J07523).