A joint probability approach using a 1-D hydrodynamic model for estimating high water level frequencies in the Lower Rhine Delta

The Lower Rhine Delta, a transitional area between the River Rhine and Meuse and the North Sea, is at risk of flooding induced by infrequent events of a storm surge or upstream flooding, or by more infrequent events of a combination of both. A joint probability analysis of the astronomical tide, the wind induced storm surge, the Rhine flow and the Meuse flow at the boundaries is established in order to produce the joint probability distribution of potential flood events. Three individual joint probability distributions are established corresponding to three potential flooding causes: storm surges and normal Rhine discharges, normal sea levels and high Rhine discharges, and storm surges and high Rhine discharges. For each category, its corresponding joint probability distribution is applied, in order to stochastically simulate a large number of scenarios. These scenarios can be used as inputs to a deterministic 1-D hydrodynamic model in order to estimate the high water level frequency curves at the transitional locations. The results present the exceedance probability of the present design water level for the economically important cities of Rotterdam and Dordrecht. The calculated exceedance probability is evaluated and compared to the governmental norm. Moreover, the impact of climate change on the high water level frequency curves is quantified for the year 2050 in order to assist in decisions regarding the adaptation of the operational water management system and the flood defense system.


Introduction
In the Lower Rhine Delta in the Netherlands, the Rhine and Meuse rivers run from the east and the south into the North Sea at Hook of Holland, to the Haringvliet in the west and into the Lake IJsselmeer in the north.The area is a center of high economic activity and maritime transportation with a dense population, which is at risk of being flooded by river or by sea, or by both.The water level in this transitional area is influenced by the upstream river flows, the downstream sea level as well as the operation of the existing controllable structures.At the upstream boundary, the Rhine flow comes from rainfall-runoff and snowmelt in the Alps; the Meuse flow is mainly determined by rainfall in France and Belgium.At the downstream boundaries, the extreme still water level (excluding waves) arises from a combination of the astronomical tides and the meteorologically induced storm surge components.In this article the extreme still water level is the so-called "storm surge".Astronomical tides are driven by astronomical forces and are deterministic, while the wind induced storm surges occur stochastically, driven by meteorological forcing.To protect the delta from sea flooding, the estuary can be closed off from the sea by large dams and controllable gates and pumps.Also, along the rivers controllable structures have been constructed in order to regulate the flow and water levels.The operation of these structures was involved in the operational water management system of the Netherlands (van Overloop, 2009).
It is important to evaluate the high water level frequencies; first to evaluate the flood risk, and second, to support the future flood defense design.The resistance of the river dikes against flooding is captured in the design water level; for instance, for the economically important city of Rotterdam the design water level is regarded as the water level with the exceedance frequency of 1/10 000.However, climate change and human interventions complicate the computation of these high water level frequencies.For example, being a very significant human intervention, the change in the operational water management system has resulted in nonhomogeneous high water level observations.
The non-homogeneous extreme observations derived from climate change and human interventions cannot be used to estimate the high water level frequencies.Instead, a joint probability approach, using a 1-D hydrodynamic model, can be used to estimate the high water level frequencies (Mantz and Wakeling, 1979;Acreman, 1994;Gorji-Bandpy, 2001;Samuels and Burt, 2002;Adib et al., 2010).As a given high water level at a transition location may result from a number of combinations of sea level and upstream fluvial flow and from how the operational water management system reacts to the situation at hand, the occurrence of all these combinations together determines the frequency of the given water level.
The first application of a joint probability approach in the Lower Rhine Delta dates back to 1969.Van der Made (1969) divided three joint probability distributions for three individual categories: high sea levels and normal discharges, normal sea levels and high discharges, high discharges and high sea levels.For each category the joint probability distribution was estimated from the observations of the peak values of the sea level and the Rhine flow in the same day.
However, the above joint probability approach only considered the peak values of the sea level and the Rhine flow, and assumed the other associated variables, for example the surge duration, to be pre-determined.These associated variables also play an important role in the operational water management system, and therefore influence the high water level frequencies in the delta.For example, after the closure of the Maeslant barrier and the Haringvliet dam, the water level within the delta depends on the fluvial flow, the storage capacity and the closure duration (Zhong et al., 2012).Recent research (De Michele et al., 2007;Wahl et al., 2012) has shown how to include more variables in the probabilistic analysis of the hydrodynamic boundary conditions.In this paper, more associated variables will be taken into account in the high water level frequency estimation.1-D hydrodynamic numerical models can be applied to assess the changing water level frequencies caused by different changes, for example changes in mean sea level rise, peak river flows and the construction and operation of controllable hydraulic structures.The detailed models can produce accurate water levels at transitional locations by examining the interaction of sea level, fluvial flow and infrastructure operations, but their computational time restrains the application of those models in a Monte Carlo simulation method.For example, only a very limited number of sampling scenarios were used by Mantz and Wakeling (1979) and Samuels and Burt (2002).However, for accurate Monte Carlo simulations a very large number of stochastic sampling scenarios are necessary.For that reason, in this research, a strongly simplified 1-D hydrodynamic model of the Lower Rhine Delta is applied.
Future climate change will affect the high water level frequencies in the Lower Rhine Delta.For the Rhine flow, climate change is expected to increase winter precipitation with earlier snowmelt (Middelkoop and Kwadijk, 2001), which will lead to an increase in the frequency and magnitude of extreme Rhine flows (Hooijer et al., 2004;Pinter et al., 2006;Linde et al., 2010).For mean sea levels along the Dutch coast, a range of 0.15 to 0.35 m rise until 2050 and a range of 0.35 to 0.85 m rise until 2100, corresponding to the reference year of 1990, are commonly used extrapolation values (van den Hurk et al., 2006;Second Delta Commission, 2008).In fact, the relative mean sea level rise will be larger when taking mean land subsidence due to glacial isostasy and subsoil compaction into consideration.The effects of climate change on the characteristics of the wind induced surge along the Dutch coastline was investigated and no evidence of significant changes was detected (Sterl et al., 2009), and, hence, we can assume that the characteristics are not influenced by climate change.Although there are still inherent uncertainties in the prediction of climate change on the hydraulic boundary conditions within climate change scenario studies, it can be assumed that the future changes in high water level frequencies can be assessed by applying an appropriate climate change scenario.In this paper estimates of mean sea level rise and increases of peak Rhine discharge in the future scenario of 2050 are included to assess the future high water level frequencies.The results can assist in adapting the operational water management system to control the negative effects of climate change in the Lower Rhine Delta.
This paper aims to assess the high water level frequencies in the Lower Rhine Delta and to identify the exceedance frequencies of the design water level.The paper starts with the introduction of the method and an illustration of the model, followed by the joint probability analysis.The results and discussion are presented, followed by the conclusions and future research recommendations.

Method
Applying the joint probability approach using a 1-D hydrodynamic model to estimate the high water level frequencies in the estuary delta is illustrated in Fig. 1.The importance sampling Monte Carlo simulation method will help generate a large number of stochastic scenarios as inputs for the above model.The outline of the method is as follows: -from historical observations, three categories are selected in order to separate different combinations of sea water levels and Rhine flows that may cause high water levels in the Lower Rhine Delta; -for each category a corresponding joint probability distribution is estimated;

Nat
-a large set of stochastic scenarios are generated by means of the importance sampling Monte Carlo simulation for each category (10 000 per category); -1-D hydrodynamic model including infrastructure operations is simulated with the scenarios of the three categories and result in the same number of high water levels (10 000) at transitional locations in the delta; and -the high water level frequency curves are computed.
In addition, the exceedance probabilities of the design water levels at Rotterdam and Dordrecht are calculated from the resulting peak water levels of the simulations.

The 1-D simplified hydrodynamic model of the Lower Rhine Delta
The Rhine Delta is a system of inter-connected rivers, canals, reservoirs, and adjustable structures as can be seen in Fig. 2 (left).A strongly simplified 1-D hydrodynamic model was used to simulate the complex flows (van Overloop, 2009) The hydrodynamic characteristics of the delta are mainly governed by discharge of the rivers Rhine (Lobith: node 14 in Fig. 2), Meuse (Borgharen: node 1) and by the water level at the sea boundaries (Hook of Holland: node 36 and Haringvliet: node 29).In the model calculations, the sea level at Haringvliet is assumed to be the same as at Hook of Holland.As the focus of the research is on the Lower Rhine Delta surrounding the economically important cities of Rotterdam and Dordrecht, the other sea level boundaries in the North (node 30, 33, 35), which do not affect Rotterdam and Dordrecht, are set to 0 m m.s.l.during the running of the model.At Rotterdam (node 24) and Dordrecht (node 22) the high water level frequencies are estimated.In addition, the dike heights along the rivers are assumed to be high enough so that no overflowing or breaching can occur.

The joint probability analysis
This section starts with the description of the observation data, followed by the estimation of the joint probability distributions for three categories.The data is shown in Table 1.
The behavior of the River Rhine, Meuse and the sea conditions may change considerably over the longer periods due to artificial or natural causes.Commonly, three popular statistical tests are employed to check whether a trend exists in an observation time series.The Mann-Kendall test can be applied to assess the significance of trends in hydrometeorological time series such as stream flow, temperature and precipitation (Mann, 1945;Kendall, 1975).The Spearman's rho test can also be used to detect monotonic trends in a time series (Lehmann, 1975;Sneyers, 1990).The Wilcoxon rank sum test can test if abrupt shifts exist in a time series (Wall, 1986).No significant trends or shifts have been detected with these three tests in the annual maximum series of sea level and Rhine and Meuse discharges, except that a least squares linear regression suggests a gradual increase of 0.20 m mean sea level rise per century.This result is in  The Rhine Delta is a system of inter-connected rivers, canals, reservoirs, and adjustable 7 structures as can be seen in Fig. 2 (Left).A strongly simplified 1-D hydrodynamic model 8 was used to simulate the complex flows (van Overloop, 2009).The model used a large 9 grid size of 20 km and constant bed slopes.The large water bodies, such as the IJsselmeer, 10 were modeled as reservoirs with a level-area description.line with previous research (van Gelder, 1996;Zhong et al., 2012).
The division into three categories is based on thresholds of the peak surge residual and the peak of Rhine flow occurring in the same day: 1.00 m in Hook of Holland and 6000 m 3 s −1 at Lobith.The surge residual is the difference of the observed sea level and the predicted astronomical tide level.This threshold value for the peak surge residual is chosen for two reasons: first of all, this value is related to the operation of the Maeslant Barrier.The peak surge residual of 1.0 m, coinciding with the high astronomical tide level and a high Rhine flow, could make the Rotterdam water level exceed the critical value of 3.0 m m.s.l.(the decision level of the closure of the barrier).Secondly, the threshold of 1.0 m was applied in the estimation of the frequency of the wind induced surge peak level (Bijl, 1997).The threshold of 6000 m 3 s −1 for Rhine discharge is determined by means of three reasons: first of all, this value is related to the operation of the Maeslant Barrier (Bol, 2005).Secondly, this value is related to the floodplains inundated along the lower Rhine branch.A discharge slightly exceeding 6000 m 3 s −1 was more or less assumed to be the critical value which resulted in the highest floodplains inundated (Kwadijk and Middelkoop, 1994).Thirdly, the threshold value of 6000 m 3 s −1 was applied by Chbab (1995), with the generalized Pareto distribution to estimate the frequencies of high Rhine flows.In this article the application of this threshold, as well as the fitted generalized Pareto distribution function, leads to a Rhine design discharge (with a probability of 1/1250 per year) of 15 250 m 3 s −1 , which is comparable to commonly used values.
The selected events from 1939 to 2009 in Fig. 3 are used to estimate the joint probability distributions of three categories.The biggest sea flooding in 1953 is missing from the website database.Instead, the information of this surge residual comes from Gerritsen (2005), in which the peak surge residual is 3.00 m and the duration is 50 h.7 design discharge (with a probability of 1/1,250 per year) of 15250 m 3 /s, which is comparable to commonly used values.The selected events from 1939 to 2009 in Fig. 3 are used to estimate the joint probability distributions of three categories.The biggest sea flooding in 1953 is missing from the website database.Instead, the information of this surge residual comes from Gerritsen (2005), in which the peak surge residual is 3.00 m and the duration is 50 hours.(Vrijling and Bruinsma, 1980;Praagman and Roos, 1987).This method was introduced and further validated for the station of Hook of Holland.

Storm surge and normal Rhine flow
The historical events of storm surges coinciding with normal upstream flows are shown in Category I of Fig. 3.The probability distribution of the storm surges in the Eastern Scheldt was estimated by separating the astronomical tide and the wind induced storm surge (Vrijling and Bruinsma, 1980;Praagman and Roos, 1987).This method was introduced and further validated for the station of Hook of Holland.
From a statistical point of view, the occurrence of the astronomical tide component is independent of the occurrence of the wind induced storm surge component at the mouth of the Lower Rhine Delta.However, these two components can interact with each other when they propagate into the delta (de Ronde, 1985).Their nonlinear interaction generally increases the surge height at a rising astronomical tide and decreases the surge height at a high astronomical tide (Bijlsma, 1989).Quantifying the nonlinear effect is beyond the scope of this study.For the sake of convenience, it can be assumed that the wind induced storm surge is independent of the astronomical tide as indicated in Fig. 4.
These surge residual curves are taken into the probability analysis with two parameters: peak surge residual h smax and surge duration T s .The probability distributions of these two parameters are applied to simulate many pseudo surge residual curves with an appropriate shape function.The astronomical tide curves can also be simulated by the same logic.As a result, the simulated surge residual curves and the simulated tide curves can be linearly combined into the simulated sea level curves.
In order to estimate the surge curve in Hook of Holland, 300 extreme surge residuals in Category I are analyzed in Fig. 3.The observed peak surge residuals and associated durations are plotted in Fig. 5. Their linear correlation coefficient is 0.0474, and therefore they are assumed to be linearly independent.For a surge residual curve at Hook of Holland, the peak surge residual and duration are generated and constrained by complex physical factors like the offshore surge, the shallow water depth, the interaction between tide and in Category I of Fig. 3.The probability distribution of the storm surge 23 Scheldt was estimated by separating the astronomical tide and the wind 24 surge (Vrijling and Bruinsma, 1980;Praagman and Roos, 1987).Th 25 introduced and further validated for the station of Hook of Holland.26 27  These surge residual curves are taken into the probability analysis with two param 12 surge residual h smax and surge duration T s .The probability distributions of 13 parameters are applied to simulate many pseudo surge residual curves with an 14 shape function.The astronomical tide curves can also be simulated by the 15 As a result, the simulated surge residual curves and the simulated tide curves can 16 combined into the simulated sea level curves.17 18 300 extreme surge residuals in Category I are analyzed in Fig. 3 in order to estima 19 curve in Hook of Holland.The observed peak surge residuals and associated d 20 plotted in Fig 5 .Their linear correlation coefficient is 0.0474, and therefore they 21 to be linearly independent.For a surge residual curve at Hook of Holland, the 22 residual and duration are generated and constrained by complex physical fact 23 offshore surge, the shallow water depth, the interaction between tide and surge, et 24 from a statistical point of view, the assumption of independence between the 25 residual and the duration is acceptable.26 The surge residual curve can be approximated by a squared cosine function.I 30 comparison between the observed surge residual curves and the simulated cu 31 symmetric cosine function of six extreme surge events agrees with this 32 assumption.In Fig. 7 the surge residual curve from the big sea flooding in 1953 33 a symmetric curve (Gerritsen, 2005).34 35 The design surge residual curve can be derived from the observed surge residu 36 surge, etc.However, from a statistical point of view, the assumption of independence between the peak surge residual and the duration is acceptable.
The surge residual curve can be approximated by a squared cosine function.In Fig. 6 the comparison between the observed surge residual curves and the simulated curves by the symmetric cosine function of six extreme surge events agrees with this reasonable assumption.In Fig. 7 the surge residual curve from the big sea flooding in 1953 also showed a symmetric curve (Gerritsen, 2005).
The design surge residual curve can be derived from the observed surge residual curves.
where h smax stands for the peak value of the surge residual, and its unit is m; T s is the duration of the surge residual, and its unit is hours.Here we assume that the surge peak occurs when t is 0. The generalized Pareto distribution and the Weibull Distribution fit the distributions of peaks and durations of these selected surge residuals, respectively.In this research, all parameters of distributions are estimated by the maximum likelihood method.
The semi-diurnal astronomical tide in Hook of Holland, is almost symmetric, and can therefore be approximated by a s T where h smax stands for the peak value of the surge residual, and its unit is m; T s is the duration of the surge residual, and its unit is hours.Here we assume that the surge peak occurs when t is 0.  (1) re h smax stands for the peak value of the surge residual, and its unit is m; T s is the ation of the surge residual, and its unit is hours.Here we assume that the surge peak urs when t is 0. General Pareto Distribution and the Weibull Distribution fit the distributions of ks and durations of these selected surge residuals, respectively.In this research, all ameters of distributions are estimated by the Maximum Likelihood Method.semi-diurnal astronomical tide in Hook of Holland, is almost symmetric, and can efore be approximated by a sinusoidal-curve and modeled as a periodical fluctuation sinusoidal curve and modeled as a periodical fluctuation of the water level h a with a period of 12.4 h and with amplitude of h HW − h LW .Where h HW is the high tide level; h LW is the low tide level; their unit is m m.s.l.; u is the time shift between peaks of tide and surge.Figure 8 shows that the simulated tide level from the sinusoidal function represents the tide well.
As a consequence of the assumed independency of the tide and the wind induced storm surge, the time shift between peaks u fits a uniform probability density function.
10 of the water level ha with a period of 12.4 hours and with amplitude of hHW -hLW.Where 1 hHW is the high tide level; hLW is the low tide level; their unit is m MSL; u is the time shift 2 between peaks of tide and surge.Fig. 8 shows that the simulated tide level from the 3 sinusoidal function represents the tide well.4 5 2 ( ) sin( ( )) 2 12.4 2 As a consequence of the assumed independency of the tide and the wind induced storm 10 surge, the time shift between peaks u fits a uniform probability density function.Time shifts u larger than 12.4 h are irrelevant, so considering a symmetrical shape, the probability density function of u becomes In conclusion, the storm surge water level is where h 0 is mean sea level.ompanying low Rhine and Meuse flows can be assumed to be constant during the rges' period, which is supposed not to influence the water levels in Rotterdam in alculation.ompanying low Rhine and Meuse flows can be assumed to be constant during the urges' period, which is supposed not to influence the water levels in Rotterdam in alculation. is estimated by the one year data of high astronomical tide levels that is derived from the harmonic analysis of water level observations (see Fig. 9).In Fig. 10 the low tide level (h LW ) is approximately a linear function of h HW .
The probability distribution of the associated normal upstream flow can be estimated by the accompanying dailyaverage Rhine and Meuse flows.The Gaussian copula function represents the dependence structure between Rhine discharge Q r and Meuse discharge Q m , where the marginal distributions fit log-normal distribution for Q r and the gamma distribution for Q m .Figure 11 indicates that the Gaussian copula presents well the dependence between daily-average Rhine and Meuse flows.The Gaussian copula dependence structure as well as the marginal distributions is well applied to simulate the upstream flows for Category I where the few occurrences of Rhine flows exceeding 6000 m 3 s −1 are maximized at 6000 m 3 s −1 .
The accompanying low Rhine and Meuse flows can be assumed to be constant during the storm surge period, which is not supposed to influence the water levels in Rotterdam in the model calculation.surge residual in Fig. 3) is applied to detect the wind induced storm surge event 6 average the annual-average number of events of 4.38 is chosen.The number o 7 occurring in one year fits a Poisson distribution and the parameter  is 4.38.The 8 GPD process can be transformed to the GEV distribution for annual maxima (show 9 appendix).The detailed information is available from Smith (2004).The probabil 10 refers to the exceedance probability of a specific * R h in one year.11 where I is an indicator function and h R is calculated from the specific input par 12 using the hydrodynamic model.In a event in Category I, the upstream Rhine 13 independent of the storm surge (Dantzig et al., 1960;Jorigny et al., 2002).14

High Rhine flow and normal sea water level 15
The events of high Rhine flows coinciding with normal sea water levels are s 16 Category II of Fig. 3.This category focuses on this kind of combinations which 17 extreme water level in Rotterdam.It is assumed that the wind induced surge com 18 can be discarded when the peak level of the surge residual is lower than 1.0 m.T 19 in this kind of combinations the astronomical tide is the only component to be co 20 in the sea water level.21 22 The high Rhine flows come from large scale storm depressions which probably al 23 about the associated high Meuse flows.The Gumbel copula is applied to desc 24 dependency, as it exhibits a stronger dependency in the positive tail.The as 25 Meuse flows are selected at the same day when the Rhine peaks occur.A Genera 26 Distribution and a Lognormal Distribution fit the selected Rhine and Meus 27 respectively.28 Equation ( 5) computes the exceedance probability of a certain Rotterdam water level h * R in one year for Category I.The peak over threshold (POT) method (1.0 m for the peak surge residual in Fig. 3) is applied to detect the wind induced storm surge events and on average the annual-average number of events of 4.38 is chosen.The number of events occurring in one year fits a Poisson distribution and the parameter λ is 4.38.The Poisson-GPD process can be transformed to the GEV distribution for annual maxima (shown in the appendix).The detailed information is available from Smith (2004).The probability p(h * R ) refers to the exceedance probability of a specific h * R in one year.
where I is an indicator function and h R is calculated from the specific input parameters using the hydrodynamic model.In a event in Category I, the upstream Rhine flow is independent of the storm surge (Dantzig et al., 1960;Jorigny et al., 2002).

High Rhine flow and normal sea water level
The events of high Rhine flows coinciding with normal sea water levels are shown in Category II of Fig. 3.This category focuses on this kind of combinations which result in extreme water levels in Rotterdam.It is assumed that the wind induced surge component can be discarded when the peak level of the surge residual is lower than 1.0 m.Therefore, in this kind of combination the astronomical tide is the only component to be considered in the sea water level.The high Rhine flows come from large scale storm depressions which probably also bring about the associated high Meuse flows.The Gumbel copula is applied to describe this dependency, as it exhibits a stronger dependency in the positive tail.The associated Meuse flows are selected from the same day when the Rhine peaks occur.A generalized Pareto distribution and a log-normal distribution fit the selected Rhine and Meuse flows, respectively.
where the relationship between the Gumbel copula parameter α and Kendall's tau τ is also shown.α is estimated as 1.7158; F r , is the marginal distribution of high Rhine flow; F m , is the marginal distribution of the associated Meuse flow; The Chi-square (χ 2 ) test is used to determine the goodness of fit between observed data with expected values derived from the Gumbel copula.The calculated value of χ 2 being 27.8 is far below the critical value of 61 for 47 degrees of freedom at the significance level of 0.5 %.In addition, the good fit is shown in Fig. 12.
High Rhine and Meuse flow curves can be generated by the design hydrographs (Parmet et al., 2002a, b) multiplied by the ratio between the generated values and the design peak values.
Equation ( 7) shows the exceedance probability of a certain Rotterdam water level h * R in one year for Category II.

Storm surge and high Rhine flow
The very limited number of observations of the joint high surge residual and high Rhine flow events in Category III is not appropriate for estimating the joint probability distribution.A rather simple method is introduced.In the historical observations, 9 events occurred in Category III and therefore it is assumed that the occurrence probability of a combination event in Category III is 9/70 per year.The marginal distributions of the peak surge residual, the surge duration, the astronomical tide and the high Rhine flow are the same as Category I and II, respectively.In a combination event, the high peak surge residual is assumed to be independent to high Rhine flow (Dantzig et al., 1960;Jorigny et al., 2002).
The Chi-square ( 2 ) test is used to determine the goodness of fit between observ 6 with expected values derived from the Gumbel copula.The calculated value of  7 27.8 is far below the critical value of 61 for 47 degrees of freedom at the sign 8 level of 0.5%.In addition, the good fit is shown in Fig. 12. 9 10

Monte Carlo simulation
A large number of boundary stochastic scenarios need to be generated based on the joint probability distribution for each category.Then the 1-D model can run using these scenarios and the outputs are the same number of peak water levels at locations of interest in the Lower Rhine Delta.The resulting series of peak water levels as well as the accompanying occurrence probabilities can be transformed into the high water level frequency curves in the delta.Due to the important economic role only the results at Rotterdam and Dordrecht are presented.Importance sampling is applied to reduce the number of samples in the Monte Carlo simulation but still get sufficiently accurate estimations (Glynn and Iglehart, 1989;Roscoe and Diermanse, 2011).In the Monte Carlo simulations, the exceedance probability P f of a specific h * R is simply taken to be n f /n, where n f is the number of samples which lead to h R ≥ h * R , and n is the total number of generated samples.In importance sampling, the number of samples which lead to h R ≥ h * R increases largely because boundary inputs are not generated from their original probability distributions, but from alternative distributions which focus on exceedance of the critical water level at Rotterdam h * R .In this study, we have used normal distributions for the most important input variables h smax , T s , Q r (high Rhine flow) and Q m (high Meuse flow), centered around the values that lead to h * R .Note that for different h * R , the corresponding normal distributions are different in order to locate around the area leading to h R ≥ h * R .The changes in distributions need to be compensated for.
where P f (h * R ) is the exceedance probability of a specific Rotterdam water level h * R ; n is the total number of samples; I is an indication function inside which the input parameters are generated from the distributions g, f stands for the original probability density distributions of related variables and g is the corresponding normal density distributions.
In the importance sampling method only input parameters referring to h smax , T s , Q r (high Rhine flow) and Q m (high Meuse flow) applied the new normal distributions instead of the original distributions.The other input parameters were sampled from their original probability distributions.Generally there are no upper bounds for these normal distributions.
To estimate the high water level frequencies and exceedance probabilities of the design water levels, 10 000 boundary events were generated with the importance sampling method and the model output were 10 000 maximum water levels at Rotterdam and Dordrecht for each of the three categories.In order to test whether 10 000 simulations was enough to get consistent results, another two groups of 10 000 simulations were generated to compare the difference.These differences were found negligible.
To estimate the effect of climate change on the high water level frequencies, an increase of the peak discharge of the Rhine as well as an increase of the mean sea level in the scenario of 2050 is studied.The mean sea level rise is assumed to be 0.35 m (van den Hurk et al., 2006) and the peak Rhine discharge increases by 10 % in reference to the year 2000 (Jacobs et al., 2000).In a second set of the Monte Carlo simulations, the input boundary conditions valid for scenario 2050 are generated by simply re-scaling the present boundary variables.

Exceedance probability of the present design water level
The design water level is crucial for the design, construction and maintenance of the flood defense system.According to the Dutch law, the design water level in Rotterdam is regarded as the water level with the exceedance frequency of 1/10 000; the design water level in Dordrecht is with the exceedance frequency of 1/2000.The present design water level is 3.6 m m.s.l. for Rotterdam and 3.0 m m.s.l. for Dordrecht (Ministerie van Verkeer and Waterstaat, 2007).
The exceedance probabilities of the design water levels estimated in our method are a little lower than the official values based on Hydra-B (Ministerie van Verkeer and Waterstaat, 2007).The results are listed in Table 2.
From Table 2 the exceedance probability of the design water level is 0 for Rotterdam in Category I.The result indicates that in current conditions the delta area can be protected from storm surges until the year of 2050.This agrees with the design standard of the Maeslant Barrier, one key hydraulic structure at the mouth of the delta.Note that the result depends on the assumption that the operation of the storm surge barriers and dams at the mouth of the delta does not fail during the storm surge events.It is believed that taking the failure of the operation into consideration could result in a higher exceedance probability in Rotterdam and Dordrecht for Category I. Detailed research on the failure probabilities of the hydraulic structures at the mouth is urgently required.
The exceedance probability in Category II is also 0 in Rotterdam for both the present and the year 2050.The exceedance probability in Dordrecht is also very low.The results indicate that high fluvial flow (Rhine and Meuse flow) has very limited influence on the extreme water levels of the downstream of the delta.Instead, the extreme upstream fluvial flow could easily result in the breaching or overflowing in the Dutch Upper Rhine Delta, which agrees with the nearcatastrophic floods of 1993 and 1995 (Engel, 1997).
The exceedance probability in Category III is still far lower than the official standard 10 −4 in Rotterdam for the present and the year 2050, while the exceedance probability in Dordrecht is higher than the official standard 1/2000 for the year of 2050.Moreover, the sum of the exceedance probability in three categories shows that the Dutch Lower Rhine Delta complies with the required norm for flood safety, except for the Dordrecht in future climate scenario of 2050.
The results show that the exceedance probability of the design water level is much higher in Category III than in Category I and Category II.It indicates that the combinated events of the storm surge and the high Rhine flow become the main source of floods in the Lower Rhine Delta.The high water level frequency curves derived from these combinated events can be drawn in Rotterdam and Dordrecht.

High water level frequency curves
The high water level frequency curves derived from the combination events of storm surges coinciding with high Rhine flows are shown in Fig. 13.
The future high water level frequency curves (the dash lines in Fig. 13) are about 0.2 to 0.4 m higher than the present curves (the solid lines in Fig. 13) in Rotterdam and Dordrecht.It indicates that climate change will lead to more extreme events which increase the high water level frequency in the future.The differences between the present and future high water level frequency curves are quantified in order to provide an indication for the further adaptation of the   The future high water level frequency curves (the dash lines in Fig. 13) are about 0.2 to 0.4 m higher than the present curves (the solid lines in Fig. 13) in Rotterdam and Dordrecht.It indicates that climate change will lead to more extreme events which increase the high water level frequency in the future.The differences between the present and future high water level frequency curves are quantified in order to provide an indication for the further adaptation of the operational water management system and the flood defense system.The present high water level frequency curve in Rotterdam shows that the exceedance probability of 3.0 m MSL is below 10 -4 .It is attributed by the controlled structures as indicated in Fig. 2 and the operational water management behind then.During the combination events of the storm surge and high Rhine flow, the Maeslant Storm Surge Barrier and the Haringvliet sluices are closed in order to prevent the sea water from propagating into the delta, and after the closure the water level in Rotterdam and Dordrecht will increase due to the Rhine and Meuse flows coming into the delta.The mouth of the delta is open again to discharge fluvial water into the sea after the storm.The closing decision level is 3.0 m MSL in Rotterdam and 2.7 m MSL in Dordrecht.To avoid extreme high water levels from the events in Category III, the construction of new structures needs exploration and the present operational water management system requires adaptation in the future.

Conclusion and future research
This paper presents the application of the joint probability sampling approach coupled with a simplified 1-D hydrodynamic model to assess the exceedance probability of the present design water level in Rotterdam and in Dordrecht.The high water level frequency operational water management system and the flood defense system.
The present high water level frequency curve in Rotterdam shows that the exceedance probability of 3.0 m m.s.l. is below 10 −4 .It is attributed by the controlled structures as indicated in Fig. 2 and the operational water management behind then.During the combination events of the storm surge and high Rhine flow, the Maeslant Storm Surge Barrier and the Haringvliet sluices are closed in order to prevent the sea water from propagating into the delta, and after the closure the water level in Rotterdam and Dordrecht will increase due to the Rhine and Meuse flows coming into the delta.The mouth of the delta is open again to discharge fluvial water into the sea after the storm.The closing decision level is 3.0 m m.s.l. in Rotterdam and 2.7 m m.s.l. in Dordrecht.
To avoid extreme high water levels from the events in Category III, the construction of new structures needs exploration and the present operational water management system requires adaptation in the future.

Conclusion and future research
This paper presents the application of the joint probability sampling approach coupled with a simplified 1-D hydrodynamic model to assess the exceedance probability of the present design water level in Rotterdam and in Dordrecht.The high water level frequency complies with the required norm for safety in the present.However, the threats of high water levels exceeding the design water level are still exposed to both cities at a low probability mainly due to the combination events of storm surges and high Rhine flows.
The new method enables assessment of high water level frequencies in a changing environment with associated effects from climate change and operation of the infrastructures.In the future climate change will lead to more extreme events and increase the high water level frequency in the Lower Rhine Delta.Moreover, the future development of local economy and urbanization will increase the flood induced damage when floods occur (te Linde et al., 2011).Therefore the adaption measures are urged.The adaptation of the present operational water management system was proposed by van Overloop et al. (2010).The method in the article, based on the Model Predictive Control method (in brief MPC), can be applied to investigate the effect of the adaption measure on reducing the high water level frequency in the delta.
In future research, the failure probability of the operation of these controllable hydraulic structures should be further incorporated.In addition, the statistical uncertainty in the joint probability approach needs to be investigated.
Fig. 1.Flowchart of the methodology Fig. 1.Flowchart of the methodology.
s.l.) 1939-2009 time unit is the same as the above sea level Lobith Rhine discharge (m 3 s −1 ) 1901-2009 daily-average discharge Borgharen Meuse discharge (m 3 s −1 ) 1911-2009 daily-average discharge Note: the source of these data is the Rijkswaterstaat website: http://www.rijkswaterstaat.nl/waterbase.Left) Description of the Rhine delta with existing operated structures and (Right) 4 overview of the simplified 1-D hydrodynamic model (van Overloop, 2009).5 6 Fig. 2. (left) Description of the Rhine delta with existing operated structures and (right) overview of the simplified 1-D hydrodynamic model (van Overloop, 2009).
Fig. 3. Selected events: Category I, storm surge and normal Rhine flow; Category II.high Rhine flow and normal sea water level; Category III.storm surge and high Rhine flow 4.1 Storm surge and normal Rhine flow The historical events of storm surges coinciding with normal upstream flows are shown in Category I of Fig. 3.The probability distribution of the storm surges in the Eastern Scheldt was estimated by separating the astronomical tide and the wind induced storm surge(Vrijling and Bruinsma, 1980;Praagman and Roos, 1987).This method was introduced and further validated for the station of Hook of Holland.

Fig. 3 .
Fig. 3. Selected events: Category I, storm surge and normal Rhine flow; Category II, high Rhine flow and normal sea water level; Category III, storm surge and high Rhine flow.

Fig. 4 .
Fig. 4. Variation with time of the extreme still water level. 11 residual (m) Surge duration (hours) 27 Fig. 5.The peak surge residuals and associated durations 28 29
Fig. 6. (Left) The historical surge residual curves and the simulated design curves; (Right) the correlation coefficient squared R 2

Fig. 6 .
Fig. 6. (left) The historical surge residual curves and the simulated design curves; (right) the correlation coefficient squared R 2 .

. 6 .
(Left) The historical surge residual curves and the simulated design curves; (Right) the correlation coefficient squared R 2 (Left) The surge residual curve of the largest flood in 1953 and the design simulated curve; (Right) the correlation coefficient squared R 2

Fig. 7 .
Fig. 7. (left) The surge residual curve of the largest flood in 1953 and the design simulated curve; (right) the correlation coefficient squared R 2 .

Fig. 9 .0
Fig. 9. High tide level fitted by the Normal Distribution

Fig. 9 .0
Fig. 9. High tide level fitted by the normal distribution.
Results from a graphical based goodness fit of the Gaussian Copula simul 2 3 Equation (5) computes the exceedance probability of a certain Rotterdam water 4 in one year for Category I.The Peak over Threshold (POT) method (1.0 m for 5

Fig. 11 .
Fig. 11.Results from a graphical based goodness fit of the Gaussian copula simulation.

Fig. 12 .
Results from a graphical based goodness fit of the Gumbel Copula simula 11 12 High Rhine and Meuse flow curves can be generated by the design hydrographs 13 et al., 2002a; Parmet et al., 2002b) multiplied by the ratio between the generated 14 and the design peak values.15 16 Equation (7) shows the exceedance probability of a certain Rotterdam water lev 17 one year for Category II.

Fig. 12 .
Fig. 12. Results from a graphical based goodness fit of the Gumbel copula simulation.
Fig. 13.The high water level frequency curves due to the combination events of storm surges and high Rhine flows (category III) in a. Rotterdam.b.Dordrecht The future high water level frequency curves (the dash lines in Fig.13) are about 0.2 to 0.4 m higher than the present curves (the solid lines in Fig.13) in Rotterdam and Dordrecht.It indicates that climate change will lead to more extreme events which increase the high water level frequency in the future.The differences between the present and future high water level frequency curves are quantified in order to provide an indication for the further adaptation of the operational water management system and the flood defense system.The present high water level frequency curve in Rotterdam shows that the exceedance probability of 3.0 m MSL is below 10 -4 .It is attributed by the controlled structures as indicated in Fig.2and the operational water management behind then.During the combination events of the storm surge and high Rhine flow, the Maeslant Storm Surge Barrier and the Haringvliet sluices are closed in order to prevent the sea water from propagating into the delta, and after the closure the water level in Rotterdam and Dordrecht will increase due to the Rhine and Meuse flows coming into the delta.The mouth of the delta is open again to discharge fluvial water into the sea after the storm.The closing decision level is 3.0 m MSL in Rotterdam and 2.7 m MSL in Dordrecht.To avoid extreme high water levels from the events in Category III, the construction of new structures needs exploration and the present operational water management system requires adaptation in the future.

Fig. 13 .
Fig. 13.The high water level frequency curves due to the combination events of storm surges and high Rhine flows (category III) in (a) Rotterdam.(b) Dordrecht.
The characteristics of the high tide level (h HW ) at Hook of Holland can be captured in a normal distribution which 11 erived from the harmonic analysis of water level observations (see Fig 9).In Fig. w tide level (h LW ) is approximately a linear function of h HW .

Table 2 .
Exceedance frequency of the design water level.