WAVE MODELLING IN A COMPLEX ESTUARY: STUDY IN PREPARATION OF FIELD MEASUREMENT CAMPAIGN EEMS-DOLLARD ESTUARY

This paper considers the Eems-Dollard estuary in the north of the Netherlands, which is part of the shallow shelf sea the Wadden Sea. This estuary is a highly complex area with tidal flats and deep channels and is characterised by an offshore directed wind, posing a challenge to wave prediction models. As little measurements are available, a measurement campaign is set up to provide field data for verifying and improving these wave models. The paper presents the locations that are most suited for wave measurements in the estuary, insight in the performance of the phase-averaged numerical wave model SWAN, and insight in the processes that play a role in the area around the corner of the Eemshaven. Furthermore, it presents insight into the reliability and applicability of SWAN in this highly complex area. An analysis of propagation effects is performed, as well as a comparison between the SWAN version as used for the Dutch dike safety assessment and the newest version, used for development, which includes the state of the art parameterisations of the physics. Furthermore, modelling is done for a schematised version of the area around the corner of the Eemshaven, considering several different model settings. Large differences occur in the results between the two SWAN versions. These differences are studied in more detail, and the causes of these differences are identified.


INTRODUCTION
Wave conditions at flood protection structures are important in the design of these structures. These conditions are often determined by using numerical wave models as measured data are not always available. It is of large importance that these wave conditions are accurately predicted by the applied model. Applying such a model in a highly complex estuary which is part of a shallow sea with tidal flats and wetlands, where the wind acts in an offshore direction, and in which no measurements are available, might mean reaching or even surpassing the model's limits of applicability.
The focus of this paper lies on the Eems-Dollard estuary in the north of the Netherlands, which is part of the Wadden Sea, a shallow shelf sea with several barrier islands, deep tidal channels, shallow tidal flats and wetlands, see Figure 1 and Figure 2. Multiple extensive studies, some including field measurements, have been performed before in the Wadden Sea, to gain a more detailed understanding of this complex area. Some noteworthy ones are Eslami Arab et al. (2012); Groeneweg et al. (2009); Van der Westhuysen et al. (2012); Groeneweg et al. (2015); Groeneweg & Van Vledder (2005); Van Dongeren et al. (2011). These studies mainly focused on the tidal channels and tidal deltas at the transition from the North Sea to the Wadden Sea. The Eems-Dollard estuary ( Figure 2) is an even more complex area. What makes the situation even more complex, is that almost no studies have been performed in the estuary and that almost no measurements are available inside the estuary. More specifically, the focus of this study lies on the dike section between the Eemshaven (the main port in the area) and the town of Delfzijl, as indicated by the red line in Figure 2. This dike protects a large part of the province of Groningen from flooding. The focus of this study is on the extreme conditions, as used for the dike safety assessment. These conditions correspond to return periods of 10,000 to 100,000 year.
Since the area is rather complex, and no measurements are available around the Eemshaven, it is difficult to determine the wave boundary conditions at this dike section. Furthermore, the (extreme) storms in the area tend to start out in the southwest and then turn to the northwest, but they never cross north to become eastward directed. This means, that at the dike, the wind is offshore directed during extreme conditions (as indicated by the arrow in Figure 2), which would indicate that wave heights tend to be low. Instead, relatively high waves are predicted for this location, probably due to the presence of a deep main channel in the area (indicated by the number 1 in Figure 2), and also a secondary channel that runs close to the dike (the 'Bocht van Watum', number 2 in Figure 2). The phase-averaged numerical wave model SWAN (Simulating WAves Nearshore, Booij et al. (1999)) is used for the dike safety assessment in the Netherlands. The version that is used (version 40.72ABCDE, the 'assessment version') predicts large, onshore directed waves at the dike, of order H s = 2 m. SWAN predicts that the waves turn around the corner of the Eemshaven and become onshore directed. However, since there are no measurements available, this means that the wave boundary conditions during (extreme) storms that are currently predicted by SWAN are without direct validation nearshore. This uncertainty in the wave conditions leads to uncertainties in the required crest level of the dikes and could potentially lead to millions of euros of unnecessary dike reinforcements. To the authors' knowledge, the performance of SWAN in predicting the wave conditions at the dike in such a complex area, consisting of deep tidal channels, shallow flats, and characterised by an offshore wind, has not been assessed before. Hence, the question is, whether using the model in this highly complex area means reaching the limitations of the underlying modelling concepts of SWAN, whether the predictions by SWAN in this area are accurate or not, what the causes are of this apparent turning of the waves around the corner of the Eemshaven, and whether the waves remain as large as the model predicts in reality as well.
This paper has three main goals. The first is to determine the locations that are most suited for wave measurements in the area, to be able to validate SWAN and gain more insight in the physics that play a role. For this first goal, an analysis of propagation effects is done. The second goal is to gain insight in the reasons why such large and onshore directed wave heights are predicted by SWAN at the dike and to gain insight in the processes that play a role in the area around the corner of the Eemshaven. The third is to assess the reliability and applicability of SWAN in this highly complex area. To reach these last two goals, a comparison between the SWAN version as used for the Dutch dike safety assessment (40.72ABCDE), and the newest version, used for development (41.10AB) and which includes the state of the art parameterisations of the physics, is performed. Furthermore, modelling is done for a schematised version of the area around the corner of the Eemshaven.
This study is made in preparation of a 12 year long field measurement campaign in the area, which aims to improve the reliability of the wave boundary conditions by improving the models, using a combination of numerical modelling and field measurements. Wind, water levels, currents, waves, wave run-up and wave overtopping will all be measured, by using a large range of different instruments, such as wave buoys, X-band radars, ADCPs, wave overtopping tanks and LIDAR.

LOCATIONS OF INTEREST
For setting up the first field measurements, the focus is on two locations, Uithuizerwad (UHW) and the Double Dike (DD, named after a dike reinforcement project in the area), see Figure 2. Uithuizerwad is located outside the estuary, to the west of the Eemshaven, where measurements of wind, water levels and waves are available at an existing measurement pole. The Double Dike is part of the dike section of interest, and is located inside the estuary south of the Eemshaven. Here, no measurements are available yet. Both locations are under oblique to very oblique wave attack, but at Uithuizerwad this wave attack is direct and coming from the Wadden Sea whereas at the Double Dike it is indirect. Furthermore, the wave conditions at Uithuizerwad are shallow water conditions, caused by a kilometre-long shallow area in front of the dike, located above mean sea level. The Double Dike lies sheltered from direct wave attack by the corner at the Eemshaven. The wave conditions there will not be depth limited, caused by a combination of a lower bottom level and higher water levels. The relatively deep (up to -10 m+NAP, Normaal Amsterdams Peil, Dutch ordnance level roughly corresponding to MSL) secondary channel Bocht van Watum runs relatively close to the Double Dike, and there is a shallower area located around mean sea level directly in front of the dike of several hundreds of metres wide.

Quantification of propagation effects
A quantification of propagation effects is made with SWAN, to determine the locations where waves should be measured to gain more insight in the physics that play a role in the area, as well as to be able to validate SWAN in more detail in a next step of this study.
The 40.72ABCDE version of SWAN is used for these stationary calculations with constant and uniform wind fields for storm conditions, corresponding to storms occurring on average once every 1, 10 and 100,000 year (the latter corresponding to the design conditions). These conditions are given in Table 1. The model settings as used for the Dutch safety assessment are applied, see Table 2. These settings were determined and validated in previous studies, e.g. Gautier (2010). The calculations use five telescoping nested curvilinear grids, with the one-way coupling from coarse to higher resolution grids by 2D-spectra. The grids are shown in Figure 3. No wave boundary conditions are applied on the outermost grid as local wave generation is sufficient for our purposes. This set-up is based on the set-up as used in Adema et al. ( , 2015. The DCSMv6 (Delft Continental Shelf Model version 6) model (Wenneker & Gautier, 2016) is split into three parts, an outer part (grid 1, blue), an inner part around the United Kingdom and Denmark (grid 2, red), and another refined part around the Dutch coast (grid 3, yellow). Furthermore, a refined version of the Kuststrook model (Wenneker & Gautier, 2016) is used (grid 4, magenta), and a sixteen times more detailed part of that model in the Eems-Dollard (grid 5, green), resulting in minimum grid cell sizes of approximately 10 by 20 m. The bathymetry is based on the bathymetry of the DCSM and Kuststrook models, as well as on the six-yearly bathymetrical survey data, and is shown on the left side of Figure 4.     QUAD DIA approximation  QUAD Bottom friction JONSWAP (Hasselmann et al., 1973) FRIC JON CON JONSWAP (Hasselmann et al., 1973) FRIC JON CON Depth-induced breaking Van der Westhuysen (2010) BREA WESTH alpha=0.96 pown=2.5 bref=-1.3963 shfac=500 Battjes & Janssen (1978); Salmon & Holthuijsen (2015) BREAK CON or BREAK BKD Triad wave-wave interactions LTA method (Eldeberky, 1996) TRIAD trfac=0.1 LTA method (Eldeberky, 1996)

Refraction
Turning rate cθ based on depth-gradients and using first-order backward differences New cθ implementation, based on gradient in wave number and using second order central differences (The SWAN team, 2017) Limiter on excessive directional turning based on Dietrich et al. (2013) Quadruplet and action limiter

Comparison SWAN assessment (40.72ABCDE) and development (41.10AB) versions
To determine why such large and onshore directed wave heights are predicted by SWAN at the dike and to gain more insight in the governing physics, a comparison is made between the safety assessment version (40.72ABCDE, cases A) and the development version (41.10AB, cases B) of SWAN. Furthermore, calculations are performed where the safety assessment version of SWAN is used for grids 1 to 4, and the development version for grid 5 (cases C). In this way, it is ensured that the differences between the models in and around the estuary are considered only. This comparison uses the same model grids, as well as the same storm conditions, as used in the quantification of propagation effects (Table 1). The (default) settings for SWAN version 41.10AB are also given in Table 2. Hence, each case (A, B or C) consists of three different calculations, one for each of the three storm conditions.

Modelling of area around the corner of the Eemshaven
To analyse the mechanisms and physics that play a role at the corner of the Eemshaven and in turning the waves around the corner in more detail, a schematised model of this is made, measuring 1.2 km by 2 km and using a rectilinear grid. Three versions of the model are used, where the first completely consists of a deep area, representing the main channel (-12.5 m+NAP). The second and third consist of a part of this deep area in the north and a shallow area, representing the area in front of the dike (0 m+NAP) to the south of this deep area. In between these two areas, running from west to east, there is a transitional slope of either 1:10 (second model) or 1:30 (third model), corresponding to the slopes as found in the area. The second model is shown on the right in Figure 4. The southern part of the western boundary (of the shallow area) is closed, representing the dike. On the north-western boundary and northern boundary, a JONSWAP spectrum is imposed, with conditions corresponding to the design conditions. These conditions are given in Table 3. Note that the whole model is rotated by approximately 45° in the counter-clockwise direction compared to reality. This means that the directions of the boundary conditions are rotated by 45° as well, see Table 3. The main focus of this analysis lies on gaining insight in the physics that play a role, e.g. by considering the 2D-spectra and source terms.

Quantification of propagation effects
The wave parameters at Uithuizerwad and the Double Dike resulting from our calculations with the assessment version correspond to the results of the Dutch dike safety assessment, with less than 2% difference in H m0 at Uithuizerwad and a maximum of 8% difference in H m0 at the Double Dike, which means that the results are largely representative of the Dutch safety assessment situation. Again large and onshore directed waves are found at the Double Dike during design conditions. The largest differences between the calculations performed here and the safety assessment results occur in the T m-1,0 wave period, with a maximum difference of 9% at Uithuizerwad and a maximum of 14% at the Double Dike. These differences are mainly caused by the different boundary conditions as used in the safety assessment, as well as the exclusion of the effect of flow fields in the present calculations.
The mean wave direction vectors for the 1 year condition are shown on the top left of Figure 5. For this case the waves become approximately alongshore directed at the Double Dike. The results for the 10 year (return period) and the design conditions are very similar, for these cases the waves become obliquely onshore incident. With the 1 year storm the mean wave direction turns less towards the dike, because of the more western wind direction. However, the trends are similar for all storm conditions, where the strongest turning of the mean wave direction occurs around the corner at the Eemshaven, at the shallower parts directly in front of the dike, and at the transitions from deep channels to shallow areas. From the absolute spatial gradient in the mean wave direction (in rad/m), which is plotted in the right panels of Figure 5, the locations of strongest wave turning can be confirmed. Once again, the behaviour is quite similar for all three storm conditions. The strongest turning occurs around the corner of the Eemshaven and at the transition from the secondary channel to the shallower area in front of the dike. There is also some wave turning further to the North, caused by a shallower area next to the main channel. Furthermore, there is an area more towards the east where turning occurs, which is caused by a deep inlet onto the large shallow flat. The same trends are found when the propagation terms (radiation stresses, SWAN energy propagation terms) are considered (not shown in the figures).
From these results, advice can be given on where to measure waves in the area to gain more insight in the physics affecting the waves, as well as to validate SWAN in more detail. As the strongest turning occurs around the corner of the Eemshaven and at the transition from the Bocht van Watum to the shallow areas in front of the dike, it is recommended to measure at these two locations, preferably with X-band radar, as well as with (arrays of) directional wave buoys.

Comparison SWAN assessment (40.72ABCDE) and development (41.10AB) versions
Large differences occur in the wave parameters between the two model versions, as shown in Figure 6 for the wave parameters during design conditions. An up to 30% lower H s is found at Uithuizerwad and the Double Dike, as well as an up to 20% lower T m-1,0 wave period when using the development version. However, these differences could already occur more offshore, and affect the waves already in the North Sea or Wadden Sea, and not necessarily in the estuary itself or close to the dike. The cases C, using the safety assessment version of SWAN for grids 1 to 4, and the development version for grid 5, provide insight into this. In this way, it is ensured that the local differences in and around the estuary are considered only. The wave conditions in large parts of the estuary are equal or almost equal for the cases with the development version on all grids (cases B) and the combination of the assessment version on grid 1-4 and the development version on grid 5 (cases C), see Figure 6. The main differences occur in the mean wave direction at location 3 to 6, where the mean wave directions between cases B and C differ by several degrees. The causes of these differences are unknown. Otherwise, the wave conditions are almost identical. As mentioned before, however, large differences occur between these cases and the cases where the assessment version was used for all five grids (cases A). The main difference being that larger low frequency waves are able to penetrate into the estuary with the assessment version. Low frequency waves are defined here as offshore generated North Sea wind waves with frequencies lower than 0.1 Hz.
The 2D-spectral boundary conditions for the boundaries of grid 5, resulting from grid 4, are largely similar for both the development and assessment versions of SWAN for the 1 and 10 year storms. However, for the design conditions, the boundary conditions are more extreme with the safety assessment version of SWAN. The northern boundary is located at relatively deep water in the North Sea. For case C with design conditions, the more extreme boundary conditions as found with the assessment version are applied to the boundaries of grid 5, which uses the development version.
The differences in the boundary conditions between cases B and C disappear quickly as the waves propagate into the domain, where these cases already show the same wave height before the waves reach the Wadden Sea. After this boundary area, cases B and C show largely the same wave conditions, which shows that despite the different boundary conditions, the wave conditions inside the Wadden Sea and Eems-Dollard estuary become largely equal when using the development version, and are thus largely determined by the model settings which influence the waves locally. This is confirmed by the very different wave conditions for case A, despite the use of the same boundary conditions as case C.
The 2D wave spectra during design conditions were considered at four locations for cases A, B and C. The locations are visible in Figure 7 and the spectra are shown in Figure 8 as polar spectra with contour lines at the levels 0.99, 0.9, ½, ¼, 1/8, 1/16, 1/32 and 1/64 times the local maximum energy density. In the North Sea, at location 1 in Figure 7 and Figure 8, the differences in wave heights and spectra are mainly caused by the differences in boundary conditions, with some differences in wind growth and whitecapping, which are somewhat more pronounced with the assessment version, and the depth-induced breaking, which is relatively somewhat stronger with the development version.

Figure 7. Ursell number (U=Hm0Lp 2 /d 3 ) during design conditions and the 2D spectra locations. Location 1 is located in the North Sea, location 2 lies on the tidal delta, location 3 is located just after the tidal delta in the main tidal channel, and location 4 is located further inside the main tidal channel.
The first differences between the two SWAN versions become visible at the transition from the deeper North Sea to the shallower Wadden Sea (location 2 in Figure 7 and Figure 8), where dissipation occurs. This reduction in wave energy is relatively stronger and occurs more rapidly for the development version than for the assessment version. This is caused by the differences in the parameterisations of the physics, mainly the depth-induced breaking. On the tidal delta, the difference between the two models in significant wave height is already around 1 m. The shape of the spectrum becomes different as well, with more westerly energy and less north-easterly low frequency energy in cases B and C.
When considering the 2D spectra from the assessment version around the transition from the North Sea to the Wadden Sea over the tidal delta (locations 2 and 3 in Figure 7 and Figure 8), largely the same behaviour is observed as in Van Nieuwkoop & Groeneweg (2014), where the assessment version of SWAN was compared to wave buoy measurements. Most of the wave energy travels over the ebbtidal delta, where the waves refract towards the shallow flats and where energy is dissipated, after which part of the wave energy enters the channel and refracts back. Here, it is also observed that inside the main channel, part of the low frequency wave energy comes from the western directions, which means that energy travels around the western side of the tidal delta, bending around the delta due to refraction, and then entering the channel. This westerly wave energy was not observed in the measurements in Van Nieuwkoop & Groeneweg (2014), but it was observed in their SWAN results. In the measurements, the directional distribution of wave energy had only one peak, with an even distribution on both sides. Van Nieuwkoop & Groeneweg (2014) present the non-linear wave interactions as described in Groeneweg et al. (2015) as an explanation for this. These non-linear interactions directionally broaden the spectrum, and are not included in SWAN. Similar effects were observed in Dusseljee et al. (2014). These nonlinear interactions likely play a role, as large values of the Ursell number (U=H m0 L p 2 /d 3 ) are found around the tidal delta and flats, indicating that nonlinear effects cannot be neglected (Figure 7).
Inside the deep tidal channel, the wave height is 1 to 1.5 m lower with case B (locations 3 and 4 in Figure 7 and Figure 8). There is more westerly energy than with the assessment version, a result of the waves travelling around the tidal delta in SWAN. Less of the original low frequency energy is left, which is caused by the differences in the refraction parameterisations.
In SWAN, several limiters are used for the calculation of the refraction. The assessment version of SWAN uses a quadratic frequency dependent refraction limiter (see Table 2), which limits the lowfrequency refraction for frequencies up to 0.2 Hz. This limiter was included because, under certain conditions, too little low frequency energy could penetrate the channels, compared to measurements (Van Vledder & Koop, 2009). This limiter, however, is not based on any known physics. In SWAN 40.72ABCDE, the turning rate c θ is based on depth-gradients and is calculated using first-order backward differences. This method of calculating the turning rate is rather inaccurate at coarse grids with steep bottom slopes. Therefore, a limiter to prevent excessive turning is necessary. In 40.72ABCDE, this limiter is based on limiting the bottom slope, similar to the one in WaveWatch III (Tolman, 2009). However, practical experience has shown that this limiter is not robust and effective enough (Dietrich et al., 2013;The SWAN team, 2017). The development version of SWAN does not use the frequency dependent refraction limiter of Van Vledder & Koop (2009). Furthermore, this version uses an updated version of the c θ -term (The SWAN team, 2017). In version 41.10AB, the turning rate is formulated in terms of the phase velocity, which yields more accurate results, even for coarse grids. Furthermore, the numerical implementation is more stable, as it uses central differences instead of the first order backward difference scheme as used in the assessment version. Hence, the use of a limiter to prevent excessive wave turning is less necessary with this method. However, with this method, an error can still occur for exceptionally large bottom slopes (The SWAN team, 2017). Therefore, a limiter is still included to prevent excessive turning. In this version, a more robust limiter is used, based on the CFL-criterion, which also prevents excessive frequency shifting (Dietrich et al., 2013).
The influence of the frequency dependent refraction limiter is visible in the 2D spectra in the deep tidal channel, where more low frequency energy is found, also for other directions, than with the development version. If this frequency dependent refraction limiter is applied with the development version, this same larger penetrative behaviour of the low frequency energy is found, but does not fully explain the differences between the two models.
Furthermore, in the tidal channel the spectrum is more symmetrical for cases B and C, with a more even distribution of wave energy on both sides of the peak direction. This high frequency shape of the spectrum corresponds better to the measurements as described in Van Nieuwkoop & Groeneweg (2014), but the larger amount of low frequency energy in case A (even with the refraction limiter which is not based on any known physics) corresponds better to the low frequency part of the spectrum in the measurements.
The apparent underestimation of low frequency energy in SWAN was also studied in Eslami Arab et al. (2012), where several hypotheses were tested. Wave driven currents due to wave breaking on a slope, nonlinear effects due to steepening of the waves on a slope, and too steep bathymetry gradients for SWAN to accurately model all could not explain the differences between SWAN and measurements. It was recommended to further study the wave reflection on steep slopes as a possible explanation.
Concluding, caution is necessary when applying SWAN in such a complex area, and caution is required when specifying the model settings. The largely different results between the two model versions shows that improvements in the understanding of the physics and the implementations of the physics in the model are still necessary. Therefore, the measurement campaign in the area is highly necessary, as well as more detailed modelling with e.g. SWASH and validation and application of the most recent developments in SWAN (ST6 physics (Rogers et al., 2012), Xnl quadruplet wave-wave interactions (Van Vledder, 2006)) in this area.
It is found that part of the differences in the wave conditions between the two SWAN versions are already caused offshore. Part of the differences are also caused locally, when the waves enter into the Eems-Dollard estuary around the corner of the Eemshaven, as well as when they travel onto the shallow flats in front of the dike. These effects are studied in the next section.

Modelling of area around the corner of the Eemshaven
To determine these local differences between the SWAN versions, and to gain more insight into the mechanisms and physics that play a role in the turning of the waves around the corner of the Eemshaven, the schematised model of this area was made. For the first model without a slope, using the safety assessment version of SWAN and during design conditions (so with an offshore directed wind), wave heights in the area close to the dike are approximately 0.8 to 1 m, see Figure 9. T m-1,0 wave periods lie between 2.7 and 4 s. This is much smaller than what was found from the large scale calculations in the previous section (H s around 2 m, T m-1,0 5-6 s). The mean wave direction lies around 335°N.
The locations for the case with a 1:10 slope are visible in Figure 9 and the 2D spectra for both the assessment version and development version are plotted in Figure 10. From the 2D spectra and the associated source terms, several effects can be observed. The further into the shadow zone, the less original (low frequency) wave energy is left. Only the most northern originating wave energy of the original (offshore) wave spectrum can reach this area, the rest is all filtered out by sheltering ( Figure  10). As visible in Figure 9, this sheltering of the waves causes wave height gradients to appear due to numerical diffusion, which resemble diffraction patterns, even when diffraction is not included in the SWAN model. The second effect consists of the quadruplet wave-wave interactions, which cause a secondary peak at twice the original peak frequency to appear in the shadow zone, as well as spreading of energy to different directions at the higher frequencies, both to the more northern and western directions ( Figure  10). The third effect is the wind growth, which causes the higher frequency peak to grow and causes energy to appear at a broader range of directions, spread around the wind direction in a ±90° sector ( Figure 10). These effects of sheltering of wave directions and wave-wave interactions cause the mean wave direction to turn more towards the north, and cause part of the wave energy to become onshore directed, even with an offshore wind direction and without diffraction in the model. For the case without any slope, the bottom friction, triads and surf-breaking do not play a role, because of the relatively deep water conditions everywhere.
In the cases with a 1:10 or 1:30 slope, there is somewhat more whitecapping at the original peak frequency at locations around to the slope, caused by steepening of the waves. Furthermore, the bottom friction is stronger, caused by the slope and shallow area. Around the slope, the triad interactions play a minor role, but on the shallow flat they disappear again. Even for the cases with a slope, the refraction only plays a minor role. The maximum turning that occurs on the slope is only a few degrees, and this refraction causes the somewhat more westerly wave energy to be able to reach the shadow zone as well. The combined positive effects of more energy reaching the shadow zone caused by refraction and negative effects of the whitecapping and bottom friction result in wave heights in the shadow zone that are several cm larger than for the case without a slope, a slightly larger wave period (0.2 s difference) and a slightly more northern mean wave direction (1-2° difference). Furthermore, the triad wave-wave interactions around the slope lead to a slightly more pronounced secondary high frequency peak in the wave spectra.
The same trends hold for locations more offshore (eastward), further from the shadow zone. The further offshore, the more of the original wave energy can reach these areas, resulting in larger wave heights, due to which also the triad interaction, bottom friction and whitecapping play a larger role, and where also depth-induced breaking starts to occur. This leads to larger differences more and more offshore (towards the east) between the cases with and without a slope. However, outside the region of the slope, the differences between the cases with and without slope quickly decrease again to become very small. The differences between the 1:10 and 1:30 slope cases are minimal as well.
Thus, overall the influence on the boundary conditions at the dike of a west-east slope of 1:10 or 1:30 and a shallow area at 0 m+MSL are minimal. For all these cases the Ursell numbers are low as well (<1), which indicates little influence of nonlinear effects.
Using the development version of SWAN leads to similar wave heights offshore, but lower wave heights in the shadow zone (up to 15 cm or up to 20% lower), see Figure 10. Compared to the assessment version, there is now slightly more (original) low frequency energy (f < 0.2 Hz) present in the shadow zone, which leads to a somewhat larger main spectral peak and wave period (0.3 to 0.5 s or up to 10% larger). The quadruplet interactions as schematised by the Discrete Interaction Approximation  add less energy to the secondary higher frequency peak, and add energy to higher frequencies than with the assessment version. Both the westerly and northerly high frequency tails of the spectrum are less pronounced, which combined with the larger low frequency peak leads to a slightly more northern mean wave direction (1°-3° difference).
When comparing the development version without and with a slope, the behaviour corresponds to what was found for the assessment version. A somewhat larger H m0 of up to 5 cm is found in the shadow zone for the cases with a slope, as well as a slightly larger T m-1,0 period (up to 0.4 s) and more northerly mean wave direction (3-4° difference). All these effects are caused by the refraction, which allows more low frequency energy to reach the area.
Thus, part of the differences in the wave conditions between the two SWAN versions were already caused offshore, but part of the differences are also caused locally, when the waves enter into the Eems-Dollard estuary around the corner of the Eemshaven and when they travel onto the shallow flats in front of the dike. These local differences are again caused by the different parameterisations of the physics between both models.

DISCUSSION
During the calculations, it was observed that SWAN has difficulties converging in this highly complex area, and that even after the maximum amount of 80 iterations as set in the Dutch safety assessment settings, varying amounts of computational points had not yet converged. This mainly influences the mean wave direction and directional spreading, as they converge slower than e.g. the wave height and wave period. Therefore, a stricter stop criterion was applied, where the calculation is stopped when the relative change in H s from one iteration to the next is less than 0.003and the curvature of the iteration curve of H s normalised with H s is less than 0.002 in 99.5% of the wet grid points. The maximum amount of iterations was set to 160. This new stop criterion improved the convergence sufficiently.
In the calculations comparing the two versions of SWAN in the Wadden Sea and Eems-Dollard estuary, another effect was observed in the results of the assessment version. Spatial alongshore oscillations with a length scale of 1 to 2 kilometre occur in the wave parameters along the dike around Uithuizerwad, mainly in the T m-1,0 (up to 10%), but also in the H m0 (up to 5%). These oscillations were found in the actual Dutch safety assessment results as well. To find the cause, a schematised model was made of an inlet between two islands in the Wadden Sea. This effect turns out to be related to the sweep mechanism in SWAN that is used for curvilinear grids. With rectilinear/regular grids, this effect does not play a role, and therefore is something to keep in mind when deciding on the use of (curvilinear) grids in these complex areas.
As mentioned before, the wave heights as found with the schematised model were much lower than with the large scale model. This could e.g. be due to errors in the rotation of the schematised model compared to reality or due to the simplifications in the bathymetry. In a next step of this study, the channel Bocht van Watum which runs close to the dike will be included in the models as well, which could also be a cause of the lower wave heights as found in the schematised model.
The assessment version of SWAN uses 36 directional bins of 10°. This coarse directional resolution might lead to removal of certain wave directions by sheltering which in reality would still be able to penetrate. This was tested by using 180 directional bins of 2° each, which indeed gives a smoother behaviour of the wave parameters and wave spectra.
By default, SWAN uses the SORDUP numerical scheme for the wave propagation. In case the BSBT scheme is used, which shows more numerical diffusion, the wave height gradients are smoothed much more, resulting in larger wave heights at the dike. This again shows that the model settings need to be carefully picked in a highly complex area like the one considered here.
Diffraction effects may play a role around the corner of the Eemshaven. In SWAN, a phasedecoupled refraction-diffraction approximation by Holthuijsen et al (2003) can be used. The approximation is expressed in terms of the directional turning rate of the individual wave components in the 2D wave spectrum and is based on the mild-slope equation for refraction and diffraction, omitting phase information. Using the approximation leads to slightly smoother wave height patterns, but does not lead to a largely different wave height at the dike. However, this does not mean that diffraction does not play a role in reality, since the phase information is missing. Calculations with a numerical model like SWASH, which is able to fully determine the diffraction effects, should be used to determine if diffraction effects play a role in reality. Furthermore, this again is an argument that the field measurements in the area are necessary.

CONCLUSIONS
This paper determined the locations that are most suited for wave measurements in the Eems-Dollard estuary, to gain insight in the reasons why large onshore directed wave heights are predicted by SWAN with offshore directed winds. It also provided information on the processes that play a role in the area around the corner of the Eemshaven and at the tidal inlet, and an assessment was made of the reliability and applicability of SWAN in this highly complex area.
The area considered here is highly complex, and no data for validation is available yet. This is highly necessary and will be filled in by the upcoming field measurements. Caution is necessary when applying SWAN in such a complex area, as the complex bathymetry and offshore wind conditions cause the underlying modelling concepts of SWAN to reach their limits. Furthermore, caution is necessary when picking both the settings for the physics and the numerics, because slightly different settings can lead to very different results. One has to be aware of the limits of applicability of the model, in conjunction with the assumptions within the model. The preferred measurement locations were determined by considering the propagation effects. These preferred locations are the corner of the Eemshaven and the transition area from the deep channel Bocht van Watum towards the shallow areas in front of the dike.
Large differences occurred in the results between the assessment and development SWAN versions. These differences appeared both offshore, at the transition from North Sea to Wadden Sea, as well as more locally inside the estuary, around the Eemshaven and close to the dike. Offshore, these differences were mainly caused by the differences in the parameterisations of the refraction. The lack of agreement between both SWAN versions and previously studied field measurements shows that improvements in the understanding of the physics and improvements in the model parameterisations are necessary. Inside the estuary, several main effects in turning the waves around the corner were observed: wave sheltering, wave-wave interactions, the wind growth and numerical diffusion effects. The nearshore differences between both models were caused by differences in the parameterisations of these physical processes, as well as by the different refraction parameterisations. It is likely that locally diffraction effects play a role as well, due to the wave focusing effect in areas with strong refraction such as around the corner of the Eemshaven. Whether the waves become onshore directed in reality as well, has to be determined by further study. This can be done e.g. by the field measurements, by including the Bocht van Watum in the SWAN models, and by more detailed numerical modelling using SWASH. All these aspects will be considered in the next steps of the study.