Dependence between high sea-level and high river discharge increases flood hazard in global deltas and estuaries

When river and coastal floods coincide, their impacts are often worse than when they occur in isolation; such floods are examples of ‘compound events’. To better understand the impacts of these compound events, we require an improved understanding of the dependence between coastal and river flooding on a global scale. Therefore, in this letter, we: provide the first assessment and mapping of the dependence between observed high sea-levels and high river discharge for deltas and estuaries around the globe; and demonstrate how this dependence may influence the joint probability of floods exceeding both the design discharge and design sea-level. The research was carried out by analysing the statistical dependence between observed sea-levels (and skew surge) from the GESLA-2 dataset, and river discharge using gauged data from the Global Runoff Data Centre, for 187 combinations of stations across the globe. Dependence was assessed using Kendall’s rank correlation coefficient (τ) and copula models. We find significant dependence for skew surge conditional on annual maximum discharge at 22% of the stations studied, and for discharge conditional on annual maximum skew surge at 36% of the stations studied. Allowing a time-lag between the two variables up to 5 days, we find significant dependence for skew surge conditional on annual maximum discharge at 56% of stations, and for discharge conditional on annual maximum skew surge at 54% of stations. Using copula models, we show that the joint exceedance probability of events in which both the design discharge and design sea-level are exceeded can be several magnitudes higher when the dependence is considered, compared to when independence is assumed. We discuss several implications, showing that flood risk assessments in these regions should correctly account for these joint exceedance probabilities.


Introduction
Between 1980 and 2016, global flood losses are estimated at ∼$1 trillion, with an estimated ∼215 000 fatalities (Re 2017). These impacts are particularly hard felt in low-lying, densely populated, deltas and estuaries (Tessler et al 2015), whose location at the land-sea interface makes them naturally prone to flooding. In these regions, when coastal and river floods coincide their impacts can be worse than when they occur in isolation. These are examples of 'compound events', defined by Zscheischler et al (2018) as '...the combination of multiple drivers and/or hazards that contributes to societal or environmental risk'. In 2017, a combination of unprecedented local rainfall intensities and storm surges from Hurricanes Harvey, Irma, and Maria, led to major flooding in Houston, Florida, and numerous islands in the Caribbean, highlighting the devastation that compound floods can cause (Dilling et al 2017, Wahl et al 2018, Zscheischler et al 2018. For decision-making and management, the simultaneous occurrence of high discharge and sea-levels is important for designing flood protection infrastructure and drainage in regions where total water levels are influenced by both. If independence is assumed, the joint exceedance probability can be underestimated (e.g. Wahl et al 2015), whilst assuming full dependence could lead to overdesign. Lagged occurrences of river and coastal floods can also influence overall inundation extent and area, if the time between them is too short for initial floodwaters to recede, like during the 2011 Thailand floods (Trigg et al 2013). Therefore, hazard mapping in deltas and estuaries that ignores coastalriver interactions can underestimate flood extent and depth. Moreover, regions can become more vulnerable in the immediate aftermath of a preceding natural disaster (e.g. Budimir et al 2014, Gill andMalamud 2014).
Whilst there is a rapidly growing recognition of the importance of understanding compound events in both in the scientific (Wahl et al 2018, Zscheischler et al 2018 and decision-making communities (UNISDR 2015), until recently they received relatively little scientific attention. Locally, there are several case-studies in Europe, Australia, USA, and China demonstrating that statistical dependence exists between the frequency or magnitude of coastal floods and rainfall or discharge (e.g. Loganathan et al 1987, Pugh 1987, Samuels and Burt 2002, Svensson and Jones 2002, 2004, van den Brink et al 2005, Hawkes 2008, Kew et al 2013, Lian et al 2013, Zheng et al 2014, Klerk et al 2015, van den Hurk et al 2015, Bevacqua et al 2017. At continental scale, two studies have examined dependence between storm surge and precipitation using observed datasets; Zheng et al (2013) in Australia and Wahl et al (2015) in USA.
However, to date, dependence between surge and river discharge has not been examined globally. Therefore, whilst case studies have shown the strength of this dependence in some locations, an overview of the larger geographical scale and magnitude does not exist, and hence discussions on their relevance for policy remain anecdotal. To address this, the objectives of this letter are: (1) to provide the first assessment and mapping of the dependence between observed high sea-levels and high discharge for deltas and estuaries around the globe; and (2) to demonstrate how this dependence may influence the joint probability of floods exceeding both the design discharge and design sea-level. The letter is intended to provide a first-cut analysis, to identify hotspots for compound flooding at the global scale; it does not systematically assess causal mechanisms or the resulting socio-economic risks.

Methods
We assess dependence between sea-levels and discharge using sea-level data from the Global Extreme Sea-level Analysis Version 2 database (GESLA-2) (Woodworth et al 2017) and observed discharge data from the Global Runoff Data Base (GRDB), supplied by Global Runoff Data Centre 8 . The research involves five steps: (1) selection of stations; (2) extracting time-series of high discharge and sea-levels; (3) correcting for flowtimes between discharge gauging stations and coast; (4) assessing dependence in the resulting time-series; and (5) assessing joint probability of floods exceeding both design discharge and design sea-level. These steps are described in the following subsections.

Datasets and selection of stations
We use mean daily discharge directly from GRDB. For sea-levels, we use hourly total sea-levels from GESLA-2, from which we extract daily maxima. We also use a time-series of daily maximum skew surge, i.e. vertical difference between predicted and observed high water in a tidal cycle (which may be with time-lag), extracted following Haigh et al (2016). We then extract station combinations from GRDB and GESLA-2 that satisfy the following criteria: (1) minimum of 20 years overlapping data; (2) minimum completeness of 75% per year in both databases; (3) minimum upstream basin area least 1000 km 2 ; (4) maximum Euclidean distance of 500 km between discharge station and tide gauge; and (5) tide gauge and river basin outlet within maximum distance of one 0.5 • grid-cell from each other. For many tide gauges, this results in several discharge stations satisfying the criteria; in these cases the most downstream river stations are selected. Following this selection, there are 187 combinations of discharge stations and tide gauges (figure 1(b)), with mean length ∼39 years and median length 36 years.

Extracting time-series of high discharge and sealevels
From the daily time-series, we extract time-series of annual maxima. To assess the sensitivity of the results to different sampling methods, we also use peaks over threshold (POT) (i.e. all values selected above a given percentile threshold). For each method, we first identify annual maximum discharge values (or POT) for each year, and then select the highest sea-level (total or skew surge) within ±1 days of this event. We used this ±1 days window as we do not have information on the exact timing of the discharge events, and are interested in discharge/surge events that occur within approximately a day of each other, rather than on the same calendar day per se. Throughout this paper, this is referred to as sealevel conditional on high discharge (SL cond Q) or skew surge conditional on high discharge (SS cond Q). Secondly, we identified annual maximum sea-levels (or POTs) for each year, and then selected the highest discharge value within ±1 days of this event. Throughout this paper, this is referred to as discharge conditional on high total sea-level (Q cond SL) or discharge conditional on high skew surge (Q cond SS). For each of the time-series of annual maxima or POT, we also extract time-series of the other variable using time-lags from −5 to +5 days. For example, for SL cond Q we identify annual maximum discharge values (or POT-series) for each year, and then select corresponding sea-levels with time lags of −5, −4, −3, −2, −1, 0, +1, +2, +3, +4, and +5 days. For the POT method, we use both 95th and 99th percentiles. Events exceeding those percentiles are identified, whereby a 3 day window is used to ensure independent events (Haigh et al 2016). In other words, if discharge (or sea-level) exceeds the percentile threshold on more than 1 day within a 3 day period, this is identified as 1 event. Sensitivity was also assessed using 5 and 7-day windows; the results are very insensitive.
For each of these time-series, we calculate parameters of five marginal distributions: LogNormal, normal, exponential, Weibull, and generalized extreme value. Then, we identify the distribution best fitting each of the individual time-series, using a goodness-of-fit test comparing empirical and theoretical non-exceedance probabilities using root mean squared error. For more information on statistical methods see supplementary information available at stacks.iop.org/ERL/13/084012/mmedia.

Correcting for flow-times between discharge gauging stations and coast
As the discharge gauging stations used are up to 500 km from the coast, flow-times between gauges and coast can be several days. Therefore, we correct for this by estimating average flow-times between discharge gauging stations and the outlet in days, and then shift the daily time-series forward by this number of days. For example, if the flow-time is one day, we assume that observed discharge reaches the coastline one day later. Bankfull flow-time estimates are taken from Allen et al (2018), who applied a kinematic wave model to the global HydroSHEDS hydrography dataset, enhanced with information regarding river length, bankfull width, and bankfull depth. However, these estimates are not available for locations >60 • N. For those locations, we applied a lagged crosscorrelation analysis based on paired time-series of modelled discharge at the gauge and downstream river outlet to derive the flow-time. Discharge is modelled using the CaMa-Flood routing model (Yamazaki et al 2011) forced with specific runoff from the Watergap V2.2 WRR2 eartH2Observe re-analysis data (Müller Schmied et al 2016, Dutra et al 2017. We find the flow-time based on the lag with the highest correlation within a maximum window of 14 days. For stations where this approach returns a correlation coefficient >0.98, the number of days is either the same as in Allen et al (2018), or different by one day, for 83% of stations. For stations where the correlation coefficient is <0.98, we manually assign a 0-day flow-time as these locations are within ∼30 km of the coastline.

Assessing the dependence
We measure dependence between discharge and total sea-level/skew surge using Kendall's rank correlation coefficient (Kendall 1938). Dependence is assessed for each of the lag-times in section 2.2, and significance assessed using = 0.10 due to the relatively low number of years at some stations. We also carried out the analyses using = 0.05.
For each station combination with statistically significant dependence, we apply copula theory (De Michele and Salvadori 2003, Grimaldi and Serinaldi 2006, Nelsen 2006) to assess the dependence structure. We use the maximum pseudo-likelihood estimator (Kojadinovic and Yan 2010) to calculate parameters for three different copulas: Gumbel (upper tail dependence), Frank (no tail dependence), and Clayton (lower tail dependence). We then use these parameters to simulate 1000 pairs of ranks for each copula model. The copula model most suitable for simulating the dependence structure for each station combination is selected by comparing non-parametric tail dependence coefficients (Schmidt and Stadmüller 2006) above a threshold of 0.8, derived from the observed and simulated rank pairs. Finally, we use the Cramer-von-Mises test (Genest et al 2009) to assess goodness-of-fit between observations and simulations. In this letter we only apply copula models for which the null hypothesis is accepted that the samples are from the same underlying distribution.

Assessing the joint probability of floods exceeding design discharge and design sea-level
For station combinations with significant dependence, and a copula model whose simulated values are statistically similar to observed values, we assess the joint probability of floods exceeding the design discharge and design sea-level simultaneously or in close succession. To do this, we need an estimate of the design standard of river and coastal flood protection measures for the river stretch or coastline closest to each discharge station and tide gauge. For river flood protection, we use protection standards from the FLOPROS database (Scussolini et al 2016). For coastal flood protection, standards are derived from the DIVA database (Hallegatte et al 2013, Sadoff et al 2015). Given the high uncertainty of protection standards in FLOPROS and DIVA, we also assess the joint probability of sea-levels and discharge exceeding a return period of 10 years.
We use the marginal distributions identified in section 2.4 to estimate the probability of discharge or sea-level exceeding design levels. Then, we calculate the joint exceedance probability assuming: (1) full independence, for which the joint exceedance probability is the probability of discharge exceeding design discharge multiplied by the probability of sea-level exceeding design water level; and (2) dependence, as defined by the copula model. Finally, we calculate the factor difference in the joint exceedance probability when using the dependence compared to independence assumption.

Dependence of discharge and sea-levels globally
First, we examine dependence between discharge and skew surge using a 0-day lag-time. We find statistically significant dependence for 41 stations (22% of stations) for SS cond Q and 67 stations (36%) for Q cond SS; values are shown in figure 2. Therefore, there are more locations with statistically significant dependence for Q cond SS than SS cond Q; this is especially the case in UK and Japan. However, across all locations there is no clear signal of difference in the strength of dependence between the SS cond Q and Q cond SS cases; values are higher for SS cond Q than Q cond SS at 52% of locations. Moreover, the mean values for SS cond Q and Q cond SS are 0.09 and 0.08 respectively (no significant difference; t-test, p = 0.366).
In figure 3 we show dependence for the lag-time (−5 to +5 days) for which is highest; corresponding lag-times are shown in figure S1. When lags are included, we find statistically significant dependence between discharge and skew surge for 104 stations (56%) for SS cond Q and 101 stations (54%) for Q cond SS. The results show that it is more common for high discharge events to follow high surge events than vice versa; for SS cond Q, this occurs in 70% of locations (i.e. locations with negative lag in figure S1), and for Q cond SS, this occurs in 76% of locations (i.e. locations with positive lag in figure S1).
Comparing results for skew surge (figure 3) with total sea-levels (figure S2), geographical patterns are similar, with regional differences. For example, in western Britain our skew surge results (figure 3) show significant dependence for most stations, yet this is not the case for SL cond Q. Here, where the tidal range is large, the tidal influence on overall sea-level is stronger than in locations with low tidal range. Haigh et al (2016) showed that most extreme sea-levels around the UK are generated by moderate storm surge coinciding with high spring tides; this is probably the reason dependencies are lower for total sea-level.
We examine sensitivity to using the POT method for event selection. Results are shown for skew surge with thresholds of 95% (figure S3) and 99% ( figure  S4). Using a 99% threshold, results are very similar to those using annual maxima. Using a 95% threshold, the number of significant results decreases, as do the values. We also carried out the analyses using a confidence limit of = 0.05 to assess the significance of the values; the results for dependence for 0-day timelags (figure S5) and with time-lags (figure S6) are very similar to those using = 0.10.

Regional patterns of dependence in discharge and sea-levels globally
In this section, we discuss regional patterns of dependence, and interpret results with reference to existing literature. The discussion is limited to regions with a large number of stations, as it is not possible to identify general patterns for those regions with limited observations.
For the west coast of the USA, we find significant dependence between high discharge and skew surge (figure 3) and total sea-level (figure S2) at most stations. For two locations (La Jolla, California, in the southwest, and Neah Bay, Washington, in the northwest) we show composite analyses of atmospheric conditions on the dates of annual maximum discharge and annual maximum skew surge in figures S7 and S8 respectively. The data for the composites are derived from the NOAA-CIRES 20th Century Reanalysis Version 2c dataset (NOAA 2015), which has a 2 • × 2 • spatial resolution and four times daily temporal resolution from 1851-2014. To cover post 2014 events, we use the NCEP-NCAR reanalysis dataset (Kalnay et al 1996), interpolated from its native 2.5 • × 2.5 • resolution to 2 • × 2 • . Figures S7 and S8 are representative of atmospheric conditions in most of the studied locations along the US West Coast, and show that atmospheric conditions are similar on the dates of annual maximum discharge and skew surge. Moreover, most of the basins studied on this coast (except the Columbia River) are relatively small and steep with fast catchment response times. Together, these aspects can  . for dependence between annual maximum discharge and skew surge for: (a) SS cond Q, and (b) Q cond SS. Locations with black dots denote no significant dependence ( = 0.10). Results shown here for highest dependence ( ) over all time-lags (from −5 to +5 days). Regional panels show Japan, USA, North America, and northwestern Europe. explain the similar patterns in dependence for SS cond Q and Q cond SS. In contrast, Wahl et al (2015) found few locations with significant dependence between surge and precipitation along this coastline. Our composite analyses do show elevated precipitable water content (especially for Q cond SS) in the catchments upstream from the tide gauges (albeit relatively weak for La Jolla for SS cond Q). The west coast of the USA is typified by a topography in which mountain ranges lie close to the coast, which can lead to orographic rainfall. In Wahl et al (2015), precipitation gauges were used if they are within 25 km radius of the tide gauges. Hence, Factor difference in the joint exceedance probability of discharge and sea-level exceeding design level (with zero lags), when calculated using the dependence structure derived from the copulas compared to assuming independence, shown for (a) SL cond Q, and (b) Q cond SL. Locations with black dots denote no significant dependence ( = 0.10) or no significant copula fit. Inset histograms show the number of stations for which different change factors are calculated. Regional panels show Japan, USA, North America, and northwestern Europe. these may be located closer to the coastline than the first mountain ranges, meaning that coastal rainfall may not show dependence with sea-level, but that river discharge does. In eastern USA, our results for discharge and skew surge are largely similar to those of Wahl et al (2015). We find significant dependence for most station combinations, although for Q cond SS we find no significant dependence for stations in northeastern USA. Figure  S9 shows a composite analysis of atmospheric conditions on the days of annual maximum skew surge for Portland (Maine). The figure shows surge events driven by strong low pressure systems to the southeast of the tide gauge, causing strong easterlies and northeasterlies. Whilst precipitable water content is elevated for this location, there is no strong dependence with Factor difference in the joint exceedance probability of discharge and sea-level exceeding design level (with lags of −5 to +5 days), when calculated using the dependence structure derived from the copulas compared to assuming independence, shown for: (a) SL cond Q, and (b) Q cond SL. Locations with black dots denote no significant dependence ( = 0.10) or no significant copula fit. Inset histograms show the number of stations for which different change factors are calculated. Regional panels show Japan, USA, North America, and northwestern Europe.
the discharge of the Kennebec River (to the north). Further investigation shows that 86% of the annual maximum skew surge events occur between October and February, whilst peak discharges occur in the boreal spring. Hence, seasonality has an important influence on dependence in such rivers with strong seasonal discharge and surge characteristics; very similar composite analysis plots are found for other locations in the region.
For the Gulf Coast, Wahl et al (2015) found significant dependence for most stations. For discharge, we only find this relationship at one station combination in each case (figure 3). For the combination of Grand Isle (near New Orleans)/Mississippi River, this can be explained by catchment size (>3 000 000 km 2 ). The other basins studied here are also relatively large, except the Guadeloupe River (ca. 13 000 km 2 ) and Sabine River (ca. 24 000 km 2 ); the latter are the rivers for which we find dependence. Moreover, highest surges in the region typically occur during hurricane season, when the discharge of the region's rivers are at their lowest.
Results for western Japan show strong dependence, especially for Q cond SS and Q cond SL, with strongest dependence for either 0-day or small positive time-lag (i.e. high surge or high total sea-level leading discharge). Similarly to the west coast USA, this region is typified by steep topography and relatively short distances to the coastline, meaning that precipitation from surgecausing storms can be enhanced by orographic effects (Negri and Adler 1993) and the resulting discharge can rapidly reach the coast.
In Europe, several local studies have been carried out in the UK and the Netherlands. For stations in the western Netherlands, we find no statistically significant dependence. Klerk et al (2015) found significant dependence between surge at Hoek van Holland and discharge at Lobith on the Rhine for a 6-day lag, which exceeds the maximum lag in our study. They relate the lagged response to precipitation delay between the Dutch coast and upstream in Germany, and typical travel times of the Rhine of four days.
For the UK, where rivers are shorter, we find many sites with statistically significant dependence. On the south coast, the Environment Agency (2000) found simultaneous dependence between discharge and total sea-level at Brockenhurst, Hampshire. In this region, our 0-day lag results for Portsmouth/Stour river (west of Brockenhurst) (figure S10) also show relatively strong dependence (SL cond Q: = 0.358, p = 0.014; Q cond SL: = 0.284, p = 0.043). More broadly for Britain, for Q cond SS we find significant dependence in most regions ( figure 3(b)), except the southeast. Similarly, Svensson and Jones (2002) found significant dependence of discharge conditional to surge north of the Firth of Forth and little dependence southwards, particularly south of the Humber Estuary. They state that during winter, high-surge events are often associated with cyclones to the north of Scotland, which are associated with precipitation from a southwesterly/westerly circulation. As the area north of the Firth of Forth is not sheltered from south-westerly winds by major topographical barriers, orographic rainfall can be enhanced when these airflows meet the hilly terrain north of the Firth, leading to high discharge. To the south, southwesterly winds have already encountered the Southern Uplands, Pennines, or Welsh mountains. Moreover, in the southeast, catchments are generally permeable, and therefore respond slowly to rainfall. In western Britain, Q cond SS shows significant dependence at all stations (figure 3(b)), with highest dependence for positive lag-times (i.e. when high surge precedes high discharge; figure S1). Svensson and Jones (2004) also found significant dependence at stations along most of this coastline. Most of the station combinations studied here are south or west facing, with fast catchment response times and orographically precipitation, meaning that discharge can arrive at the coast on the same day or shortly after a high surge event. For SS cond Q, we find significant dependence along large parts of the coastline, except for several western stations and one station in northeastern Scotland. Using composite analyses, Hendry et al (2018) showed that the storms that generate high discharge at eastern UK sites differ from those generating large surges. Discharge causing storms on the East coast typically tracked over the centre of Britain, with the low-pressure system being located over central Britain at the time of peak discharge. Meanwhile storms generating surges tracked north of Scotland, with the low-pressure system located over Scandinavia at the time of peak surge.
For Australia, we find significant dependence between discharge and total sea level for the majority of stations (figure 3). In Victoria, there are several stations where the dependence is not significant. This can probably be explained by lack of dependence in the storm-causing atmospheric processes in this region, since Zheng et al (2013) found significant dependence between surge and precipitation in most areas of Australia, but not in those parts of Victoria for which we find no dependence between discharge and skew surge.

Effects on joint exceedance probabilities
An understanding of the simultaneous joint occurrence is particularly important for designing flood protection infrastructure and drainage, whereas lagged dependence can be relevant for hazard mapping and planning emergency responses. In figure 4 we show the difference in joint exceedance probability of discharge and total sea-level exceeding design level on the same day, assuming dependence versus independence. In figure  5 we show the same results for the ranked correlations with the best lag-time. We only show results where there is significant dependence and where we identify a copula model with significant goodness-offit; copula models used are shown in figures S11 and S12 respectively. The results show the importance of considering the dependence structure on the estimated joint exceedance probability. For SL cond Q, significant dependence and copula models are found for 20% of stations for simultaneous occurrence and 55% of stations for lagged occurrences. For Q cond SL, significant dependence and copula models are found for 25% of stations for simultaneous occurrence and 42% of stations for lagged occurrences. The differences in joint exceedance probabilities when assuming dependence versus independence are large. For those station combinations for which significant dependence and copula models are identified, the factor difference exceeds +2 in 68% of cases for SL cond Q and 74% of cases for Q cond SL for the zero lag analysis, and in 59% of cases for SL cond Q and 67% of cases for Q cond SL for the lagged analysis. In these locations, compound river and coastal flood events are more than twice as likely to occur than if they were independent. To assess the sensitivity of the results to the estimated flood protection standards in FLOPROS and DIVA, figures S13 and S14 show the same results, but using a return period of 10 years for both river and coastal protection. The spatial patterns are robust, though the factor differences are generally lower, since the protection standards are greater than 10 years in most of the locations studied.
When assessing flood risk, it is vital to know flood impacts across a whole spectrum of exceedance probabilities. In figure 6, we show two examples of differences in joint exceedance probability, assuming dependence versus independence, for different combinations of exceedance probabilities of the two marginal distributions (Santa Monica and Santa Clara River in California, USA; and Komatsushima and Yodo River Japan); note that the colour bars for the two sub-figures have different scales. For both examples, the (non)-inclusion of dependence has an increasing influence on the joint exceedance probability, as the marginal exceedance probabilities become lower. Hence, when assessing the influence of compound floods in a risk assessment, it is not sufficient to assume a linear shift in the exceedance probability loss curve, but the change in exceedance probability should be assessed along the entire curve. Moreover, this shows that differences in joint exceedance probability for the dependence versus independence assumptions increase as the marginal exceedance probabilities decrease (i.e. with higher return periods). This is an important finding for flood risk management in a changing world where flood protection standards are generally being increased.

Outlook and concluding remarks
We provide the first global mapping of dependence between observed high sea-levels and high river discharge. Significant dependence exists at more than half of the station combinations studied. Therefore, this dependence is of broad importance geographically. This has large implications for flood risk assessments in deltas and estuaries influenced by both coastal and riverine processes. We show that the joint exceedance probability of events in which both design discharge and sea-level are exceeded can be much higher when dependence is considered. Hence, flood risk assessments in these regions should correctly account for these joint exceedance probabilities.
Whilst we do not systematically assess causal mechanisms behind the dependencies, we do identify several possible important mechanisms. Clearly, atmospheric conditions are of key importance. In some regions, high discharge and high sea-levels stem from the same atmospheric conditions-leading to similar dependencies when examining discharge conditional on high sea levels compared to sea levels conditional on high discharge-as is shown for the composite analyses for the western USA. However, in other regions high discharge stems from different atmospheric conditions than surge-causing events (as described here for northeastern Scotland), such as atmospheric rivers or frontal precipitation (De Luca et al 2017). Moreover, where discharge shows strong seasonality, this influences the probability of a high surge-causing event coinciding with the highest annual discharge. This is shown in this paper for the example of the Kennebec River and Portland in the USA. For several regions, we also show the importance of orographic effects and catchment response time. The latter can be influenced by human interventions, such as river channelisation (e.g. Zawiejska and Wyzga 2010). Other factors that could influence the timing of dependence between high discharge and sea levels include geology and land use, which influence permeability and soil moisture holding capacity, and therefore the timing and magnitude of discharge peaks (e.g. Fenicia et al 2009).
A limitation of our analysis is the spatial bias of observed data towards several geographical regions. Hence, we cannot make general statements about dependence in most of Africa, South America, and most parts of Asia. Also, most of the station combinations have a relatively short time-span. A way to address these issues could be through modelling. Physically-based and stochastic models have successfully been set-up at local scales to assess dependence between different flood variables under present-day and/or future conditions (e.g. Kew  Results and publications of the ISIMIP project could be very helpful (www.isimip.org/). Such model outputs can provide spatially and temporally comprehensive time-series required to consistently and systematically assess the sensitivity of dependencies between surge and discharge to various drivers, and could help in classifying the kinds of basins or regions in which compound floods are important. Moreover, they could provide data for assessing the impacts of future climate change on compound flooding, which Motakhari et al (2017) recently showed to be very important in the USA. At the same time, it would expand the current research to those regions not covered by observation stations.