Numerical model set-up

Storm tides, tides, and waves are simulated using MIKE 21 FM HD, a coupled, depth averaged, hydrodynamic/wave model developed by the Danish Hydraulic Institute (see ref. 13). The model domain covers the entire North Sea and part of the adjacent North Atlantic and accounts for large scale meteorological and hydrodynamic effects13. In the coastal region, seabed topography at ~15 m resolution was obtained from the Schleswig-Holstein Agency for Coastal Defense, National Parks and Marine Conservation (LKN-SH). In the remaining domain, bathymetry at 30 arc-second intervals was obtained from the General Bathymetric Chart of the Oceans (GEBCO), which is produced by the British Oceanographic Data Centre (BODC). At the open boundary the model is forced with astronomical tides varying in time at a total number of 127 points along the domain. Tide levels have been calculated using the MIKE internal tide model. MSL rise at the boundary is forced using historical observations17 and RCP projections18,19 for present-day and future conditions by 2100 added on to the observed SLR, respectively. Waves and storm tides for all scenarios are produced by applying continuous wind and pressure fields from reanalysis data for the years 1970–201316. This time period was chosen due to the availability of high quality in-situ data but also to include the 1976 storm tide, the largest event ever recorded in many places of the German Bight. Model runs were output every 10 minutes. For every scenario, the 99.7th percentile water level exceedances are then estimated at ~500 m increments along the ~470 km North Sea coastline in Schleswig Holstein.

SLR projection

We consider the median of three different SLR projections by 2100 associated with the RCP4.518, RCP8.518, and RCP8.5HE19 (high end). All projections represent the effects of thermosteric and halosteric density changes, the response of the ocean to wind and pressure forcing, changes in ocean mass (Greenland and Antarctic ice sheets, glaciers, and groundwater), and glacial isostatic adjustment (see the Fifth Assessment Report of the IPCC8 for details). In addition, the RCP8.5HE projection includes rapid ice melt in the Antarctic, a plausible but more extreme sea level rise scenario29 that should nonetheless be considered from a coastal decision-making and management point of view. At a central point in the German Bight, these projections suggest a mean SLR of 0.54 m, 0.71 m, and 1.74 m by 2100 under RCP4.5, RCP8.5, and RCP8.5HE, respectively.

Extreme value statistics

Extreme value statistics are used to infer magnitudes of both storm tides and waves at specific ARI’s. We employ the Peak Over Threshold (POT) method4,30 and fit the following generalized Pareto distribution (GPD) to a ranked list of independent events exceeding a specified threshold of simulated high water peaks,

where c is the location (threshold) parameter, b is the scale parameter, a is the shape parameter and the threshold of exceedances is x. The parameters are estimated using the Maximum Likelihood Estimation (MLE) method4, with the threshold level of 99.7th percentile yielding consistent and stable results in the German Bight30. A declustering scheme based on the extremal index4 ensured that data were independent. Wave heights are described using the GPD but also a range of other common distribution functions including the Lognormal, Normal (Gaussian), Exponential, Weibull, and Generalized Extreme Value (GEV) distribution (see e.g. ref. 4). The best fitting distribution is assessed by calculating the minimum RMSE between the theoretical and empirical wave distributions.

The two univariate marginal distributions of storm tides and waves where then applied to assess their joint magnitudes and frequencies. We used Archimedean Copulas to describe the dependence between the two marginal distributions31, and hence the bivariate ARI’s. Specifically, for each SLR scenario and coastal grid point, we first obtain coincident samples of peak storm tides and wave heights in a window that is ± 120 minutes from the predicted high tide (Fig. 4s). The marginal distributions for storm tides (Fig. 2b) and waves (Fig. 2c) are obtained using univariate analysis. The dependence of storm tides and waves in our modelled data sets are then assessed using Kendall’s rank correlation. The correlation coefficient then becomes an input parameter in our copula analysis. Next, three types of copulas (Gumbel-, Clayton-, and Frank Copula) are evaluated and the model with the minimum RMSE between the parametric and the empirical copula32 is retained to estimate bivariate ARI’s. To qualitatively assess whether results are reasonable, 10,000 random events are generated from the parametric copula and the marginal distributions and compared to the numerical model data for consistency (see e.g. ref. 21 for an example).

Next, bivariate contour lines for the ARI 100 are calculated, resulting in a family of possible combinations which have the same recurrence interval. For example, in Fig. 2a, a small storm tide (300 cm) with large waves (185 cm) has the same historical ARI as a large storm tide (450 cm) with small waves (~100 cm; see HIS case). Wave run-up height (see next section) is then calculated for each ARI 100 event (elbow shaped contour line in Fig. 2a). The maximum overall height (i.e. storm tide plus wave run-up height) is assumed to be design relevant (see the red dot in Fig. 2a) and differences between the HIS and scenario runs indicate SLR induced changes in design heights.

Run-up

Dikes are constructed to withstand the impact of extreme water levels and waves. Potential dike failures result from several mechanisms, including overflow induced by elevated water levels and dike breaching caused by wave overtopping. In practice, dikes are built to withstand the wave run-up height R u,2% , the vertically measured distance which is exceeded by 2% of all incident waves20. Along the German North Sea coast, we assume all dikes consist of smooth embankments with a 1:6 slope (see Fig. 5a), following current recommended best practice. However, in reality the dike slopes slightly vary along different coastal stretches, and would need to be considered in a site-specific assessment.

Though extensively studied, uncertainty and bias is still found when wave run-up formulas are compared to physical test results (see e.g. ref. 20; and references therein). Here, we focus on relative wave run-up height changes from different SLR scenarios to minimize the effect of these uncertainties. Our assessment is based on a formula provided in ref. 20 describing the wave run-up on smooth and straight slopes (assuming all waves to attack perpendicular and in relatively deep water at the dike toe but without any wave breaking in front of the dike20), where the relative wave run-up height is calculated as

where H m0 is the wave height at the toe of the structure. The breaker parameter ξ m−1,0 is defined as

where α is the outer dike slope and the deep water wave length L 0 is given by