The co-incidence of storm surges and extreme discharges within the Rhine–Meuse Delta

The Netherlands is a low-lying coastal area and therefore threatened by both extreme river discharges from the Meuse and Rhine rivers and storm surges along the North Sea coastline. To date, in most flood risk analyses these two hazardous phenomena are considered independent. However, if there were a dependence between high sea water levels and extreme discharges this might result in higher design water levels, which might consequently have implications for flood protection policy in the Netherlands. In this study we explore the relation between high sea water levels at Hoek van Holland and high river discharges at Lobith. Different from previous studies, we use physical models forced by the same atmospheric forcing leading to concomitant and consistent time series of storm surge conditions and river discharge. These time series were generated for present day conditions as well as future climate projections and analysed for dependence within the upper tails of their distribution. In this study, dependence between the discharge at Lobith and storm surge at Hoek van Holland was found, and the dependence was highest for a lag of six days between the two processes. As no significant dependence of the threats was found for cases without time lag, there is no need for considering dependence in flood protection and policy making. Although future climate change is expected to lead to more extreme conditions in river discharges, we cannot conclude from this study that it will change the magnitude of the dependence for extreme conditions.


Introduction
Flooding is one of the most frequently occurring and damaging natural hazard phenomena. Flood losses, along with other impacts, can be caused by different flood processes such as coastal, river, pluvial and groundwater flooding. To protect areas against flooding, flood protection measures such as levees, reservoirs and storm surge barriers have been built. Within the Netherlands this has led to flood protection for coastal and river flooding for events up to return periods of 10 000 years and 1250 years respectively.
The probabilistic design and safety assessment of flood protection along the Rhine-Meuse delta in The Netherlands accounts for the combined occurrence of storm surge in the North Sea, as well as high river discharges in the Rhine and Meuse Rivers. Nevertheless they are assumed to be statistically independent, which may result in an underestimation of the actual risk. Especially with respect to future climate scenarios, this statistical dependence may become an increasingly relevant issue, for instance in the Dutch Delta Programme, which has the objective to define a strategy that will 'protect the Netherlands from flooding and ensure adequate supply of freshwater for generations ahead' (The 2011Delta Programme 2010. To support decisions on a strategy for flood defences, changes in storm and precipitation regimes need to be accounted for, including potential changes in statistical dependence.
A recent, relatively small scale, event with coinciding pluvial flooding and storm surge conditions in the province of Groningen (The Netherlands) (van den Hurk et al 2014), has triggered a debate on whether storm surge conditions and high discharges in the Netherlands are dependent and whether this may lead to higher flood hazards. This is of particular interest in the Rhine Delta area, the economic centre of the Netherlands (including the main port of Rotterdam, see figure 1), where flooding can potentially have disastrous consequences. The probability of co-incidence of such events across the Netherlands was first studied by Van den Brink et al (2005) and Kew et al (2013). These studies mainly relied on meteorological data: Van den Brink et al (2005) used ECMWF weather forecasts (Molteni et al 1996), while Kew et al (2013) focused on possible changes in co-incidence probability in future climate and used the ESSENCE 17member global climate model ensemble (Sterl et al 2008). The study with the ESSENCE ensemble reveals is that there is a significant relationship between rainfall events over the Rhine basin and NNW wind conditions that generally lead to storm surges. Similar results were found recently, based on analysis of observations of storm surge levels and precipitation on stations along the coastline of Australia (Zheng et al 2013(Zheng et al , 2014.
Different from earlier work, we study in this letter, the conditional probabilities of occurrence of high discharges in the Rhine and high storm surge levels along the Dutch coast with a hydrological and hydraulic model cascade. The rationale to use models rather than observations is that we also investigate the effect of possible future climate trends in the co-incidence of storm surge extremes and high discharges. The model cascade estimates time series of flows in the Rhine River, as well as water levels at the North Sea coastline. Hence we estimate the real hazard processes rather than meteorological proxies as done in earlier studies (Van den Brink et al 2005, Kew et al 2013. By doing so we include memory effects of hydrological processes, travel times, attenuation and local effects, which are not accounted for if using meteorological fields only. We study these probabilities under present-day conditions as well as for a set of climate change scenarios. Due to the long run times of the model cascade it was only possible to do a calculation for one climate projection member per scenario, so only 30 years per case. Therefore we study 90% quantiles rather than 99% quantiles like Kew et al (2013). However if an increasing dependence is visible for 90% quantiles it is likely that this will also be the case for 99th and higher quantiles, considering the behaviour of most models for bivariate variables (Joe 1997, Frees and Valdez 1998, Nelsen 1999, Diermanse and Geerse 2012).

Model set up
To estimate the joint occurrence of high flows and extreme storm surge levels, we use the Climate Knowledge Facility (CKF). CKF consists of a set of hydrological, hydrodynamic, water quality, control and coastal models that allow the computation of combinations of hazardous phenomena or water system interactions between e.g. coast and river. The CKF system is further described in supplement 1, available at stacks.iop.org/ERL/10/035005/mmedia. In this study, we use daily discharge outputs from the hydrological HBV-96 (Lindström et al 1997) at Lobith, where the Rhine enters the Netherlands (see figure 1) and outputs from the Delft Continental Shelf Model (DCSM) (Zijl et al 2013) for the water level at Hoek van Holland. The HBV runoff model contains routines for snow accumulation and melt and soil moisture storage. Lateral water transport in the Rhine sub-basins, as well as transport in the Rhine River itself is simulated by the HBV routing module (Lindström et al 1997). Therefore, its response to rainfall accounts for antecedent water storage conditions within the basin (e.g. in soil moisture, groundwater and snow pack) as well as travel times through the river which introduces lags between meteorological events and hydrological events. The HBV sub-basins are displayed in figure 1, left panel.
The DCSM model effectively translates the effect of atmospheric pressure and prolonged wind fetch across a large domain on surge levels and combines these effects with water level anomalies due to tide. As surge levels are nonlinearly dependent on the water level without any consideration of wind, the inclusion of the tidal component is essential. The DCSM mesh is also displayed in figure 1, left and right panel.

Model performance
We assessed the model performance of both models by comparison of daily average discharges at Lobith and daily maximum water level time series at Hoek van Holland with observations. For the discharge values, we used the Nash-Sutcliffe efficiency (Nash and Sutcliffe 1970). The Nash-Sutcliffe efficiency (E) puts more emphasis on the performance of high flows, since it is a quadratic-error-based measure. For the water levels, we used the Spearman rank correlation (Spearman 1987). Rank correlation is more suitable for water levels because it prevents over-emphasizing on systematic bias (which may be due to the offset between a mesh grid cell value and the observation point), and rather focuses on the similarity in temporal pattern between the observation and simulation series (Kottegoda and Rosso 2008). The temporal patterns are of more interest since our study focusses on simultaneous occurrence of two phenomena.
This results in ρ = 0.70 for the sea water level and E = 0.92 for the discharge. This means that the modelling of the discharge is very accurate, while the modelling of the sea water levels is only reasonable. Possible reasons for this are a relatively coarse representation of wind fetch within the DCSM model due to the relatively coarse grid sizes used, the use of somewhat different drag relations in ERA-interim compared to DCSM, and the coarse spatial resolution of the ERAinterim forcing. Visual inspection revealed that the peaks found in the historical data are also found in the model results, so the extreme events are modelled reasonably. Another reason for the relative inaccuracy is that the data is aggregated to daily maxima. To further investigate model performance during extreme events, figure 2 shows empirical Gumbel and generalized extreme value distributions for the annual extreme values of both the model realizations and the measurements. The graph shows that the extreme values are well represented in the model results.
Due to the focus on time lag between values in the remainder of this paper, it is important to assure that the timing of discharge events is correctly simulated by the models, as discharge events in the Rhine are subject to travel times of several days. Therefore the arrival times of high river discharges at Lobith are compared against observations. From the arrival times of the 200 highest measured discharges, 158 are also in the highest 200 discharges found in the model, and all of them are in the highest 500. Thus it can be concluded that the arrival times of high discharges are quite well represented by the model. The atmospheric forcings for the different scenarios were generated using change factors as described in Bakker (2014). Thus the effect climate on wind forcing was included in the model. Sea level rise was added in a post-processing step. More information on the climate change scenarios can be found in supplement 2.

Assessing asymptotic behaviour and tail dependence
A measure that can be used to assess whether extreme values of bivariate stochastic variables (X, Y) show dependence towards the tail is (Ledford and Tawn 1997): where p is in the interval (0, 1) and F x and F Y are marginal distribution functions for paired realizations x and y of X and Y. If X and Y are asymptotically independent χ P is equal to If χ ≠ 0 P this implies that the variables X and Y are asymptotically dependent.
The value of χ P is a measure for the rate of dependence: if χ = 1 P it means that if variable X is extreme there is a 100% probability that the co-occuring value of Y is extreme as well., whereas if χ = 0.1 P this probability is 10%. Coles et al (1999) propose to use Chi ( χ) and Chibar ( χ ) plots for assessing asymptotic behaviour between the two extremes. These are given by: respectively. For asymptotically dependent variables it can be derived from equations (3) and (4) that χ > 0 P (with higher values meaning a stronger dependence) and χ = 1 P . To assess how well these measures can detect asymptotic dependence from observations we carried out an experiment in which 500 synthetic series of 5000 bivariate samples were generated from the threshold-excess logistic model of Tawn (1988). This bivariate model can be used to model asymptotically dependent variables (Zheng et al 2013), and is described by the following bivariate distribution function: , where X and Y are random variables, x and y are potential realizations of X and Y, G is the bivariate distribution function of X and Y, x 0 and y 0 are thresholds above which function G is valid and α is the model parameter that determines the rate of asymptotic dependence between X and Y. x and y are independent if α = 1, and fully dependent if α = 0. Table 1 gives derived estimates of χ P and χ P based on the 500 simulated series for different values of α. It can be seen that with the series length of 5000, it should be possible to detect dependence that can be described by an α up to at least 0.95. This shows we are able to detect even a relatively weak asymptotic dependence (α = 0.95) with the measures Chi ( χ) and Chibar( χ ) from a series of 5000 observations.   The dependence of two variables X and Y can also be assessed by comparing the upper quantiles of the variable Y with the same quantiles of variable Y, conditional on realizations of X. More concrete: the following values y 1 and y 2 are derived and, subsequently, compared: 2 With y 1 and y 2 being the thresholds for the quantile p of variable Y, and x is an increasing threshold for variable X. This method was for instance used by Geerse (2013). If X and Y are independent, y 1 = y 2 , whereas if they are asymptotically dependent it would result in y 2 ⩾ y 1 .

Dependence of sea water levels and Rhine discharges under present-day climate conditions
For all analyses the data was sampled for the winter half year, as including summer data may result in spurious dependence. A check to assess whether the number of remaining joint events is significant compared to the null hypothesis of independence was also carried out by Zheng et al (2013), results for cases with 6 day lag and without lag are shown in figure 3.
This shows that there is a significant number of joint events in the dataset for a lag of the discharge of six days, as the number of joint events in the dataset (red dot) is much higher than for the randomly sampled case (blue histogram). Therefore it seems that for the asymptotic dependence of water levels and discharges the time lag between peaks of the two random variables is very important. Figure 4 shows a  comparison between different percentiles of discharge (non)-conditional on the sea water level for the winter half year, for different time lags. This figure shows that there is especially a strong relationship between the two random variables for a discharge lag of six days. This corresponds to the physical reality of weather systems moving from the North Sea in South Westerly direction resulting in rainfall in the Netherlands and Germany, combined with the typical travel time of flood waves through the Rhine of about four days resulting in a lag of the proportion of six days. Figure 5 shows Chi and Chibar-plots for the lag of six days in the reference situation (blue line) and a comparison with the independent case (red line). In this plot some tail dependence is detected (blue line), as χ p ( ) and χ p ( ) are outside the upper confidence interval of the independent scenario (red dashed lines). From the graphs however, no conclusion can be drawn on whether this dependence is asymptotic. The maximum value found for χ P is approximately 0.1, which, if the BLTE-model were fit would result in an α of about 0.9.
We subsequently applied the method described by equations (6) and (7) using the same discharge lag of 6 days. This resulted in figure 6 for model results (left panel) and measurements (right panel). From this figure a difference in absolute values for the water level is clearly observed when comparing model realizations and measurements. This is most likely due to the limitations of the wind fetch schematization in the  . Realizations for both model realizations and measurements for combinations of sea water level and discharge during the winter half year. The dash-dotted lines represent unconditional 50th and 90th percentile values of the complete discharge series. The solid black and blue lines represent discharge values for the same percentiles, but conditional on water levels at Hoek van Holland larger than the value given on the x-axis. These lines include 95% confidence bounds based on bootstrap realizations of the respective datasets. storm surge model. However, for both measurements and model results a significant dependence between storm surges and Rhine discharges can be observed for extreme values, as in the figure, the black and blue lines and their corresponding triangles, as well as their confidence intervals, are clearly above the dash-dotted lines.

Influence of climate scenarios
The similarity of results between the measured time series and the model realizations gives confidence in the use of the model realization in the future climate conditions. When considering the influence of climate scenarios on the dependence of sea water levels and discharges it is especially interesting whether this dependence increases for more extreme climate change scenarios. Figure 7 shows the conditional and unconditional 90% quantile discharges for different climate scenarios.
From the figure it appears that the difference between conditional and unconditional discharge does not increase significantly (i.e. the mean of the conditional discharge in the reference scenario is within the range of the confidence intervals for the climate scenarios). This implies that climate change does not induce a more frequent simultaneous occurrence for both hazards, and therefore that the influence of climate change on coinciding high discharges and sea water levels is insignificant. The same figure but using absolute discharge values is shown in supplement 3.

Discussion
Historically joint probabilities of occurrence of sea water levels at Hoek van Holland and river discharges at the Rhine have not been taken into account in determination of design water levels (de Quay 1967). Recent research by Geerse (2013) and Kew et al (2013) suggests that there is a relationship between high sea water levels at the North sea and river discharges at the Rhine based on statistical analysis (Geerse 2013) and meteorological data analysis as a proxy for the North Sea water levels and Rhine river discharges (Kew et al 2013). In this paper a similar relationship is found but now based on physical models, simulating the discharge and water level processes themselves, and forced by meteorological data. Although the dependence between the two processes increases slightly in higher quantiles, no asymptotic behaviour towards more extreme quantiles can be observed. The relationship found only holds when a time lag between the water level and discharge series is introduced. A time lag of six days (i.e. moving the discharge series forward by six days) gives the highest dependence. The findings in this paper are also in line with the findings by Geerse (2013), who analysed historical data for the same locations and also observed an increase in water levels due to dependence. The DCSM v5 model used for generating high sea water levels at Hoek van Holland was able to generate peaks in accordance with the available historical data. However the height of the peaks was not always represented correctly, as the wind forcings used in this study may be biased with respect to the forcing used for model calibration. In detail, such differences in forcing may be due to difference in spatial resolution between the wind forcings used in this study and the forcing used for model calibration, and a somewhat different drag relation formulation between ERAinterim and DCSM. In general however the occurrences of storm surges are correctly represented in time, which makes the model suitable for this study. Compared to the study by Geerse (2013) the use of physical models also enables analysis of climate change effects on coinciding extremes.
Where Kew et al (2013) use proxy functions for extreme conditions, the model cascade used in this paper provides a better estimate of the relationships between the two parameters including possible shifts in time due to the memory within the modelled physical systems. Comparison of model results and measurements has shown that the timing of events is correctly represented in the models. Time lag between occurrence of extremes has been found to be very important for their relationship. In fact, the dependence between high sea water levels and high river discharges is strongest for a discharge lagged by 6 days with respect to the water level series at Hoek van Holland. In terms of policies and climate adaptation strategies a situation with this time lag is hardly relevant, but it shows nevertheless that the phenomenon of coinciding extreme discharges and sea water levels exists in this area. For smaller time lags of around 1 day less dependence is found, in accordance with Geerse (2013).
Regarding the quantiles chosen, the 90% quantiles used in this study, are hardly relevant for analysis of extreme events in the Netherlands, as safety standards are well above this level. However, given the limited period of the model runs (30 years), doing the analysis for higher quantiles was not feasible. Given the behaviour found here and in other studies (Geerse 2013, Kew et al 2013, it is however likely that the conclusions drawn on the basis of these quantiles can be extended to higher quantiles, in fact, dependence may even increase towards more extreme quantiles (Joe 1997, Diermanse andGeerse 2012). This could be further substantiated by using more samples, e.g. by the use of an ESSENCE ensemble (Sterl et al 2008) ran through the same model cascade.
In the comparison of the climate scenarios no significant change in dependence of sea water levels and river discharges is observed. The influence of the changed wind forcings on resulting water levels is very small; the main changes are due to sea level rise and changing precipitation in the Rhine basin leading to higher winter discharges in the future. There is however no visible influence of climate change on coincidence of extreme discharges and sea water levels. In the delta approach used to estimate future climate time series, changes in long time scale persistence were not accounted for, although climate models do resolve changes in such time scales to a certain extent (Johnson et al 2011). We therefore hypothesize that if such long-term persistence (e.g. inter-monthly and yearto-year variability in weather patterns) were resolved in the time series generation, more significant changes in the correlation structure would be found. Repeating the experiment with transformed time series that do accommodate such long-term climatological persistence may therefore be an important consideration for future study.

Conclusions
In general it can be concluded that high sea water levels at Hoek van Holland and high Rhine discharges at Lobith show significant dependence as was also shown in earlier studies (Geerse 2013, Kew et al 2013. However the dependence only appears with a time lag between the considered extremes: only if the sea water levels are related to a discharge at a later point in time, the dependence between discharge peaks and water levels at Hoek van Holland becomes stronger, with a maximum found at six days. For cases important for policy analysis, such as the case without any time lag between the two extremes, the relationship is rather weak. This makes the dependence found hardly relevant for policy making, as peaks do not tend to arrive in the area of interest at the same time. In smaller water systems, these time lags between high water levels and discharges may be much smaller such as found for the recent events in the North of the Netherlands as described in another contribution to this issue (van den Hurk et al 2014), as well as precipitation at the Australian coast (Zheng et al 2013(Zheng et al , 2014, where a clear relation between coincidence and the distance of rain gauges from the coast was found. These events emphasize the importance of further studying joint probabilities of extreme events, in particular where smaller scale systems are considered.
It should be noted that in this study, lower quantiles were used than those that match flood protection standards due to the limited sample size. A larger dataset is needed in order to study even higher tail quantiles.
Although future climate change leads to more extreme conditions in river discharges, our study could not identify significant change in the magnitude of the dependence between discharges and water levels during extreme conditions. We recommend repeating the analysis with a better representation of long time scale persistence within the used meteorological fields.