Introduction

Coastal and estuarine flooding can have devastating and long-term impacts, and global flood losses are estimated to be over USD 1 trillion since 1980 (McGranahan et al. 2007; Munich Re 2017). Globally, low-lying estuaries at risk of flooding support essential services including transport and energy infrastructure, water supply, and access (i.e. ports and harbours), with 21 of the world’s 30 largest cities and 339 million people located next to estuaries (Edmonds et al. 2020). When flooding occurs, defences and critical thresholds can be exceeded, whereby infrastructure can be damaged with little or no notice and lead to significant consequences (Environment Agency 2020). These consequences of flooding include human casualties and fatalities, damage to property and infrastructure, and degradation of natural spaces that are important for physical and mental wellbeing (Martin et al. 2020). In the UK, coastal flooding is rated as the second highest risk of civil emergency in the UK, after pandemic influenza (HM Government 2020), and has an annual cost of up to £2.2 billion for flood management and emergency response (Penning-Rowsell 2015).

Estuaries are located at the land-sea interface, and flooding can be driven by a combination of high sea-levels and river discharge (Ward et al. 2018). High sea-levels can occur due to astronomical high spring tides and can be further exacerbated when they co-occur with storms generating large surges and waves at the coast. Storms can also generate heavy precipitation and lead to high fluvial and pluvial flows. High sea-levels are known to increase upper estuary flooding through the so-called ‘backwater effect’, where fluvial waters are impeded due to high sea-levels (Maskell et al. 2013; Ahn et al. 2018; Guo et al. 2020; Zhang et al. 2020). Recent research has highlighted that these drivers rarely occur independently (Bray and McCuen 2014; Hendry et al. 2019; Robins et al. 2021). Therefore, the consequences of estuarine flooding can be exacerbated when oceanographic and fluvial drivers occur concurrently or in close succession (Bevacqua et al. 2020), and this is a phenomenon called ‘compound flooding’. Hurricane Harvey (2017), Irma (2017), and Maria (2017) are recent examples of compound events, where sustained storm surges coincided with unprecedented heavy rainfall and river discharge to cause flooding (Wahl et al. 2015). Water levels in southeast Texas were up to 3 m above predicted levels during Hurricane Harvey due to a short phase lag in maximum runoff, onshore wind stress, spring tides, and an extended storm surge (Wing et al. 2019; Valle-Levinson et al. 2020). A change in climatic trends and land-use, notably rapid urbanisation (Sebastian et al. 2019), and strong surge-river dependence (Santos et al. 2021) exacerbated the impacts of Hurricane Harvey. The co-occurrence or lagged occurrence of high sea-level and river discharge can influence total water levels and the depth and extent of estuarine inundation and must be considered in decision making processes.

A strong dependence between the magnitude of sea-level and extreme river flows indicates likelihood of compound flooding and has been found on the US coast, along the coast of Portugal, Madagascar, Taiwan (Couasnon et al. 2020), and for the UK in the South and West (Svensson and Jones 2004; Hendry et al. 2019; Couasnon et al. 2020). For the UK, the north-eastwards track of North Atlantic low atmospheric pressure systems during autumn and winter causes concurrent high sea-levels and precipitation in south/west-facing estuaries (Haigh et al. 2016). Elevated sea-levels and river discharges can be generated from the same storm through, for example, low atmospheric pressure systems generating a storm surge along with high precipitation (Bilskie and Hagen 2018), e.g. Storm Desmond, 4–6 December, UK (Matthews et al. 2018). Dependence is weaker on the east coast of Britain as storm surges and heavy precipitation events are driven by different storm patterns and characteristics (Svensson and Jones 2002). However, past analyses were limited by using daily mean river flow data where the coarse temporal resolution may not adequately capture extreme flow magnitudes on Britain’s west coast (Robins et al. 2018). Robins et al. (2021) examined the likelihood, interactions, timing, and behaviour of drivers of compound flooding events at sub-daily scales in two contrasting estuaries in Britain. The Dyfi estuary is a small, steep catchment on the west coast of Britain; half of the 937 skew-surge events which exceeded the 95th percentile, used as a threshold to identify events likely to contribute to flooding, occurred within a few hours of a similarly extreme fluvial peak. The Humber estuary is a larger, shallower catchment on the east coast where extreme fluvial and skew-surge peaks were less frequent; here only 15% of events co-occurred (Harrison et al. 2021; Robins et al. 2021). Robins et al. (2021) concluded that smaller catchments and those with a flashy hydrological regime require analysis of river flows at sub-daily scales.

Sub-daily river discharge data can be used to resolve hydrograph behaviour, and can be analysed with surge residual, skew surge, or total water level to explore dependence between sea-levels and river discharge (Lucey et al. 2021). Storm surges have a stronger link with meteorology than total water level, which is largely influenced by the variability of astronomical tides (Haigh et al. 2016), and are best for identifying the causal processes of flooding in estuaries (Svensson and Jones 2004). Surge residual is commonly used to analyse dependence between extreme storm surge and river discharge (Nasr et al. 2021; Santos et al. 2021); however, the residual can misrepresent the surge as even a small difference in the timing of the predicted and actual tide creates an ‘illusory’ surge (McMillan et al. 2011; Paprotny et al. 2016). Skew surge provides more certainty than the residual as it is independent of tides, which is important in regions such as the UK, where shallow shelf seas lead to tide-surge interaction processes that can affect sea-levels (Prandle and Wolf 1978; Williams et al. 2016). Therefore, skew surge can be a more reliable indicator of the meteorological component of sea-level (Batstone et al. 2013). The skew surge is the difference between the predicted astronomical high tide and nearest observed high water and has been used to analyse dependence across Europe (Paprotny et al. 2018; Hendry et al. 2019; Camus et al. 2021) and the USA (Moftakhari et al. 2019; Jane et al. 2022). It is clearly important to understand causal processes of compound events, but it is equally important to understand in practical terms whether flooding could occur in estuaries. It can also be useful to consider the magnitude of tides, which have a large influence on total water levels at the coast, particularly in the UK (Haigh et al. 2016). Compound events driven by storm surges and high river discharge will often only cause flooding during high spring tides (Hendry et al. 2019; Bevacqua et al. 2020; Rulent et al. 2021); therefore, analysis of total water levels at sub-daily timescales is also needed.

Compound events can be defined by identifying the statistical dependence or number of joint extremes between drivers such as storm surge and river flow, or total water level and river flow (Paprotny et al. 2018; Wahl et al. 2015; Ward et al. 2018; Bevacqua et al. 2020; Hendry et al. 2019; Couasnon et al. 2020). Regional and global studies of compound flooding events have assessed upper tail dependence using copula theory to characterise bivariate joint distribution (Bevacqua et al. 2020; Paprotny et al. 2018; Ward et al. 2018), or used conditional sampling pairs one extreme flood driver with another value, so that at least one (if not both) values are extreme (e.g. annual maxima or peaks over threshold) (Wahl et al. 2015; Moftakhari et al. 2017). A peak over threshold method can provide stable estimates of the compounding potential for high discharge and storm surge events (Jane et al. 2022), but has been shown to produce lower correlation coefficients between variables than a block maxima approach (Camus et al. 2021). Data and sampling methods are a subjective choice and can influence and support understanding the likelihood of compound coastal flooding (Lucey and Gallien 2021).

Given the widespread potential for compound flooding, there is a need to understand how the timing and magnitude of the interacting drivers affect compound hazards in estuaries to support forecasts and warnings, emergency response, and long-term management plans (Penning-Rowsell et al. 2000). Estuaries that are susceptible to compound events require accurate hazard mitigation strategies that consider multiple drivers at appropriate temporal and spatial scales, and these should be accounted for in flood warnings and timely evacuation orders to minimise impacts to coastal communities.

This research aims to extend analysis of the likelihood and timings of high sea-level and river discharge events at sub-daily timescales for different estuaries and catchment types across Britain. The research uses a 30-year period (from 1984 to 2013) of 15-min instantaneous output frequency sea-level and river-discharge data to (i) highlight where across Britain sub-daily scale co-occurrence analysis of flood risk is needed; (ii) identify the best methods to identify joint occurrence of high storm surges and high river discharges around the coast of UK; (iii) identify estuaries that are susceptible to compound events where compound flooding could occur; and (iv) identify seasonal trends in sea-level and river discharge which could potentially increase flooding at certain times of year. The results identify estuaries where sub-daily analysis of compound events is needed to resolve higher frequency fluvial-surge events, and can be applied to direct future research to understand how changing river and sea-level climates may influence susceptibility of these estuaries to compound events and potential flooding.

Methods

Data

Observed total water level (TWL) data were used to represent sea-level and obtained from 27 class-A tide gauges, which form part of the UK National A-Class Tide Gauge Network, from the British Oceanographic Data Centre (BODC) (Fig. 1 and listed in Table S1). Data was obtained for each tide gauge at hourly frequency from 1984 to 1992 and at 15-min frequency from 1992 to 2013. Pre-1992 data was interpolated to 15 min to align with the resolution of the river gauge data. TWL data that is flagged in the records as improbable, null, or interpolated were discarded to ensure only accurate observations of TWL were included in the analysis. Missing tide gauge data resulted in gaps in the time-series record for TWL, between 79 and 98% coverage (Mumbles, Barmouth, and Sheerness have lowest coverage), but ensure only reliable and observed values were used in the analysis. Some tide gauges are not close to a corresponding estuary, e.g. Barmouth tide gauge is located up to 20 km away from Dyfi estuary. To test whether these tide gauge records are a good representation of regional water level and storm surge characteristics at the estuary mouth, we assessed modelled hindcasts TWL records from the Coastal Dataset for the Evaluation of Climate Impact (CoDEC) (Muis et al. 2020). These are generated from the Global Tide and Surge Model (GTSM) at a 2.5-km coastal resolution, forced with ERA5 climate reanalysis. We selected the modelled grid point closest to the mouth of each estuary and compared this to the measured records. Up to 40 cm (root mean squared error) was seen in magnitude of total water level when tide gauge data and modelled water level from the estuary mouth simulated by CoDEC were compared (shown in Fig. S1 in Supplementary Material). Mean phase difference between observed total water level at the nearest tide gauge and modelled total water level at the estuary mouth for each gauge is between 0 and 1 h, which is unlikely to substantially influence how extreme river-sea-level events are identified and associated lag times (shown in Fig. S2 in Supplementary Material).

Fig. 1
figure 1

Location of 126 river gauges (orange circles) and 27 tide gauges (blue stars) used in this analysis. Red lines indicate pairing between tide gauge and river gauges. Main and secondary UK rivers are shown in grey (Ordnance Survey Open Rivers 2021). Further details of each gauge are provided in Table S1

Total water level was linearly detrended by subtracting the annual mean from each year to remove the effects of historical mean sea-level rise from the time series (Coles 2001; Stephens et al. 2020; Luxford and Faulkner 2020). Skew surge (S) was calculated from the detrended TWL, as the difference between maximum TWL and maximum predicted based on astronomical tidal constituents for every 12.42-h tidal cycle, regardless of phase lags between the two maxima (Mawdsley and Haigh 2016). This resulted in one S value per 12-h tidal cycle (approximately 706 per year) that represents the whole tidal cycle, and the time of maximum TWL was assigned to each S value to give an indication of when the hazard occurred. S is used as an indicator of the meteorological component of sea-level to understand the causal processes of compound flooding in estuaries. The residual surge was not used here as it can misrepresent hazards as the maximum non-tidal residual does not necessarily occur at the time of high water.

River discharge measurements (Q) at 15-min instantaneous output frequency were obtained from the Environmental Agency (EA), Natural Resources Wales (NRW), and the Scottish Environment Protection Agency (SEPA) for 265 rivers across Britain. Q measurements were taken from the most downstream, non-tidal gauges on each river. River gauges were isolated which recorded (i) 30 years of discharge from 01 January 1984 to 31 December 2013 (with no data gaps) and (ii) maximum discharge during this period that exceeded 50 m3/s to identify rivers that could potentially contribute to flooding. This set 30-year period (1984 − 2013) was selected to ensure consistency between gauges and to maximise the number of gauges with full data included in the analysis. These criteria isolated 126 gauges, located in 126 distinct estuaries, for our analysis and these are shown in Fig. 1 and listed in Table S1 in the Supplementary Material. These 126 gauges give good spatial coverage across Britain, apart from in southeast England where several river gauge records were shorter than 30 years or no peaks in the record exceed 50 m3/s. A fast Fourier transform and Lomb–Scargle periodogram (Lomb 1976) was applied to the 126 river records to confirm no tidal signal was present in the river discharge records. A 24-h running mean was applied to the Q record as a smoothing operation. River gauges were paired with the nearest tide gauge on the same coastline the estuary discharges on to (see Fig. 1).

Identifying Co-Occurrence Events

Five methods were used here to identify extreme values in the S and Q records, and their relative timings (Methods 1–5). Daily mean and daily maximum Q are commonly used variables to represent the magnitude of river discharge each day, and Methods 3–5 identify peaks in the Q and S record. A list of variables are provided in Table S2.

Method 1: Daily Means

Daily mean Q values were calculated from the 15-min instantaneous Q record for each river gauge record. Then daily mean values which exceeded the 3-month moving mean + 50th percentile threshold in each record were extracted for the analysis. This threshold was set to identify values above the mean which may contribute to flooding. As most storm surges last longer than a tidal cycle, the two S values per day were averaged. One daily value was used as sometimes only one S occurs per day. Daily mean skew surge which exceeded the 3-month moving mean + 50th percentile threshold was also extracted for the analysis. This threshold was set to identify S which may also contribute to flooding. If neither daily mean Q nor S exceeded their threshold then the values were discarded from the analysis so that only extreme values were considered. Daily mean Q values are available to download from the National River Flow Archive (https://nrfa.ceh.ac.uk/) and often used in co-occurrence analysis, so this method aimed to identify if this data is representative of river behaviour.

Method 2: Daily Maximums

Daily maximum Q values were calculated from the 15-min instantaneous Q record, and daily maximum values which exceeded the 3-month moving mean + 50th percentile threshold were extracted for the analysis. The corresponding daily maximum S value which exceeded the 3-month moving mean + 50th percentile threshold was extracted. If neither daily maximum Q peak nor S exceeded their threshold then the event was discarded from the analysis so that only extreme values were considered.

Method 3: Q Peaks Over Threshold (POT)

Peaks in the 15-min instantaneous Q record were identified, and the peaks which exceeded a threshold were extracted for the analysis—their peak magnitudes and timings denoted as PQ and TQ, respectively. The threshold was defined as the 3-month moving mean + the addition of different Nth percentiles calculated from river discharge at each gauge. Increasing N creates a higher threshold; fewer number of Q peaks are selected, so that only the most extreme peaks in river discharge remain. In the Supplementary Material, we present a sensitivity test (Fig. S3) where we varied the threshold based on N = 10th through to N = 95th percentiles, that justifies our choice presented here. Following the sensitivity analysis, the 3-month moving mean + 50th percentile of each record was set as this threshold, ensuring that at least 100 peaks were isolated in each record so that the subsequent correlation analysis was robust (see Fig. 2). This method identified all the largest and most prominent peaks in each record, which could potentially lead to compound events, i.e. not just the largest on a given day, or the most prominent or independent peaks, but also those that occurred close together from clustered or prolonged storms. Smaller magnitude peaks in the summer, which could potentially lead to flooding (Marsh and Hannaford 2007; Macdonald et al. 2010), were also selected due to seasonality in the 3-month moving mean. Some ultrasonic gauges used in larger estuaries resolve more peaks than on other gauges, and the relative thresholds are different for each catchment due to variability in baseflow and extremes; therefore, we have taken an approach (N = 50th percentile) that ensures at least 100 events are captured in each gauge. There was no upper limit set on the number of peaks that were selected; if the peak exceeded the threshold, then it was included in the analysis.

Fig. 2
figure 2

a River discharge (blue line) and isolated peaks (PQ and TQ, red dots) above the threshold 3-month moving mean (red solid line) + 50th percentile (red dashed line) and S (black cross) with the 3-month moving mean (black solid line) + 50th percentile (black dashed line) at the 66,011 Conwy river gauge. b Example of corresponding extreme S (PS and TS, black circle) in a 20.08-h hydrograph window (orange shading) around the maximum Q peak

The coincident maximum S value was then identified within a hydrograph window around each TQ. S values which exceeded the 3-month moving mean + 50th percentile threshold were extracted for the analysis. If no S value exceeded the 3-month moving mean + 50th percentile threshold within the hydrograph window, then the river peak was discarded from the analysis. Their magnitudes and timings denoted as PS and TS, respectively. The S threshold was set to a 3-month moving mean + 50th percentile which isolates positive S values just above 0 m. The hydrograph window is unique to each gauge and calculated here by pairing the gradient of the rising hydrograph limb of each selected peak with one of 30 normalised, idealised gamma curves, each representing a different hydrograph shape, which have a known duration and magnitude (Robins et al. 2018). A two-parameter gamma distribution is used to generate synthetic hydrographs to more accurately calculate duration flow exceeds a threshold for clustered river discharge peaks, peaks where baseflow does not return to zero, or peaks where data is missing on the rising or falling limb of the observed hydrograph (Patil et al. 2012; Robins et al. 2018; Cidan and Li 2020). A larger gradient indicates a steeper, flashier hydrograph, and this hydrograph would be paired with idealised hydrograph 1 (see Fig. 3a). The maximum storm gradient on the rising limb of each observed hydrography was calculated, representing the flashiness of each storm. Since the idealised hydrographs were normalised so their integral equalled 1, these were scaled up to match the magnitude/duration of each PQ exceeding the 3-month moving mean threshold. This method ensures that the hydrograph window is such that the S and Q events could actually overlap and hence interact. The average event duration across all events in each record was calculated as the hydrograph window to generate a representative duration. For example, the mean hydrograph duration and subsequent hydrograph window was 20 h 5 min at gauge 66011 (Conwy) (see Fig. 2b) with the 10th percentile duration calculated as 5 h 45 min, and the 90th percentile duration calculated as 33 h 45 min. The mean hydrograph duration at gauge 8006 (Spey) was 45 h 25 min, with the 10th percentile duration calculated as 17 h 14 min and the 90th percentile duration calculated as 75 h 51 min. There was variance in the hydrograph windows as a function of the variability in the magnitude of river discharge peaks at all gauges, and this method provided an indication of the hydrograph duration to support analysis of the likelihood of compound events occurring at each gauge. This method identified co-occurrence events with both extreme river discharge and S, which could result in flooding.

Fig. 3
figure 3

Schematic to show methods for calculating a DQ from idealised hydrographs 1 and 30 that exceeded the 50th percentile, b DTWL that exceeded 50th percentile and ce TQ-TWL and DQ+TWL. Pink shading indicates duration of Q exceeding 50th percentile, and orange shading indicates duration of TWL exceeding 50th percentile

Method 4: Top 500 Q

The largest 500 peaks in each (smoothed) Q record were identified to represent extreme river behaviour. A storm de-clustering algorithm was used to ensure top 500 Q peaks are independent of each other (Haigh et al. 2016). The storm de-clustering method discards all other smaller magnitude peaks around the Q peak that are within a storm length, based on the hydrograph window. For example, if the hydrograph window is 24 h, then smaller peaks 12 h either side of the Q peak is removed from the record. For these 500 extreme Q events, the maximum S value in the hydrograph window, from Method 3, was then identified. All S are included in this method, not just the extreme S above a threshold, to ensure as close to 500 co-occurrences are identified to build a dataset that is representative of the types of co-occurrences at each gauge, e.g. not just extreme Q and extreme S, but also extreme Q and less extreme S. This will give a correlation result that is representative of the dependence between Q and S given Q is extreme at each gauge. If no corresponding S was in the record due to data gaps, then the river peak is discarded from the analysis.

Method 5: Top 500 Skew

This method was similar to Method 4, but instead the largest 500 S values in the TWL record were first identified. For these 500 extreme skew events, the largest peak value of the (smoothed) Q record within the hydrograph window, about S, was identified. As with Method 4, this procedure generated a dataset that was representative of the dependence between S and Q, given that S is extreme. Methods 4 and 5 would likely identify many similar events, but not necessarily identical, which could lead to differing results.

Dependence and Extreme Co-Occurrence

Two metrics were then used to assess co-occurrences in extreme S and Q across 126 river discharge-tide gauge pairs. Firstly, dependence between S and Q was calculated, for each of the five data sets (Methods 1–5) described in “Does Daily Mean Discharge Data Capture Compound Events?” section, using Kendall’s rank correlation τ (Kendall 1938) which captures non-linear relationships and identifies the chance of the most extreme Q and S values coinciding (i.e. joint severity). Secondly, we calculate the mean number of compound events per storm season (1 June–31 May) to show, historically, how often Q and S are both extreme at each gauge. This metric, termed ‘annual mean compound events’, shows the average number of times extreme drivers occur together enough to co-occur, but does not account for the state of the tide or whether flooding occurred. Selecting events from June to May captures a full winter storm season (December–February) when largest peaks are most likely to occur. The number of extreme events per storm season is calculated for each storm season in the record and averaged over the 30-year record. For this second metric, an additional threshold was applied to each analysis to define extreme events as those where only Q and S values that both exceed the 95th percentile are included. This ensures that only the most extreme coincident Q-S pairs are considered for analysis. For Methods 4 and 5, only events where both drivers are extreme are included. The sensitivity of dependence and annual mean compound events to: (i) each dataset and (ii) each method, was analysed to identify locations that may be susceptible to compound events.

S–Q Dependence Using Daily Mean Values

The Kendall’s rank correlation τ and annual mean compound events calculated from Method 1 (daily means) were compared with results from Method 3 (POT Q, N = 50th) at each gauge to identify locations where a daily mean is and is not representative of river discharge behaviour. Method 1 was taken no further in the analysis as daily mean is not representative of river discharge variability.

Top 10 Locations

Of the 126 river gauges analysed, the top 10 in terms of strongest correlations and most annual mean compound events for Methods 2–5 were identified. The locations of the top 10 strongest correlations and annual mean compound events were compared between methods. This identifies if different methods produce strongest dependence and most annual mean compound events at the same or different gauges.

Strong Compound or No-Compound Events

Dependence and annual mean compound events were used to identify locations where compound events were most and least likely and could lead to flooding, based on whether Methods 2–5 calculated similar results. The following thresholds were applied:

  1. i)

    Strong compound events: Gauges where all methods calculated Kendall’s rank correlation τ over 0.15, or gauges where annual mean compound events exceeded 3.

  2. ii)

    Strong no-compound events: Gauges where all methods calculated Kendall’s rank correlation τ less than 0.05, or gauges where annual mean compound events less than 1.

  3. iii)

    Result varies: Gauges where some methods calculate strong compound events and some calculate strong no-compound events for dependence and annual mean compound events respectively.

The threshold to represent a ‘strong correlation’ is appropriate in the context of these results and the wider field of research. Similar thresholds are set in other research on compound flooding to suggest strong correlation. For example, Svensson and Jones (2004) set a threshold of 0.1 to represent a strong dependence between river discharge and storm surge in the UK. The range of correlation coefficients presented in this research also align with those presented in Hendry et al. (2019), which range from 0.1 to 0.35 on the west coast of UK, and 0.0 to 0.1 on the east coast. All correlation coefficients are statistically significant at 95% (as discussed in Supplementary Information Fig. S1). Altering the threshold would alter the results, and the gauges presented as being susceptible to compound events; however, this is an appropriate threshold within the scope of this research.

Analysis of Extreme Q and TWL Durations and Lag Times

Magnitude of peak river discharge (PQ), timing of peak river discharge (TQ), and timing of skew surge (TS) from Method 3 (POT Q, N = 50th) were used to calculate duration and timings of co-occurrences. The dependence between Q and S is important to understand the causal processes of compound events in estuaries. In practical terms, the potential for compound flooding to occur is largely influenced by the variability of astronomical tides. Therefore, the following analysis will also analyse TWL time-series, because the relative timing of PQ and PS in each 12-h tidal cycle is most likely to control if flooding will occur. The mean duration of the extreme river events, the mean duration of the extreme TWL event, the mean lag time between the events, and the mean duration of overlap between the events were calculated to better understand how event timings contribute to the likelihood of co-occurrence and flooding.

Q Duration (DQ)

The duration of peak river discharges cannot easily be calculated from each individual hydrograph event, as subsequent peaks can occur before the record has returned to base discharge levels and hence obscures the record. As described in “Does Daily Mean Discharge Data Capture Compound Events?” section, DQ was calculated using normalised, idealised gamma curves. The duration that each PQ exceeds the 50th percentile on the idealised hydrograph was calculated, which represents the duration Q exceeded the mean to represent severity of the storm. The duration of each storm was obtained by scaling up the duration of the paired gamma curve to the magnitude of peak river discharge, and then averaged across all events in the 30-year record to produce a representative DQ value for each gauge.

TWL Duration (DTWL)

The duration that TWL exceeded the 50th percentile was calculated and used as an indicator of the duration high-water levels exceed mean sea-level in each estuary (see Fig. 3b). The duration TWL exceeds the 50th percentile is likely to increase if a storm surge alters the height of tidal high water or an extreme river discharge occurs and is therefore representative of high-water levels. For each S event isolated in Method 3 (POT Q, N = 50th), the corresponding magnitude and timing of TWL, PTWL, and TTWL respectively, were identified. The time-series during the rising and falling limbs of the event were linearly interpolated from 15 to 1 min, and the duration that TWL exceeded the 50th percentile was calculated. These durations were then averaged across all events in the 30-year record for each gauge to produce a representative DTWL value for each gauge.

Lag Time (TQ-TWL)

The elapsed time between \({\mathrm T}_{\mathrm Q}\) and \({\mathrm T}_{\mathrm{TWL}}\) was calculated for each co-occurrence event in each record (see Eq. (1) and Fig. 3b). The mean lag time across all events in each 30-year record was then calculated to give a representative TQ-TWL value for each gauge.

$${\mathrm T}_{\mathrm Q-\mathrm{TWL}}={\mathrm T}_{\mathrm Q}-{\mathrm T}_{\mathrm{TWL}}$$
(1)

Duration of Overlap (DQ+TWL)

A series of equations were used to determine if an overlap of Q duration (DQ) and TWL duration (DTWL) was likely to occur. Equation (2) firstly identifies if an overlap is likely to occur, then Eqs. (3)–(8) calculate if there is a total or partial overlap of peaks and the duration of overlap. Equations (3)–(8) assume that Q and TWL are symmetric around the peak so they can be applied across all catchments.

Equation (2) is used to determine if there is or is not an overlap. If TQ-TWL exceeds DQ and DTWL then no overlap occurs:

$$\left(\frac{{\mathrm D}_{\mathrm Q}}2+\frac{{\mathrm D}_{\mathrm{TWL}}}2\right)<{\mathrm T}_{\mathrm Q-\mathrm{TWL}}$$
(2)

If an overlap occurs, then Eqs. (3) and (4) calculate if this is a total overlap. If DQ or DTWL is substantially longer than the other, the duration of overlap is assigned as the shorter of the two durations.

$$\frac{{\mathrm D}_{\mathrm{TWL}}}2>\left({\mathrm T}_{\mathrm Q-\mathrm{TWL}}+\frac{{\mathrm D}_{\mathrm Q}}2\right)$$
(3)
$$\frac{{\mathrm D}_{\mathrm Q}}2>\left({\mathrm T}_{\mathrm Q-\mathrm{TWL}}+\frac{{\mathrm D}_{\mathrm{TWL}}}2\right)$$
(4)

The following example calculates duration of total overlap where DQ is 18-h, DTWL is 7-h, and TQ-TWL is 4-h. DQ exceeds the right-hand side of the equation, indicating Q overlaps the full duration of DTWL creating an DQ+TWL of 7 h.

Equations (5) and (6) determine if a partial overlap between DQ and DTWL occurs, and Eq. (7) calculates what the duration of overlap should be assigned as.

$$\frac{{\mathrm D}_{\mathrm{TWL}}}2<\left({\mathrm T}_{\mathrm Q-\mathrm{TWL}}+\frac{{\mathrm D}_{\mathrm Q}}2\right)$$
(5)
$$\frac{{\mathrm D}_{\mathrm Q}}2<\left({\mathrm T}_{\mathrm Q-\mathrm{TWL}}+\frac{{\mathrm D}_{\mathrm{TWL}}}2\right)$$
(6)

Note than Eqs. (5) and (6) are similar to 3 and 4, but the duration on the left-hand side of the equation is shown as being smaller than the right-hand side. The overlap is calculated as shown in Eq. (7).

$$\frac{{\mathrm D}_{\mathrm Q}}2+\left(\left(\frac{{\mathrm D}_{\mathrm{TWL}}}2\right)-{\mathrm T}_{\mathrm Q-\mathrm{TWL}}\right)$$
(7)

As an example, in Eq. (8), DQ is 5 h, DTWL is 10 h, and TQ−TWL is 4 h:

$$\frac{5}{2}+ \left(\left(\frac{10}{2}\right)- 4\right)$$
(8)
$$2.5+ 1=3.5$$

These durations were then averaged across all events in the 30-year record for each gauge to produce a representative DQ+TWL value for each gauge.

Locations Susceptible to Potential Compound Events

Results from the above dependence, annual mean compound events, and duration analyses were combined to identify locations that have potential to flood. Locations were identified based on a series of thresholds that Kendall’s rank correlation τ and annual mean compound events must exceed at each gauge and an overlap of Q and TWL must occur. The results also highlight that when the thresholds are changed then different locations are identified as susceptible.

Seasonality of Compound Events

The final stage of the analysis identified seasonal trends in TWL and Q to understand when co-occurrences may happen at each gauge. The top 10 PQ and top 10 PTWL each year (1984–2013) were isolated, and then plotted to identify which months they occurred in (Fig. 4a and b). The monthly occurrence of top 10 Q and TWL are then added together (Fig. 4c), and the variance of the sum of monthly extreme occurrences calculated for each gauge. The variance was then plotted for each gauge, to identify where seasonality may cause more extreme Q and TWL to coincide. Figure 4 shows variance at gauge 66011 (Conwy) where most top 10 Q occurred in January and December and maximum tides occurred around the spring/autumn equinoxes. Gauge 91002 (Lochy) is shown in Supplementary Material Fig. S4, where maximum total water levels do not occur at the equinox but coincide with maximum Q in winter months. Skew surges at Lochy make up a larger proportion of total water level as tidal range (mean high water spring—mean low water spring) is small 2.95 m, compared with other tide gauges in UK which see a tidal range over 4 m.

Fig. 4
figure 4

Seasonality of a top 10 annual Q; b top 10 annual TWL; and c sum of Q occurrences at 66011 (Conwy) and TWL occurrences at Llandudno tide gauge

Results

The results presented here identify spatial patterns throughout British estuaries in the dependence between extreme skew surge and extreme river discharge behaviour and the number of seasonal co-occurrences of these events. This has enabled us to identify estuaries particularly susceptible to compound events. Results are based on 126 river gauges and 27 tide gauges, covering 30 years (1984–2013) at 15-min temporal resolution.

Does Daily Mean Discharge Data Capture Compound Events?

The first objective of the research was to establish if daily mean river discharge is of sufficient temporal resolution to accurately represent river behaviour for compound events analyses across Britain. Most previous analyses (e.g. Hendry et al. 2019) used daily mean river discharge data. Sixty-three gauges have a hydrograph window, as described in “Does Daily Mean Discharge Data Capture Compound Events?” section, of less than 12 h, and 26 gauges have a refined hydrograph window between 12 and 24 h (Fig. 5a). Out of the 126 gauges, 89 (68%) have a hydrograph window less than 24 h, therefore daily mean data cannot accurately resolve the magnitude of peak Q (PQ) as they will be smoothed out. Method 1 will underestimate the number of annual mean compound events. Figure 5b shows that river gauges with a hydrograph window less than 24 h occur throughout Britain. The differences in Kendall’s rank correlation τ and annual mean compound events between Methods 1 and 3 are shown in Supplementary Material Fig. S5. Method 3 produces stronger correlation and more annual mean compound events than Method 1 in gauges across Britain. Method 1, which uses daily mean Q, is therefore not used further in this analysis as it will not accurately capture river hydrograph behaviour.

Fig. 5
figure 5

a Frequency of hydrograph window. b Hydrograph window for each river gauge. c Relationship between hydrograph window and catchment area

Longest DQ occurs at gauges 54,057 (Severn, 125 h) and 6007 (Ness, 51.6 h) (Fig. 5c). Gauge 6007 (Ness, 51.6 h) has a smaller catchment area but longer DQ because the rising limb of most of the hydrographs is not steep (the hydrograph is damped by influence of Loch Ness and storage of water as snow in Scottish catchments), but the magnitudes are large therefore shallower and longer duration idealised hydrographs are matched to this gauge. Outliers include gauge 39,001 (Thames, 77.8 h) which has larger catchment area but smaller durations due to fewer peaks selected in Method 3. A combination of high baseflow in larger estuaries and ultrasonic gauges used at these locations, which are sensitive to capture high frequency variability, means the relative thresholds are high so fewer peaks exceed this. Gauge 55023 (Wye, 93.87 h) has a longer DQ for the size of the catchment.

Spatial Variation in Dependence and Annual Mean Compound Events

The second objective of this research was to identify the best methods to identify joint occurrence of high storm surges and high river discharges around the coast of UK. Compound event dependence (Kendall rank correlation, τ) and annual mean compound events maps, based on analysis Methods 2–5, are shown in Fig. 6. τ is calculated from the full record of PQ and PS, and annual mean compound events represent a mean value across 30 years, with a range of annual values around the mean, to represent how often the drivers are both extreme from June to May each year.

Fig. 6
figure 6

Kendall rank correlation τ ad and annual mean compound events eh for each of the four methods to select extreme compound events

Each of the four methods produced a different spatial pattern of dependence and annual mean compound events. Method 2 (daily maximum Q and S) shows a clear east–west split throughout Britain in τ, with strongest τ up to 0.27 at gauge 79002 (Nith) on the west coast, indicating greater chance of the most extremes coinciding. The highest number of annual mean compound events for Method 2 occur on the west coast of Scotland, with 9.58 at gauge 89005 (Lochy), 8.79 at gauge 89007 (Abhainn a’ Bhealaich), and 8.48 at gauge 84011 (Gryfe). There were few occurrences per season in SW-England, despite the strong correlations, and on the east coast. Two gauges record no annual mean compound events when S exceeds 95th at gauge 54057 (Severn) and 39001 (Thames).

Method 3 (POT Q) also shows an east–west split in τ, with strongest τ up to 0.32 occurring consistently on the west coast of Britain. The strongest correlation is recorded at gauge 76007 (Eden), where there is a greater chance of extreme Q and extreme S coinciding. Most annual mean compound events occur in NW-England and W-Scotland; up to 8.89 extreme S co-occur with peak Q at gauge 72004 (Lune), and there is variance around this mean value with a minimum of 2 and a maximum of 13 occurrences per season. At gauge 8006 (Spey), 1.82 seasonal cooccurrences are recorded, with a minimum of 0 and a maximum of 4. Fewer annual mean compound events are recorded in SW-England, despite the high correlations. No annual mean compound events are recorded at gauge 54057 (Severn) and 39001 (Thames). These gauges have a very high baseflow, and no PQ and PS exceeds the 3-month moving mean + 95th percentile threshold within the hydrograph window.

Method 4 (top 500 PQ) shows weaker correlations in gauges across north of Britain, apart from gauge 70,004 (Yarrow, τ = 0.35). Stronger correlations occurred in SW-England, but most gauges showed weak correlations. At some gauges, e.g. 5480 (Roding) and 52007 (Parrett), not all top 500 PQ have an associated PS value within the hydrograph window. Gauges with a hydrograph window less than 12 h show that PQ are not always going to be associated with a skew surge value that occurs every 12 h at the time of tidal high water. The top 500 PQ will be identified, but only 100 PQ may have a PS. Other gauges, including 22001 (Coquet) has 490 occurrences between PQ and PS but is weakly correlated (t = 0.04) as most PS are non-extreme. A greater number of annual mean compound events are seen in the north of Britain, predominantly on the west coast but also some gauges show higher occurrences on the east coast, e.g. gauge 27009 (Ouse) with 8.03 occurrences. Higher occurrences are recorded at gauge 54057 (Severn) and 39001 (Thames) than with other methods; this method is not dependent on Q exceeding a higher threshold.

Method 5 (top 500 PS) also shows stronger correlations in SW- and NW-England, similar to Methods 2 and 3. All gauges record between 475 and 500 co-occurrences between PQ and PS with this method to calculate correlations. This is because Q is at 15-min resolution, and it is possible to identify maximum PQ in the hydrograph window around PS. Higher annual mean compound events occur in W-Britain, and gauge 39001 (Thames) does not generate any compound events because of no corresponding Q exceeding the 95th percentile.

Statistical Agreement Between Co-occurrence Methods

Here we demonstrate the differences that occur due to the different analysis Methods 2–5. We plotted the dependence τ from each method against each other in Fig. 7, with RMSE, R2, and p value calculated between the two datasets. The plots are ranked in the order of agreement (R2) (a–f: weakest to strongest). Method 5 (top 500 PQ) and Method 4 (top 500 PS) show the weakest agreement with dependence produced from Method 2 (daily max.) to Method 3 (POT Q, N = 50th). This is most likely due to the larger datasets used in Methods 4 and 5 which generate weaker correlations overall. Method 2 (daily max.) and Method 3 (POT, N = 50th) show best agreement and smallest error; these methods identify similar PQ in the records but different methods for selecting PS.

Fig. 7
figure 7

Kendall rank correlation τ calculated from each method plotted against each other, ranked in the order of R2

Figure 8 shows annual mean compound events for each method plotted against each other and ranked in the order of R2. Method 4 (top 500 PQ) generates weak R2 and large RMSE when compared with other methods. Method 2 (daily max) and Method 3 (POTQ) show reasonable agreement; there may be instances where there are two peaks in a hydrograph in 1 day, and a daily maximum is only able to identify one of them. Method 5 (top 500 PS) shows good agreement with Method 2 (daily max) and Method 3 (POT Q). Identifying annual mean compound events in Method 5 is also dependent on identifying PQ above the 3-month moving mean + 95th percentile threshold, which therefore generates more similar results than Method 4 (top 500 PQ).

Fig. 8
figure 8

Annual mean compound events calculated from each method plotted against each other, ranked in the order of R2

Top 10 and Strong Compound Events

The third objective of this research was to identify estuaries that are susceptible to compound events. We classified each gauge based on whether there was an agreement in τ and annual mean compound events between the different analysis methods (Figs. 9 and S7). If high correlation and frequent annual mean compound events were calculated across all four methods, then this indicates strong compound event occurrence. Weak correlation and few annual mean compound events between all methods indicate low occurrences of compound events. If one method did not produce a high or low result, then it is classified as ‘result varies’. The colour of each circle in Fig. 10 identifies gauges where the top 10 strongest τ or highest annual mean compound events was the same across four, three, two, or just one method.

Fig. 9
figure 9

Gauges across Britain which show (i) strong compound (filled red circles); (ii) no-compound flooding (grey circles); or (iii) varied results (white circles) and appear in the top 10 across (i) 4 methods (dark red outline); (ii) 3 methods (red outline); (iii) 2 methods (orange outline); or (iv) 1 method (yellow outline) for a Kendall rank correlation τ and b annual mean compound events

Fig. 10
figure 10

Mean duration of a Q (DQ) and b TWL (DTWL) exceeding 50th percentile; c mean lag time (TQ-TWL); and d mean duration of overlap (DQ+TWL) from Method 3 (POT, N = 50th). Refer to schematic in Fig. 3 for details on how each parameter is calculated. Note that the colour scale indicates shorter durations in red

Figure 9a indicates that the methods used here identified strong compound events in gauges in W-Britain, with a higher proportion of red circles indicating high correlations are calculated here across the four methods. Strong no-compound events occur in NE-Britain, with 10 gauges showing consistently weak τ. Six gauges appear in the top 10 strongest correlations for three out of the four methods: 4001 (Conon), 6007 (Ness), 54,057 (Severn), 55,023 (Wye), 76,007 (Eden), and 79,002 (Nith). Gauge 54,057 (Severn) has strong correlations but low annual mean compound events. In the Severn, Methods 2–5 each identify over 400 instances when PQ and PS exceed 50th percentile within a longer hydrograph window (125 h), and these events are well correlated. However, the additional threshold for calculating annual mean compound events means that there are substantially fewer times that PQ and PS exceed the 95th percentile and extreme events are not captured. Gauge 37,010 (Blackwater) in SE-England shows correlation occurs in the top 10 for Method 4 (top 500 Q). This is because out of 500 PQ, only 91 occur with an extreme skew surge in the 17.9-h hydrograph window, and those 91 are strongly correlated.

Figure 9b also shows that strong compound events (red-filled circles) were more widely spread across Britain when considering annual mean compound events; however, gauges which appear in the top 10 highest number of annual mean compound events in two or more methods are all on the west coast of Britain. Two gauges have a top 10 result for annual mean compound events which appears across all four methods (dark red); gauge 72,004 (Lune) and 89,003 (Orchy). Eight gauges show strong compound events on the east coast of England and Scotland, and three gauge have a top 10 result for annual mean compound events. These gauges have a hydrograph window that exceeds 20 h, therefore more PS are likely to be identified. Hence these could be combination events rather than compound (i.e. occurring at the same time but not linked to the same storm).

Duration of Drivers of Compound Events

To further explore the third objective of this research, to identify estuaries that are susceptible to compound events, data from Method 3 (POT Q) were used to quantify the durations of Q (DQ) and TWL (DTWL), the lag time between the peaks of these events (TQ-TWL), and the duration of overlap (DQ+TWL) (Fig. 10). Method 3 was selected for further analysis as the 15-min frequency Q more accurately captures PQ and allows for more detailed analysis of hydrograph duration. Refer to schematic in Fig. 3 for details on how each parameter is calculated. Note that the colour scale indicates shorter durations in red.

A shorter DQ indicates a flashier hydrograph (Fig. 10a): There are no clear spatial trends because DQ is more closely linked to catchment area and flow regime, with a longer Q duration occurring in larger catchments (Fig. 5c). Further to this, a larger catchment or greater magnitude of peak Q generates greater variance around mean DQ. DQ at gauge 37,001 (Roding, max Q: 75 m3/s, area: 303 km2) is 1 h 18 min, and the 10th percentile is 20 min and 90th percentile is 3 h. Whereas DQ at gauge 55,023 (Wye, max Q: 805 m3/s, area: 4010 km2) is 44 h 17 min, and the 10th percentile is 20 h 18 min and 90th percentile is 76 h 6 min.

Figure 10b shows duration of TWL exceeding 50th percentile (DTWL) for the nearest tide to TQ. The durations of TWL are linked to local tidal dynamics and characteristics of the tide captured by each tide gauge paired with the river gauges. DTWL in the Severn Estuary has a short duration where there is a strongly asymmetrical tide which causes rapidly rising flood tide, e.g. 54,032 (Severn, 6.4 h). DTWL varies between gauges using the same tide gauge based on the river events selected in the record and the specific high tide linked to the river event. Longer durations occur in NE-Scotland, e.g. gauge 6007 (Ness, 7.7 h), and NW-Scotland, e.g. 84,013 (Clyde, 7.5 h). DTWL shows some variance around the mean at each gauge, with some greatest variance at Clyde (paired with Milport tide gauge), where the 10th percentile is 6 h 17 min and 90th percentile is 8 h 45 min. Smaller variance is shown at Conwy (paired with Llandudno tide gauge) where the 10th percentile is 6 h 8 min and 90th percentile is 7 h 5 min.

The lag time (TQ-TWL) varies spatially from less than 1 h, e.g. 37,006 (Can, 40 min) to up to 36 h at 54,032 (Severn) (Fig. 10c). A shorter lag time indicates that Q and S occur in quick succession of each other, whereas a longer lag time indicates that co-occurrence between Q and S is less likely. TQ-TWL at each gauge will show variance within the duration of the hydrograph window, which could be up to 50 h in some larger catchments. Figure 10d shows gauges where an overlap between Q and TWL is likely to occur, and the estimated duration of the overlap. This is calculated based on TQ and TTWL and the lag times; for example, a shorter lag time and longer Q and TWL durations will cause an overlap. A longer duration of overlap indicates peak river discharge and total water level occur together over a longer period, and could contribute to flooding. Longer durations of overlap occur in NW-England and W-Scotland, whereas shorter durations of overlap occur in SW-England and SE-England.

This analysis was also completed using durations that Q and TWL exceed the 95th percentile (results not shown); however, these durations were too short with the associated lag times, and no gauges showed overlap. This indicates it is very unlikely that the most extreme flows and water levels would co-occur.

Locations Most Susceptible to Compound Events

The following analysis continues to use results from Method 3 (POT Q) and combines results from the previous sections to identify locations that are most susceptible to compound events based on the Kendall’s rank correlation τ, annual mean compound events, and potential for Q and TWL to overlap based on duration and lag time.

Figure 11 shows locations where τ > 0.2, annual mean compound events > 3, and average Q and S duration likely overlap (e.g. if average Q and S duration both exceed the lag time then co-occurrence is likely). Twenty-six gauges, listed in Table 1, are most susceptible to compound events based on these criteria. The potential for compound flooding to occur is more likely with a higher tau as the most extreme peak flows are more likely to co-occur with the most extreme surges. Most gauges are located on the west coast, with just one exception in NE-England. Gauge 72004 (Lune) has the highest number of annual mean compound events (8.9), and these have also been split into the mean number of winter (December, January, February) and summer (June, July, August) occurrences. Gauge 74001 (Duddon) has the highest number of summer occurrences indicating that compound events could happen at any time in the year.

Fig. 11
figure 11

Gauges that are susceptible to co-occurrence events based on criteria of Kendall’s rank correlation τ > 0.2, annual mean compound events > 3, and overlap is likely

Table 1 Gauges susceptible to compound events based on criteria of Kendall’s rank correlation τ > 0.2, annual mean compound events > 3, and overlap is likely

The following result shows that different gauges are found to be potentially susceptible to compound events when the criteria are changed. The τ identifies locations that are potentially likely to experience flooding due to compound events, as a higher τ indicates most extreme peak flows are more likely to occur with the most extreme surges. Figure 12 shows gauges that are susceptible to compound events focus when the τ is not used as a criteria, where annual mean compound events > 3 and Q and S likely overlap.

Fig. 12
figure 12

Gauges that are susceptible to co-occurrence events based on criteria of annual mean compound events > 3, and overlap is likely

Table 2 lists 46 gauges that show high annual mean compound events and potential for overlap between Q and S. These gauges show a very high number of annual mean compound events (Table 2); the datasets recording peak river discharges and corresponding extreme S are larger indicating there is more potential for flooding, but there are also more events that will not cause compound events as shown by weak τ. Gauges susceptible to compound events are in NW- and N-Britain, with the most gauges clustered in NW-England and W-Scotland. Gauges on the east coast may represent higher occurrences of combination events.

Table 2 Gauges susceptible to compound events based on criteria of annual mean compound events > 3, and overlap is likely

Seasonality of Q and TWL Contributing to Likelihood of Co-occurrence

The fourth objective of the research was to identify seasonal trends in sea-level and river discharge which could increase the potential for compound flooding at certain times of year. The variance of the sum of monthly occurrences of most extreme TWL and Q is shown in Fig. 13. Gauges in W-Scotland, including 91002 (Lochy), 84013 (Clyde), and 4001 (Conon) shows highest variance. The gauges in W-Scotland which show high variance all correspond to the Millport tide gauge. As seen in Fig. S5 in Supplementary Material, the top 10 TWL and Q both most frequently occur in January and December. The weather systems in W-Scotland appear very seasonal, and frequently occur in January and December to cause both large storm surges and peak river discharges. This is reiterated in previous results which show NW- and W-Scotland see high annual mean compound events. Storm surges also contribute a larger portion of TWL here, due to the smaller tidal range. There is more potential for compound events to happen in this region in winter months.

Fig. 13
figure 13

Variance of the sum of monthly frequencies of most extreme observed peak river discharges and most extreme observed total water level (surge + tide)

SW-England and S-Wales, which have high Kendall’s rank correlation τ in Method 2 (daily max.) and Method 3 (POT Q), have low variance as the top 10 TWL are controlled by the largest tides occurring at the equinox in spring and autumn. Weaker variance in SW-England and S-Wales indicates less seasonality and monthly variation. Seasonality has less influence on when extreme Q and TWL happen through the year, and some sites may show constantly high occurrences of extreme Q and TWL. A stronger Kendall’s rank correlation τ in SW-England and W-Wales indicates that when a storm does occur, it is more likely it could increase the potential for flooding.

Discussion

Compound events in estuaries arise as a result of high sea-levels caused by storm surges and high astronomical tides occurring at the same time as, or in close succession to, high river discharge (Svensson and Jones 2004; Ward et al. 2018). The combination of extreme sea- and river-levels can exceed critical thresholds and pose a worst-case flooding hazard to coastal communities, infrastructure, and ecosystems, compared with situations where the two drivers occur separately. It is important to understand where and when co-occurrence between sea- and river-level has happened in the past to understand how trends may alter under future climate change and land-use change scenarios (Robins et al. 2016). In this paper, we have shown that extreme river flow durations in most catchments across Britain are indeed < 24 h and therefore require a new analysis based on sub-daily resolution data. We have shown that different methods of data selection (e.g. hourly vs daily Q, or POT Q) and different metrics for presentation of results (correlation or annual mean compound events), lead to different findings. The results presented here show that daily mean Q data generally do not capture river variability. Daily mean Q data smooths the record, and the magnitude of PQ (magnitude of peak discharge values) are underestimated and not identified in catchments with flashy behaviour or a DQ shorter than 24 h. Daily mean Q data may be suitable for identifying PQ in some larger estuaries (e.g. Severn, 54,057; Thames, 39,001), which have Q durations longer than 24 h but should be avoided in estuaries where hydrograph duration is less than 24 h.

Method 2, which used daily maximum Q, has been shown to produce results which agree with Method 3. This dataset could be used in all estuaries to capture a basic understanding of the dependence between the magnitude of PQ and PS. However daily maximum Q does not capture the specific timing of PQ and misses multiple PQ if more than one peak occurred on the same day.

Method 3, which used 15-min frequency Q, is the highest temporal resolution available and should be utilised in analysis where possible, in particular where estuary hydrograph duration is less than 24 h. Method 3 (POT Q, N = 50th) is designed to ensure at least 100 PQ are selected above high thresholds to represent river behaviour. This value of 100 PQ is set specifically to ensure that enough peaks are selected in gauges where ultrasonic gauges are used, which cause high baseflow and high thresholds, e.g. gauges 54057 (Severn) and 39001 (Thames), or where artificial influences alter the flow regime, e.g. gauge 37010 (Blackwater). The method would be strengthened by applying different thresholds for flooding for different locations in the POT Q analysis, to understand the likelihood of compound flooding occurring as a result of compound events. In practical terms, it would be useful to understand if the analysis to identify extreme peaks in Q and S and potential for compound events does lead to flooding.

The 15-min frequency Q more accurately captures PQ (magnitude of peak discharge values) and TQ (timing of peak discharge values) to allow for analysis of hydrograph window, DQ (duration hydrograph exceeds a threshold), TQ-TWL (lag time), and DQ+TWL (duration of overlap). The novel hydrograph window method is used to ensure only pairs of high discharge and skew surge or total water level that will overlap are considered in the analysis. The method could be developed to use the duration calculated for each individual PQ to identify the nearest skew surge or total water level. DQ uses the same gamma curve fitting procedure to calculate the duration that each hydrograph exceeds the 50th percentile on the idealised hydrograph. The threshold could be raised or lowered to give a representative duration that flow exceeds a more or less extreme threshold. An alternative approach is to calculate duration from the observed record, but flow does not always return to low values due to clustered events so this would not use all PQ in record. This would skew the DQ as this value would be calculated from only events which do return to zero, which can be challenging to isolate and identify in flashy catchments.

Method 3 assigned PQ to S (skew surge) and TWL (total water level). Each variable serves a purpose to better understand the drivers of compound events, and potential for compound flooding, to occur in estuaries. The skew surge was used to understand more about the chance of the most extremes coinciding and how often the drivers were both extreme. The same skew surge value may be assigned to multiple peak river discharges if these are clustered close together and selected skew surge may occur before or after PQ. Only the most extreme PQ and PS were selected, which may skew the correlation results. Pairing PQ and total water level identified more about the potential for compound flooding to have actually occurred, as the magnitude of the tide is what will have likely caused water levels to breach defences (Ward et al. 2018; Lucey and Gallien 2021). Further research could quantify the value of using skew surge compared with residual surge when identifying estuaries susceptible to compound events. Skew surge and total water level are both variables that further develop understanding of co-occurrence events in the UK, where tides are large, and methods could be applied to estuaries worldwide to further understand causal processes of compound events.

In Method 3, a clear understanding of Q and S durations, and duration of likely overlap have been used to identify gauges susceptible to compound events. The hazard assessment metrics presented here (e.g. DQ, TQ-TWL) can be used by local authorities and national agencies to identify catchments where coupled tide-surge-river flooding models can be applied to support long-term hazard planning. Accurate hazard metrics and coupled models can inform early warning systems and evacuations; action can be taken in these estuaries if a peak in river discharge due to heavy rainfall is likely to occur at a similar time of a large forecast storm surge. There is good statistical agreement in τ and annual mean compound events between Methods 2 and 3 indicating daily maximum Q and 15-min frequency Q may both be suitable to identify PQ and identify correlations and annual mean compound events, as they both capture the magnitude of extreme river discharge.

Method 4 did not always isolate a corresponding skew surge, and so the top 500 PQ were not always included in correlations. Method 5 selected maximum Q within the refined hydrograph window about S, did not necessarily capture peaks in Q time series. Each method can capture co-occurrences, and data availability and aims of the research will determine which method could be used for future analysis. This research has showed that Method 2 (daily maximum Q) and Method 3 (15-min data) are the preferred dataset to use to capture river behaviour to identify likely co-occurrences.

Like earlier studies (i.e. Svensson and Jones 2004; Hendry et al. 2019), Methods 2, 3, and 5 identify an E-W split in Kendall’s rank correlation τ and stronger correlations are seen in W-Britain. Methods 2, 3, 4, and 5 identify more annual mean compound events in NW-Britain; however, Method 4 also identifies more in E-England. Method 4 (top 500 Q) does not always generate a large dataset between Q and S, as the hydrograph window may not be sufficiently long enough to capture PS. Most events from Method 4 are not co-occurrences as correlations are weak overall, apart from exceptions in SW-England; PQ does not necessarily co-occur with an extreme PS in the hydrograph window. Method 5 (top 500 S) more consistently generates a large dataset between Q and S in a hydrograph window to give an accurate representation of river- and sea-level behaviour, and more closely resembles results from Methods 2 and 3. Events identified in Method 5, when co-occurrences are identified dependent on PS, are more likely to be co-occurrences. Storms which generate PS will likely also generate PQ.

We have identified gauges across W-Britain that are susceptible to compound events, and spatial trends in Kendall’s rank correlation τ and annual mean compound events highlight drivers respond differently on the west coast. Stronger Kendall’s rank correlation τ are evident in SW-England, but fewer annual mean compound events are seen in this region. This indicates that when there is a low-pressure system which generates a storm in SW-England, it will normally lead to peak river flows (Svensson and Jones 2004). The largest surges co-occur with the largest peak flows, but surge and peak flow are rarely simultaneously above their 95% threshold. This highlights that interpreting dependence based on correlation alone can be misleading; stronger correlation indicates the most extreme sea-levels and river discharges can occur concurrently, but it does not necessarily happen very often. High river discharges can occur more often with high skew surges in catchments with a smaller area, such as those in Devon and Cornwall in SW-England and S-Wales (Hendry et al. 2019). Tidal range in SW-England, notably the Severn Estuary, is large (12.27 m at Avonmouth) which can lead to extreme coastal water levels (Lyddon et al. 2019) and contribute to more severe compound flooding when it does occur. During very high tides, large estuaries can also become tidally locked due to the backwater effect, and the level of the incoming high tide stops the river water flowing out to the sea, e.g. Dee Estuary (Cai et al. 2014; Environmental Agency 2016). Therefore, it is particularly important to understand additional hazard assessment metrics associated with total water level (e.g. DQ, PTWL, TQ-TWL) in larger estuaries, and when these lead to compound flooding.

Strong Kendall’s rank correlation τ are generated in NW-England and W-Scotland, with some gauges also seeing high annual mean compound events. When a low-pressure system or convective storm passes through NW-England or W-Scotland, then compound events do occur. It is more useful for practical flooding purposes to use a ‘annual mean compound events’ metric, to understand how many times Q and S drivers have both been extreme and co-occurred. Catchments respond quickly to rainfall and peak river discharges are likely to occur on the same day as extreme skew surges (Svensson and Jones 2004). Fewer annual mean compound events were identified on the east coast of Britain, which has previously been identified due to East coast surge and precipitation events being driven by different storm patterns (Svensson and Jones 2002; Zong and Tooley 2003; Camus et al. 2022). Calculating annual mean compound events is a valuable and practical metric as this can identify regions where there is a higher potential for compound flooding events as they have already occurred.

Further to this, the results confirm that compound events are more likely to occur in winter across all study sites, with most notable variance in monthly frequency of extreme Q and TWL in NW-England and W-Scotland. Top 10 annual Q and TWL values coincide in December and January, indicating that extreme sea-levels and river discharges are more likely to co-occur in winter. Gauges with a smaller tidal range, e.g. Milport on the Clyde, W-Scotland where tidal range is 2.95 m, may experience compound flooding as Q and S make up a greater proportion of total water level. There is also a high likelihood of the most extreme sea-levels and river discharges co-occurring on the Cumbria coast, NW-England, in winter due to high variance in monthly extremes. In this region, heavy rain throughout the year (due to convective thunderstorms) generates large river discharges, short transmission times (due to steep catchments and impermeable geology), and a larger tidal range (7.42 m at Workington and 8.49 m at Heysham) with high surge values increase the magnitude of co-occurrence events. The River Lune (gauge 72,004) and Eden (gauge 76,007) are highlighted as particularly susceptible to co-occurrence events, which is important to note for flood managers in Lancaster and Carlisle. These are large cities in the lower reaches of the rivers, which could be susceptible to compound events and potential flooding and have previously been subjected to severe flooding, e.g. during Storm Desmond (Martin 2015; Mortimer 2015).

The results presented here identify locations that are susceptible to co-occurrence events, with the majority of estuaries vulnerable to compound flooding occurring west Britain. This research can be interpreted alongside recorded instances of flooding to explore exact Q, S, and TWL conditions lead to compound flooding. This analysis would help to explore in more detail the exact timing and magnitude of drivers which lead to flooding, and to set more accurate thresholds used in Method 3 (POT Q). Additional research to capture the state of the tide when skew surge and river discharge co-occur would also help to understand the contribution of this driver to compound flooding. Future work could also consider site-specific equations to generate a representative hydrograph and Q duration which accounts for asymmetry around the peaks in the hydrograph or total water level time-series. Further to this, temporal trends in Q and S records can be explored to see if it is likely that the maps presented here will change in the future. Changes in sea-level and storm characteristics associated with climate change (Seneviratne et al. 2012) could alter the timing and magnitude of Q and S peaks, and make co-occurrence events more likely at gauges where they are currently not susceptible (Zscheischler et al. 2018), or less likely at other locations. Future research could build on analysis by Harrison et al. (2021) which investigated changes in the compound flood hazard due to projected sea-level rise and changes in surge and fluvial discharge for the Dyfi and Humber, U.K. Future projections of sea-level and river discharge and changes in evaporation, soil moisture, and groundwater contribution could be utilised in further study to identify how timing and duration of peaks may alter the likelihood of compound flooding events across the UK. Modelling studies are needed to explore impact of intra-estuary processes and inundation because of co-occurrence events.

Conclusion

This is the first time that sub-daily analysis of sea-level and river discharge has been extended across differing estuaries and catchment types across the UK. Thirty years of historical sea-level and river flow data, at 15-min instantaneous output frequency, is analysed for 126 estuaries across to Britain to identify where sub-daily analysis is required and estuaries susceptible to compound events. Different methods of data selection and identification of peak river discharge events generates different results. The results confirm that daily mean data is not sufficient to capture high frequency variability in river discharge events in 68% of gauges analysed. Smaller catchments have a river duration less than 24 h; therefore, daily mean data does not accurately capture the magnitude and timing of peak discharge events. Daily maximum and 15-min instantaneous output frequency river discharge data used in a peak over threshold method is sufficient to capture the magnitude of peak discharge events. However daily maximum river discharge is not high enough temporal resolution to analyse the duration of peak river discharge and total sea-level, lag time, or duration of overlap between the two drivers.

All methods show greater dependence between skew surge and river discharge on the west coast of Britain, with strongest Kendall’s rank correlation τ shown in SW- and NW-England. Annual mean compound events are highest in NW-England. The most extremes in sea-level and river discharge at gauges 72,004 (Lune), 89,003 (Orchy), and 82,003 (Stinchar) occur at the same time and are well correlated, which increases the potential for compound flooding to happen. Interpreting dependence based on correlation alone can be misleading, and it is more useful to understand when compound events have happened in the past. Further to this, analysis of the relative timings of total water level (tide + storm surge) and river discharges shows that Q and S drivers can overlap for up to 4.5 h, increasing the likelihood of co-occurrence and duration for potential flooding to occur. Estuaries that are most susceptible to compound events based on the Kendall’s rank correlation τ, annual mean compound events, and potential for Q and TWL to overlap based on duration and lag time are identified. Different thresholds identify different gauges as susceptible to compound events, but overall help to identify which estuaries should be the focus of future research and hazard mitigation strategies to minimise the consequences of compound events and potential flooding impacts. Future research should consider how changing sea-levels and storm characteristics might alter the magnitude, duration, and lag times of the drivers of compound flooding events to alter which estuaries are susceptible to this hazard.