Sensitivity analysis of numerical aspects of SWAN: Activity 3.1 of SBW project Waddensea

(1)

(2)

Sensitivity analysis of numerical aspects

of SWAN

Gerbrant van Vledder, Marcel Zijlema, Erick Rogers and Jacco Groeneweg

Report May, 2007

(3)

TITLE: Sensitivity study of numerical aspects of SWAN

ABSTRACT:

The spectral wave model SWAN plays a key role in the estimation of the Hydraulic Boundary Conditions (HBC) for the primary sea defenses of the Netherlands. The SWAN model is used to transform offshore wind and wave conditions to wave conditions near the dikes. It requires specialized knowledge to study the possible shortcomings inSWAN related to the wave boundary conditions in the Waddensea. The shortcomings are both physical and numerical. The physical shortcomings of SWAN will be investigated in a future subproject. In the present subproject the focus is on the numerics. The objective of the present subproject is fourfold. Firstly, given the presently used grid and bottom schematizations as described in WL (2006b), to investigate the errors due to grid schematization and numerical discretizations. Secondly, to investigate the convergence behavior in the Waddensea. Thirdly, to investigate the required spectral resolution. Here, attention is paid to the frequency range regarding initial wave growth and the required directional resolution for the penetration of swell waves into the Waddensea. Fourthly, to investigate the requirements to curvi-linear grids from the viewpoint of wave modelling. In this way, we gained insight into the quality of the employed grid and bottom schematizations with SWAN.

The main conclusions from this study are:

Concerning the previous hindcast study by WL (2006b), grid resolutions of 100 m seem to be sufficiently fine to obtain acceptably small discretization errors;

The default convergence criterion is not robust enough and should be replaced by the curvature-based criterion; The directional resolution affects the penetration of swell into the Waddensea, but no convergence was obtained with respect to the directional resolution;

Curvi-linear grids for wave model applications should be denser in areas where physical processes cause strong variations in the wave conditions;

REFERENCES: RKZ-1697 SAP bestelnummer “4500041303 pos 20”

VER AUTHOR DATE REMARKS REVIEW APPROVED BY

1 J. Groeneweg et al. 22 November 2006 draft A.R. van Dongeren M.R.A. van Gent 2 G.Ph. van Vledder et al. 21 March 2007 revised draft A.R. van Dongeren M.R.A. van Gent 3 G.Ph. van Vledder, M. 29 May 2007 final A.R. van Dongeren M.R.A. van Gent

Zijlema, W.E. Rogers, J. Groeneweg

PROJECT NUMBER: H4803.60

KEYWORDS: Numerical analysis, SWAN, SBW, hydraulic boundary conditions

NUMBER OF PAGES: 40

CONFIDENTIAL: YES NO

(4)

Amelander Zeegat... 2—1 2.1 Introduction ... 2—1 2.2 Estimate of global error due to (geographic) grid discretization ... 2—1 2.2.1 Global errors of wave parameters computed for area1 ... 2—2 2.2.2 Global errors of wave parameters computed for area2 ... 2—3 2.3 Accuracy due to spectral grid resolution ... 2—4 2.4 Summary ... 2—7

3 Convergence behavior in Waddensea computations ... 3—1

3.1 Introduction ... 3—1 3.2 Results with original convergence criterion ... 3—1 3.3 Results with the curvature-based convergence criterion ... 3—2 3.4 Summary ... 3—3

4 Requirements to the spectral resolution... 4—1

(5)

4.3 The effect of directional resolution on swell propagation ... 4—4 4.3.1 Model setup... 4—4 4.3.2 Results ... 4—5 4.3.3 Convergence behavior ... 4—7 4.4 Summary ... 4—8

5 Requirements to curvi-linear grids ... 5—1

5.1 Introduction ... 5—1 5.2 Experience with previously used non-uniform grids ... 5—2 5.3 Summary of experience with non-uniform grids ... 5—3 5.4 Method of analysis ... 5—4 5.5 Effect of grid resolution ... 5—5 5.5.1 Convergence behavior ... 5—6 5.6 Aspect ratio... 5—6 5.7 Skewness ... 5—7 5.8 Requirements ... 5—8 5.9 Summary ... 5—8

6 Comparison of a hindcast on a curvi-linear and a rectangular grid ... 6—1

6.1 Introduction ... 6—1 6.2 Model setup ... 6—1 6.3 Dedicated curvi-linear grid for the Amelander Zeegat... 6—2 6.4 Analysis of results... 6—3 6.5 Summary ... 6—4

7 Conclusions and recommendations ... 7—1

(6)

List of Figures

2.1 Buoy locations in Amelander Zeegat in area1 , area2, area3 and area4. 2.2 Error estimates of wave parameters in area1 of Amelander Zeegat. 2.3 Error estimates of wave parameters in area2 of Amelander Zeegat.

3.1 Convergence behavior of Hs and Tm01 in Amelander Zeegat buoys obtained at grid

‘grid1’ with a mesh size of 100m.

‘grid2’ with a mesh size of 20m and constant water level and no current.

‘grid2’ with a mesh size of 20m and space-varying water level and current.

4.1 Outline of curvi-linear grid of the Dutch coast for the Waddensea area and outline of grid AZG1.

4.2 Bathymetry of Waddensea on the curvi-linear grid GridCL.

4.3 Variation of the significant wave height Hm0 and the spectral period Tm-1,0 for grid

AZG2 and for the storm event of 8 Febr. 2004 using an upper frequency of 0.85 Hz. 4.4 Variation of the spectral periods Tm01 and Tm02 for grid AZG2 and for the storm event

of 8 Febr. 2004 using an upper frequency of 0.85 Hz.

4.5 Difference in significant wave height Hm0 and spectral period Tm-1,0 for grid AZG2

and for the storm event of 8 Febr. 2004 based onSWAN computations with upper frequencies of 0.85 Hz and 1.25 Hz.

4.6 Difference in spectral periods Tm01 and Tm02 for grid AZG2 and for the storm event of

8 Febr. 2004 based onSWAN computations with upper frequencies of 0.85 Hz and 1.25 Hz.

4.7 Difference in significant wave height Hm0 and spectral period Tm-1,0 for grid AZG2

and for the storm event of 8 Febr. 2004 based onSWAN computations with upper frequencies of 0.85 Hz and 1.65 Hz.

4.8 Difference in spectral periods Tm01 and Tm02 for grid AZG2 and for the storm event of

8 Febr. 2004 based onSWAN computations with upper frequencies of 0.85 Hz and 1.65 Hz.

4.9 Variation of the significant wave height Hm0 and spectral periods Tm-1,0, Tm01 and Tm02

along the main ray in the Amelander Zeegat. Results ofSWAN computations with upper frequencies of 0.85 Hz and 1.25 Hz. Storm event of 8 Febr. 2004.

4.10 Variation of the significant wave height Hm0 and spectral periods Tm-1,0, Tm01 and Tm02

along the main ray in the Amelander Zeegat. Results ofSWAN computations with upper frequencies of 0.85 Hz and 1.25 Hz. Storm event of 8 Febr. 2004.

4.11 Spatial variation of significant wave height Hm0 and spectral period Tm-1,0 computed

on grid AZG3A using a directional resolution of 10°.

4.15 Relative difference of significant wave height Hm0 and spectral period Tm-1,0 between

the computations on grid AZG3A with directional resolutions of 2° and 10°.

(7)

the computations on grid AZG3A with directional resolutions of 2° and 4°.

4.18 Variation of integral wave parameters Hm0, Tm-1,0, Tm01, and in output point 6 as a

function of the directional resolution . Based on computations on grid AZG3A and a water level of NAP +3 m.

4.24 Location of cross sections trough curvi-linear grid AZG3A

4.25 Variation of depth and integral wave parameters along cross section 6 in grid AZG3A as a function of the directional resolution . Water level NAP +3 m and a directional distribution of 10°.

4.28 Convergence behavior of integral wave parameters for test location 14 in the Amelander Zeegat. Based on computations on grid AZG3A, a water level of NAP+3m and a directional resolution of 10°.

5.1 Bottom profile for quasi-1D-computations

5.2 Detail of curvi-linear Grid5 with a spatial resolution varying from 200 m to 50 m near the bottom slope, and a constant spacing of 50 m in y-direction.

5.3 Variation of significant wave height Hm0, spectral periods Tm-1,0, Tm01, mean wave

direction and directional spreading along the 1D transect based on computations on the quasi-1D grids 1, 2, 3 and 4.

5.4 Variation of significant wave height Hm0, spectral periods Tm-1,0, Tm01, mean wave

(8)

5.5 Convergence behavior of integral wave parameters on the quasi-1D grid 1 at output location 1.

5.6 Convergence behavior of integral wave parameters on the quasi-1D grid 1 at output location 6.

5.7 Outline and grid properties of quasi-1D grid 5, with y=50 m.

5.8 Comparison of integral wave parameters along 1D-transect at y=10km, computed on the quasi-1D grids 5 and 10, with a y=200 m and 1,000 m, respectively.

5.9 Comparison of integral wave parameters along 1D-transect at y=10km, computed on the quasi-1D grids 5 and 10, with a y=200 m and 10,000 m, respectively.

5.10 Outline, aspect ratio and depth of curved grid 25 with a maximum skewness of 8.9°. 5.11 Comparison of integral wave parameters along 1D-transect at y=10km, computed on

the curved quasi-1D grids 21 and 25, with a shift in x-direction of 0 m and 1,000 m, respectively.

5.12 Comparison of integral wave parameters along 1D-transect at y=10km, computed on the curved quasi-1D grids 21 and 26, with a shift in x-direction of 0 m and 2,500 m, respectively.

5.13 Variation of the absolute value of the gradient in significant wave height Hm0, and the

spectral period Tm-1,0 and Tm01 along the 1D-transect at y=10 km computed on the

quasi-1D grids 3 and 5.

6.1 Outline of and average cell length of curvi-linear grid AZG1. Grid lines are shown every line.

6.2 Aspect ratio and skewness of curvi-linear grid AZG1.

6.3 Spatial variation of significant wave height Hm0 and spectral period Tm-1,0 for the

storm event of 8 Febr. 2004, computed on grid AZG2A.

6.4 Spatial variation of the absolute values of the gradient in significant wave height Hm0

and spectral period Tm-1,0 computed on grid AZG2A for the storm event of 8 Febr.

2004.

6.5 Outline of and average cell length of curvi-linear grid AZG1A. Grid lines are shown every second line.

6.6 Aspect ratio and skewness of curvi-linear grid AZG1A.

6.7 Outline of and average cell length of curvi-linear grid AZG2. Grid lines are shown every second line.

6.8 Outline of and average cell length of curvi-linear grid AZG2A. Grid lines are shown every fourth line.

6.9 Outline of and average cell length of curvi-linear grid AZG3. Grid lines are shown every third line.

6.10 Outline of and average cell length of curvi-linear grid AZG3A. Grid lines are shown every fifth line.

6.11 Difference in significant wave height Hm0 and spectral period Tm-1,0 between

computations on the grids AZG1 and AZG1A for the storm event of 8 Febr. 2004. 6.12 Comparison of integral wave along the output ray in the Amelander Zeegat based on

computations on the grids AZG1 and AZG1A for the storm event of 8 Febr. 2004. 6.13 Difference in significant wave height Hm0 and spectral period Tm-1,0 between

(9)

6.18 Comparison of integral wave along the output ray in the Amelander Zeegat based on computations on the grids AZG3 and AZG2A for the storm event of 8 Febr. 2004. 6.19 Difference in significant wave height Hm0 and spectral period Tm-1,0 between

(10)

List of Tables

2.1 Mesh differences in Amelander Zeegat buoys in area1 for h = 100m and 200m. 2.2 Mesh differences in Amelander Zeegat buoys in area2 for h = 20m and 40m. 2.3 Computed (maximal) Courant numbers in Amelander Zeegat buoys.

3.1 Number of iterations obtained with different stopping criteria. 4.1 Summary of boundary conditions for the two storm events.

(11)

1 Introduction

1.1 The hydraulic boundary conditions

In compliance with the Flood Defenses Act (“Wet op de Waterkering, 1996”), the primary coastal structures must be checked every five years (2001, 2006, 2011 etc.) for the required level of protection on the basis of the Hydraulic Boundary Conditions (HBC) and the Safety Testing Regulation (VTV: Voorschrift op Toetsen op Veiligheid). TheseHBC must be derived anew every five years and established by the Minister of Transport, Public Works and Water Management.

There is a degree of uncertainty concerning the quality of the current HBC, in particular those for the Waddensea. This is because they were obtained from an inconsistent set of measurements and design values (WL, 2002a). For the rest of the Dutch coast (the closed Holland Coast and the Zeeland Delta) the wave transformation modelSWAN (Booij et al., 1999,SWAN homepage http://www.swan.tudelft.nl) was applied.

Presently, there is insufficient wave information for validation and calibration of SWAN in the Waddensea. Consequently, it is not yet possible to obtain reliable and validated HBC with this model in a tidal inlet system such as the Waddensea. Furthermore, it is uncertain to what extent swell waves penetrate into the Waddensea. Measurements near the Emmapolder in Groningen have shown that swell provides a considerable contribution to wave run-up, in the order of 30% (Seijffert, 1991).

1.2 SBW Boundary Conditions Waddensea – Project H4803

The abovementioned problem was the direct cause for the request from the subproject “Boundary Conditions”, which is part of the main project “Strength and Loading of Coastal Structures (SBW: Sterkte en Belasting Waterkeringen)”, toWL| Delft Hydraulics to formulate a Plan of Action in which a strategy would be determined to answer the principle question: “How do we arrive at reliable Hydraulic Boundary Conditions for the Waddensea for 2011?” In addition to the penetration aspect, the general suitability of theSWAN wave model in the Waddensea must be determined and the improvements required to produce reliable HBC in the Waddensea should be specified.

This Plan of Action (WL, 2006a) is the responsibility within theSBW project (SBW, 2005), because it has the task of improving the quality of the models and methods used to derive theHBC to enable the managers and experts to have sufficient confidence to use these tools for the five-yearly tests. This is achieved by performing activities such as carrying out measurements and developing methods and models to be able to transform the deep-water statistics into statistics in front of coastal structures and to consider the higher loads on, and the failure mechanisms of, these defenses.

(12)

is used to derive the Hydraulic Boundary Conditions and focuses on improving the predictions of the wave boundary conditions produced by the wave model that are required to make the transformation from deep water to coastal structures.

In July 2006 RIKZ commissioned WL | Delft Hydraulics for the project “Execution of the Plan of Action SBW Hydraulic Boundary Conditions” (in Dutch: “Uitvoering Plan van Aanpak SBW RVW Waddensea”). Ap van Dongeren is projectleader of this project (H4803).

The subproject “Sensitivity analysis on numerical aspects of SWAN” is a subproject of the aforementioned project. The project team of this subproject consists of dr. ir. M. Zijlema (Delft University of Technology), dr. ir. G. Ph. van Vledder (Alkyon Hydraulic Consultancy & Research), W. E. Rogers (Naval Research Laboratory, Stennis Space Center) and dr. J. Groeneweg (WL | Delft Hydraulics). Marcel Zijlema and Gerbrant van Vledder carried out different parts of the numerical study. Rogers performed the external review. Jacco Groeneweg was projectleader.

1.3 Sensitivity analysis on numerical aspects – Subproject

H4803.60

The SBW project (SBW, 2005) has the task of improving the quality of the models and methods used to derive the HBC. The spectral wave model SWAN plays a key role in the estimation of the Hydraulic Boundary Conditions (HBC) for the primary sea defenses of the Netherlands. The SWAN model is used to transform offshore wave conditions to wave conditions near the dikes. It requires specialized knowledge to study the possible shortcomings inSWAN related to the wave boundary conditions in the Waddensea. Various shortcomings regarding the parameterization of physical processes and numerical methods inSWAN were revealed in previous studies, such as WL/Alkyon (2002);WL/Alkyon (2004). These were summarized in Alkyon (2005) insofar as they apply to the Netherlands inland waterways. However, these physical and numerical shortcomings also apply to the Waddensea, in particular the wind input modelling.

The physical shortcomings of SWAN will be investigated in a future subproject. In the present subproject the focus is on the numerics. A sensitivity analysis of numerical aspects inSWAN was carried out with respect to the Waddensea, according to the activities defined in the Plan of Action (WL, 2006a). This includes:

1. Discretizations in geographic space: Error estimates due to geographic discretizations. 2. Convergence behavior: Study of the convergence behavior and consideration of the

default convergence criterion. 3. Grid schematization:

a) Spectral resolution: choice of resolution in directional space and frequency space, as well as cut-off frequencies and spectral tail.

b) Curvilinear versus rectilinear with a critical consideration of the presently used curvilinear grids that are designed for flow computations.

(13)

schematization and numerical discretizations. Secondly, to investigate the convergence behavior in the Waddensea. Thirdly, to investigate the required spectral resolution. Here, attention is paid to the frequency range regarding initial wave growth and the required directional resolution for the penetration of swell waves into the Waddensea. Fourthly, to investigate the requirements to curvi-linear grids from the viewpoint of wave modelling. In this way, insight was gained into the quality of the employed grid and bottom schematizations withSWAN.

To limit the scope of work of this study, the effect of currents on the requirements to the computational grids has not been investigated, except when the global errors of the spectral grid resolution (current refraction) were determined. However, neglecting current effects does not imply that they are not important.

In this report computation were made for the tidal inlet of Ameland between the islands of Terschelling and Ameland. This inlet is also used for hindcast studies with SWAN and comparison against buoy data collected along a ray in this inlet. It is noted that measured data are not used in this study.

(14)

2 Estimate of errors in wave parameters due

to grid resolutions of Amelander Zeegat

2.1 Introduction

To verify theSWAN model results of e.g., Amelander Zeegat, it is important to get insight in the accuracy of the computations. Errors can be obtained at different levels or arise from numerical computations. Examples are use of lower- or higher-order difference schemes, iteration convergence, non-orthogonality or non-smoothness of curvilinear grids, drying and flooding, assumptions in formulations of physical processes, etc.

In Sections 2.2 and 2.3 we estimate global errors in wave parameters Hs, Tm 1,0 en Tp due to

the use of

geographic grid discretization only and spectral grid discretization,

respectively, of the model of the Amelander Zeegat.

These discretizations are related to the energy transport of the action balance equation. Since source terms can be evaluated directly, numerical accuracy is not relevant for these terms, except for errors due to interpolations. In Witteveen and Bos (2004), it has been shown that interpolation errors are insignificant compared to the discretization errors. Hence, with respect to the source terms, it amounts to the quantification of the physical accuracy which is beyond the scope of this subproject.

Furthermore, it must be stressed that the value of discretization errors of the energy transport at relatively small scales in geographical space can be ignored when the source terms are included, since they are dominant (Tolman, 1991). Hence, it is expected that the effect of the source terms on the accuracy of energy transport in geographical space is negligible and therefore are not included in this sensitivity analysis.

Section 2.4 concludes some remarks with respect to the hindcasts of Amelander Zeegat.

2.2 Estimate of global error due to (geographic) grid

discretization

The global error is the error occurring in the numerical approximation. More exactly: is the exact solution of the continuous modelproblem and h is the numerical approximation

obtained with a discretized model with mesh size h. The global error is defined as

(15)

To get insight in the global error, we make computations with three different mesh sizes: h, 2h and 4h. Next, we calculate the following mesh differences in the solution of subsequent mesh sizes:

2

2 4 2

h h h

By means of the so-called Richardson correction we may obtain an estimate of order p for the global error. The estimate of the global error is

2 2 1 2 1 / / h p h h Thus, we have 2 log h log 2 h p

The order of accuracy of the numerical approximation measures the ability of the approximation to reproduce the transport terms of the action balance equation. In general, the higher the order of accuracy of an approximation, the better it is able to reproduce the transport terms.

The calculations of the Amelander Zeegat are carried out with four (nested) areas, see Figure 2.1. The basis grid used for the largest area (area1) has a mesh size of 100 m, while for the other nested areas (area2, area3 and area4) a basis mesh size of 20 m is employed. Further details on these grids can be found in WL (2006b, Section 3.4.1 and Figure 3.25). Here we repeat the calculations twice for area1 with mesh sizes 200 m and 400 m and for

area2 with sizes 40 m and 80 m.

The estimates of the global errors are made for the significant wave height Hs and periods

Tm 1,0 and Tps. Note that the so-called smoothed peak period Tps is computed more accurately

than the standard peak period Tp by means of parabolic fitting the largest three energy

density bins of the spectrum after which the maximum of the fitting is determined.

2.2.1 Global errors of wave parameters computed for area1

We consider the storm situation of January 2, 2005 at 10:00 hr. Details on model settings can be found in the WL (2006b, Section 3.4). At the incoming boundary the wave height, mean period Tm01, mean wave direction and directional spreading are given by 5.3m, 7.2s,

(16)

Figure 2.2 depicts the global errors _h and ₂_h of wave parameters Hs, Tm 1,0 and Tps at

domain area1 with basis mesh size of h = 100 m. The largest differences due to grid discretization are found between the islands Terschelling and Ameland and between the islands and the mainland. The differences become smaller as the mesh sizes become smaller. At h = 100 m, the largest difference in Hs is in order of 5 10 cm whereas the largest

difference in the period measures is about 0.2 0.5 s.

Looking more specifically, the mesh differences of Hs, Tm 10 and Tps in some wave buoys in

Amelander Zeegat have been determined (see Figure 2.1 for the buoy locations). These are displayed in Table 2.1. As expected, the smallest differences can be found near the outer boundary (AZB11 and AZB12) and are in the order of 0.1 mm and 1.0 ms for significant wave height and wave periods, respectively. The largest differences are to be near and behind the islands (AZB41, AZB42, AZB51 and AZB52) and are in the order of 10 cm and 10 ms for wave height and periods, respectively. These differences are probably due to the representation of the islands on a given grid which is different for each mesh size. Note that these differences are thus not due to discretization errors. The very large differences in AZB41 and AZB42 of the order 30-40 cm are probably due to nesting since they are close to the boundary of the nested area area2. On average, the mesh difference in Hs for mesh

size h=100 m is about 10 mm, whereas that of period measures is roughly 10 ms.

Although there is no systematic behavior in the error as function of the mesh sizes, it appears that roughly p=1. This is expected, since the local order of the applied scheme SORDUP is 2. For small coastal regions, a global error of order 1 is good enough.

h = 100m h = 200m (Hs) 105 m (Tm 10) 103 s (Tps) 104 s (Hs) 105 m (Tm 10) 103 s (Tps) 104 s AZB11 9.0 4.2 15.9 236.0 0.6 18.1 AZB12 51.0 0.3 6.0 218.0 3.3 16.0 AZB21 1074.0 12.6 54.9 443.0 9.9 40.1 AZB22 69.0 4.8 11.1 3203.0 36.1 168.9 AZB31 822.0 5.1 17.1 3940.0 12.2 42.9 AZB32 2403.0 5.4 9.9 3962.0 2.4 26.9 AZB41 38067.0 25.5 66.0 37761.7 18.7 269.0 AZB42 30117.0 14.7 6.9 36898.0 47.4 436.9 AZB51 1358.1 12.0 51.0 2853.1 90.4 219.0 AZB52 3180.9 3.9 158.1 3363.1 65.0 267.9

Table 2.1. Mesh differences in Amelander Zeegat buoys in area1 for h = 100m and 200m.

2.2.2 Global errors of wave parameters computed for area2

The results obtained with area2 are shown in Figure 2.3. Compared to the results of area1 (cf. Figure 2.2) the region with largest differences is significantly smaller. Moreover, the mesh differences are also smaller. Globally, the mesh difference in Hs for h = 20 m is within

5 cm, the error in Tm 1,0 is within 0.1 s and the error in Tps is far within 0.5 s. They are

(17)

different processes (generation, dissipation and interactions) used in SWAN are relative larger.Again, differences at the edge of the islands are due to the representation of these islands on different grids.

In Table 2.2, mesh differences of wave parameters in the wave buoys AZB21, AZB22, AZB31 and AZB32 (all in area2) are shown. Typically, the mesh difference for h = 20 m in significant wave height obtained at average 2 mm and the difference in periods is of the order 1.0 ms. For the grid with a mesh size of 40 m, the difference in Hs is at average about

1 cm, whereas the difference in period measures is roughly 10 ms at the four forementioned buoy locations. h = 20m h = 40m (Hs) 10 4 m (Tm 10) 104 s (Tp) 10 4 s (Hs) 10 4 m (Tm 10) 104 s (Tp) 104 s AZB21 10.5 9.0 2.1 104.1 119.0 35.0 AZB22 23.4 9.9 6.0 10.5 31.0 4.0 AZB31 22.5 24.9 15.9 52.2 76.0 15.0 AZB32 14.4 32.1 8.1 11.1 110.0 24.0

Table 2.2. Mesh differences in Amelander Zeegat buoys in area 2 for h = 20m and 40m.

Although the errors on the grid with a mesh size of 100 m are small (see Section 2.2.1) and a 100 m resolution may be sufficient (for area2, area3 and area4), the choice for grids with mesh size of 20 m is better because of the accurate representation of the islands and bottom topography with rapid changes. This requires a sufficiently fine grid, because in the final computations the source terms, especially surf breaking, will be taken into account.

2.3 Accuracy due to spectral grid resolution

The spectral grid discretization also influences the accuracy of the model results. Contrary to the geographic grid discretization, the Richardson correction and doubling the mesh sizes can not be applied to estimate the global errors in wave parameters due to spectral grid resolution ( , ). The reason for this is twofold:

1. The choice of and is more restricted by the physical conditions and formulations. In most cases where the wind growth is important and because of the use of the DIA technique for nonlinear 4 wave-wave interactions, =10o is the default choice (except for swell waves, see Section 4.3). Moreover, due to the use of DIA, = 0.1 appears to be the most optimal resolution.

2. The Richardson correction is only appropriate for uniform grids. The frequency bins, however, are logarithmically distributed.

(18)

x

N N N S

c c c

x (2.1)

where cx and cy are propagation velocities in geographic x and y direction, respectively, c

and c are propagation velocities in spectral and direction, respectively. We neglect the source term S since, it is not relevant in this explanation. Next, we discretize the (homogeneous) equation as follows:

0 N c N c N x cx x

with a short-hand notation for finite differencing of N in the corresponding –space. Since,

c

_x

/

x

may acts like a time derivative 1/ t (as a part of the material derivative), the equation is rewritten as

0 N t c N t c N x

Finally, we define the following dimensionless parameters:

t

c

CFL

and

CFL

c

t

Since, both formulation and interpretation of these dimensionless parameters is similar to the well-known Courant number of a simple wave-like equation, i.e. the ratio of a time step to a cell residence time, they may be called 1D Courant-like numbers for Eq. (2.1). The Courant number helps determine the size of x. A smaller value of the Courant number

causes a smaller value in x, and therefore a more accurate solution. Hence, a small (<1)

Courant number will improve accuracy.

Similarly, the 2D Courant-like numbers for both frequency and direction grids, respectively, are thus given:

y x c CFL c c x y , y x c CFL c c x y

As a rule, the lower the Courant number the more accurate is the wave propagation in geographical space. Hence, the effect of ( , ) to the global accuracy of the energy transport is related to the mesh size ( x, y). In other words, given the spectral resolution

(19)

to reduce the Courant number is the mesh size ( x, y). Note that we are not varying the

spectral space resolution because of the physical modelling constraints as mentioned before. Simulations of the storm situation of January 2, 2005 10:00hr have been carried out. Only propagation terms are included whereas all physical processes are disabled. Both the nested grids of size 100m (area1) and of size 20m (area2, area3 and area4) are used. With respect to the computation using grids of 20m we apply both constant water level and no current as well as space-varying water level and current (maximum flood) obtained from WAQUA

calculations1; see WL (2006b). Due to the presence of a current c 0, so that the Courant number CFL can be computed as well.

Table 2.3 shows the computed Courant numbers in wave buoys in Amelander Zeegat2. Note that the buoys AZB11 and AZB12 are outside the nested grids area2, area3 and area4. In case of no current (and constant water level) the Courant number CFL is relatively low and is of order 0.1. Hence, we may expect that refraction is computed accurate enough in both 100m and 20m grids in the case without current.

Buoys x = 100m (constant water level, no current) x = 20m (constant water level, no current) x = 20m (space-varying

water level and current)

CFL CFL CFL CFL AZB11 0.00877 AZB12 0.00845 AZB21 0.173 0.117 0.114 0.140 AZB22 0.155 0.066 0.0929 0.106 AZB31 0.465 0.137 1.63 1.16 AZB32 0.277 0.141 0.439 0.770 AZB41 0.172 0.089 1.32 0.199 AZB42 0.080 0.092 1.39 0.194 AZB51 0.105 0.066 0.073 0.084 AZB52 0.126 0.080 0.089 0.096

Table 2.3. Computed (maximal) Courant numbers in Amelander Zeegat buoys.

On the other hand, if currents are taken into account, the results are completely different, in particularly, for buoys AZB31, AZB32, AZB41 and AZB42. In these points, the Courant numbers CFL and CFL are of order 1.0 and 1.5, respectively. The origin of this phenomenon lies in the fact that in between the islands Terschelling and Ameland, where the aforementioned buoys are located, the spatial gradients of the current are relatively high compared to the rest of the domain, see WL (2006, Figure 3.7). As a consequence, the magnitude of both c and c is enlarged locally.

Although these Courant numbers are larger than 0.7, it is expected that the accuracy of the wave propagation in geographical space is good enough, at least in a global sense. Note that the criterion for the Courant number is a sufficient requirement, but is not necessary for the

1

For grid with size of 100m of area1, noWAQUA results are available, since currents are expected not to have an effect on offshore conditions.

2

(20)

accuracy. Hence, the computed wave conditions at AZB31, AZB32, AZB41 and AZB42 are not necessarily inaccurate.

2.4 Summary

For the hindcasts of the Amelander Zeegat as outlined in WL (2006b), a number of rectilinear (nested) areas with different resolutions are employed. Two grid sizes are chosen: 100 m and 20 m. As a result, the global errors in the significant wave height Hs, the mean

period Tm-1,0 and (smoothed) peak period Tps due to geographic grid discretization have been

quantified. For the largest area with mesh size of 100 m, the error in Hs is at average about

10 mm, whereas the error in wave period measures is roughly 10 ms and for the nearshore areas with resolution of 20 m, they are 2 mm and 1.0 ms, respectively. Hence, they are considered to be small bearing in mind that the errors due to simplification of the physical processes by means of parameterization used in SWAN are relatively larger. In hindcast studies Royal Haskoning (2003) and WL/Alkyon (2003) showed that these model errors are of the order 5-10 cm for the significant wave height and 1-2 s for the wave period. With respect to the hindcasts of the Amelander Zeegat, however, grids with a resolution of 20 m are preferred because they enable the representation of the large gradients in the bottom at relative short distances.

It must be stressed that this study is based on stationary simulations and thus provides steady-state error estimates. The above conclusions might not be valid for other type of applications like non-stationary simulations.

The accuracy of the wave propagation in geographical space due to the spectral space with given resolution has been quantified by means of the Courant-like numbers. It must be stressed that the spectral resolution can not be varied because of the physical modelling constraints inSWAN. We distinguish between computations with a constant water level and no current and computations with space-varying water level and current. In the first case, the Courant-like numbers appear in the order of 0.1 which means that the numerical accuracy of refraction is more than sufficient. However, in the second case, the Courant numbers are of order 1.0 1.5 mostly between the islands Terschelling and Ameland, due to the relatively large gradients in the current. In the rest of the domain, these numbers are of order 0.1. Since, the relatively large Courant numbers occur only very locally, we considered the computed wave propagation accurate enough, at least in a global sense. The employed spectral resolution ( =0.1 and =10o) is appropriate for future wind-sea dominated computations.

(21)

3 Convergence behavior in Waddensea

computations

3.1 Introduction

Because the action balance equation is solved implicitly, nonlinear terms within that balance have been linearized, and non-linear four-wave interactions and refraction may propagate energy from one quadrant to another, the (stationary) SWAN solution is determined iteratively (Zijlema and Van der Westhuysen, 2005). In WL (2006b) the iteration process runs from n=1 to n=30 and is terminated if the maximum number of iterations is reached or the following criteria for the significant wave height Hs and mean relative wave period Tm01

are both satisfied in at least 99% of all wet grid points (i,j), seeSWAN User Manual (2006):

m 01 . 0 ) , ( ) , ( or 01 . 0 ) , ( ) , ( ) , ( 1 1 1 j i H j i H j i H j i H j i H n s n s n s n s n s (3.1) s j i T j i T j i T j i T j i T n m n m n m n m n m 1 . 0 ) , ( ) , ( or 01 . 0 ) , ( ) , ( ) , ( 1 01 01 1 01 1 01 01 (3.2)

In general, the iterative method can be stopped if the approximate solution is accurate enough. A good termination criterion is very important, because if the criterion is too weak the solution obtained may be useless, whereas if the criterion is too severe the iteration process may never stop or may cost too much computational time. In order to assess the present criteria (3.1) and (3.2) we carry out the simulations with exactly 30 iterations. The reason for choosing 30 iterations is that no changes in significant wave height Hm0 and

almost no changes in mean wave period Tm01 have been found after 30 iterations. In these

simulations, all the physical processes, except triads, are included. With respect to wind growth and whitecapping the adapted Yan (1987) formulation and Alves and Banner (2003) expressions were employed, respectively3, as described in Van der Westhuysen et al. (2006).

3.2 Results with original convergence criterion

Figure 3.1 depicts the convergence behavior of Hs and Tm01, by showing the ratio of the

solution after n iterations and the ‘exact’ solution (obtained after niter=30 iterations), at several buoys in Amelander Zeegat obtained with domain area1 of 100 m mesh size. Apparently, at least 25 iterations are needed to get converged solutions. However, using the criteria (3.1) and (3.2), it appeared that the iterative process terminated after only 8 iterations.

3

(22)

Looking at the convergence behavior obtained with the nested area2 with mesh size of 20m, Figures 3.2 and 3.3 reveal that at least 20 iterations are needed to get stationary wave parameters for both constant water level/no current and space-varying water level/current. Using criteria (3.1) and (3.2), no more that 16 and 14 iterations, respectively, are required. Hence, we may conclude that the present criteria (3.1) and (3.2) are not strict enough to obtain accurate results after termination of the iterative procedure. It was found that the iteration process can converge so slowly that at a certain iteration level n-1 the difference between the successive iterates, H_sn H_sn1 can be small enough to meet the convergence criteria, causing the iteration process to stop, even though the converged solution has not yet been found. In particular, this happens when convergence is non-monotonic such that the process is terminated at local maxima or minima that may not coincide with the converged solution. Furthermore, it became apparent that, unlike Hs, the quantity Tm01 is not an

effective measure of convergence. It was found that the relative error in Tm01, i.e.

1 01 01 01

n n n

m m m

T T T , does not monotonically decrease near convergence, but keeps oscillating during the iteration process (see Figures 3.1, 3.2 and 3.3). This behavior is due to small variations in the spectrum at high frequencies, to which Tm01 is sensitive. This behavior is

problematic when any form of stricter stopping criterion is developed based on Tm01.

Therefore, in the improved termination criterion proposed in Zijlema and Van der Westhuysen (2005), Tm01 has been abandoned as a convergence measure and only Hs, which

displays more monotonic behavior near convergence, is retained.

3.3 Results with the curvature-based convergence criterion

The new stopping criterion of Zijlema and Van der Westhuysen (2005) not only considers the relative error (3.1) but also the second derivative or curvature of the curve traced by the series of iterates (iteration curve). Since the curvature of the iteration curve must tend towards zero as convergence is reached, terminating the iteration process when a certain minimum curvature has been reached would be a more robust break-off procedure. The resulting curvature-based termination criterion is given by

1 2 3 ( ) , 3, 4, 2 n n n n s s s s C n s H H H H n H (3.3)

where C is a given maximum allowable curvature. The curvature measure is made

non-dimensional through normalization with Hs n

. Condition (3.3) must be satisfied in at least 99% of all wet grid points before the iterative process stops. This curvature requirement is the primary criterion. The weaker criterion (3.1) is retained in addition to the stricter criterion (3.3).

The new stopping criteria have been applied with the curvature parameter of C = 0.001.

(23)

Present stopping criteria (3.1) and (3.2)

New stopping criteria (3.1) and (3.3)

x = 100m 8 30

x = 20m (constant water

level, no current) 16 19

x = 20m (space-varying

water level and current) 14 18

Table 3.1. Number of iterations obtained with different stopping criteria.

3.4 Summary

(24)

4 Requirements to the spectral resolution

4.1 Introduction

A key problem in the derivation of the Hydraulic Boundary conditions for the dikes along the Waddensea is the estimation of the amount of penetration of long wave energy through the tidal inlets. Such long waves (swell) are often characterized by a small amount of directional spreading. According to the SWAN manual (SWAN team, 2006), the directional resolution should reflect the directional spreading of the incident wave conditions. For many applications, where swell waves coexist with locally generated wind seas, choosing a swell related fine directional resolution is not practical because it increases the computational requirements considerably. In practice, however, a directional resolution of 10° is chosen. It is not yet known whether this resolution is too crude for computing the penetration of long period waves through the tidal inlet.

The spectral resolution concerns the discrete frequencies and the directions that have to be specified inSWAN. The following two requirements are investigated. The first requirement concerns the frequency range for wind wave growth. The second requirement concerns the required directional resolution for swell penetration and propagation into the Waddensea. This aspect is of special interest for the main land dikes of the Waddensea, since it is assumed that relatively long waves may contribute significantly to the wave loads on these dikes.

In this section results are presented that were computed withSWAN on curvi-linear grids for the area around the Amelander Zeegat. These grids are based on the curvi-linear grid used by WL (2006b). A section of this grid was taken and is referred to as grid AZG1. The outlines of the large and the selected curvi-linear grids are shown in Figure 4.1. The bottom topography for the large curvi-linear grid is shown in Figure 4.2.

Based on grid AZG1 refined grids were made in which the grid resolution was increased with a factor of 2 or 3, yielding the grids AZG2 and AZG3. As outlined in Section 6 of this report modifications to these grids were made in which the grid resolution was further refined within the Amelander Zeegat. These grids are referred to as the grids AZG1A, AZG2A and AZG3A. Further information on these grids is given in Section 6.

In this section two storm events were used for sensitivity studies. The two situations are taken from WL (2006b), their Table 3-5. The first situation is for 8 Febr. 2004, the second one is for 8 Jan. 2005. The wave boundary conditions were imposed on the boundaries of the curvi-linear grid in the same way as in WL (2006b).

(25)

Date Time (MET) Hm0 (m) Tm01 (s) Wdir (°) Dspr (°) U10 (m/s) Udir (°) Water level (m) NAP 20040208 22:20 6.3 7.8 315 29 15 326 2.54 20050108 18:00 6.7 7.5 274 31 19 270 2.29

Table 4.1: Summary of boundary conditions for two storm events.

4.2 Frequency and directional resolution for wind seas

4.2.1 Frequency resolution

In third generation spectral wave models the discrete frequencies are distributed geometrically such that succeeding frequencies differ by a constant factor: fi+1=(1+x)fi,

where x has a typical value of 0.1. This type of spacing is related to the efficient use of the Discrete Interaction Approximation (Hasselmann et al., 1985) for the computation of the non-linear four-wave interactions. The value of x=0.1 implies that each subsequent frequency is 10% higher than the previous. The value of 10% is recommended by Hasselmann et al. (1985) and it is based on tuning experiments of fetch limited wave growth. Van Vledder et al. (2000) showed that choosing a too fine resolution may lead to non-physical wave evolution and they suggest that the DIA only works for frequency resolutions of about 10%. However, the precise limits of x are not yet known. The determination of these limits falls outside the scope of this study.

The frequencies are distributed between a lower limit fmin and an upper limit fmax. The range

of frequencies depends on the range of peak frequencies that can be expected in the whole computational domain. Since in SWAN only one spectral resolution can be specified per computational grid, it may happen that a wide range of peak frequencies may have to be considered. This holds especially for large basins like the Waddensea where both long period (low frequency) swell waves and short (high frequency) wind generated waves may co-exist.

The frequency range should be large enough to include all relevant frequencies within a wave spectrum. This range can be different for each situation. It depends on the characteristics of the incoming North Sea waves, the wind speed and depth. The latter is of importance due to the generation of higher harmonics in shallow water by non-linear three-wave interactions. The frequency range can be related to the peak frequency of the three-wave spectrum (or to the peak frequencies of a multi-peaked spectrum).

The minimum frequency should be about 0.4 times the lowest expected peak frequency. For Waddensea applications the lowest expected peak frequency usually occurs at the incoming boundary located in the North Sea. The value of 0.4 follows directly by inspection of a Pierson-Moskowitz spectrum and looking at the decay of spectral densities below the peak frequency.

(26)

The second maximum should be about 3 times the main peak frequency in shallow areas where triad interactions play a role. In the case the non-linear four-wave interactions are computed using an exact method, a factor 6 is recommended (Van Vledder, 2006). However, application of this exact method is not part of this study.

Application of the above criterion to the peak frequencies that occur just behind the Wadden islands would result in a relatively high maximum frequency. Bottema (2004) estimates this peak frequency on the basis of parametric growth curves using the wind speed U10 and the

size of the first downwind grid cell x. The corresponding peak period is estimated as:

0.27 2 10 10 2 13.7 / p g x U T g U

For the Waddensea application it is more practical to estimate this peak frequency on the basis of parametric growth curves using the wind speed and a fetch of 1 km as input. The occurring peak frequencies can also be estimated in aSWAN trial run.

The effect of choosing different upper limits for the frequency range in SWAN was investigated by performing three computations for the Amelander Zeegat on the curvi-linear grid AZG3. The situation is for 8 Febr. 2004 and is taken from WL (2006b), their Table3-5, see also Table 4.1. This situation refers to an incident wave direction from the Northwest and a wind speed of 15 m/s and a water level of NAP +2.54 m. The computations were performed for three frequency ranges: [0.03-0.85], [0.03-1.25] and [0.03-1.65].

The results of the computation on the smallest frequency domain [0.03, 0.85] were visualized in the form of contour plots of the significant wave height Hm0 and the spectral

period Tm-1,0, Tm01 and Tm02 in the Figures 4.3 and 4.4. The arrows indicate the mean wave

direction, they are scaled with the significant wave height Hm0. Each of these period

measures emphasizes a certain frequency range within the spectrum. These figures indicate that in the computational area the range of spectral periods varies from 1 s just behind the Wadden islands to 10 s on the North Sea. The spatial distributions of the various period measures show interesting features south of the tidal inlet. The spectral periods Tm-1,0 and

Tm01 have a narrow strip with higher periods that stretch southwards into the Waddensea,

whereas this band does not exist for the spectral period Tm02. These differences show that the

influence of North Sea waves south of the tidal inlet is only noticeable for the lower frequencies in the wave spectrum.

The effect of choosing a different upper frequency fmax on the significant wave height and

(27)

The results indicate that the significant wave height Hm0 is hardly affected by choosing a

different frequency band, differences are less than 1 cm. The spectral period Tm-1,0 is least

affected, as expected, whereas the spectral period Tm02 is most sensitive to the choice of the

upper frequency fmax. The differences in parameter values become larger (compared to

choosing fmax= 0.85 Hz) as the upper frequency increases.

The results indicate that the largest differences occur for the spectral periods Tm01 and Tm02

just behind the Wadden islands. The extent of this area is limited to a few kilometers at the lee side of the island, on the ebb-tidal delta, shallow area in between the islands and to some areas close to the dikes where the effect of triads is noticeable at higher frequencies. The wave conditions in the deep channel do not seem to be affected by choosing a different frequency range. For the spectral period Tm-1,0 some relatively large differences occur in the

tidal inlet when an upper frequency fmax of 1.65 Hz is used (See Figure 4.7).

Inspection of the ray plots (Figures 4.9 and 4.10) indicates that differences only occur in the first 5 km of this ray, i.e. over the ebb-tidal delta. Further along this ray, following the main channel, the differences are very small. At most, an over-prediction of 0.25 s occurs near the dike of the main land.

These results indicate that in practice an upper frequency of 0.85 is sufficient to compute the wave conditions in the Waddensea. Since the frequencies are geometrically distributed, addition of a few extra discrete frequencies only results in a small increase in computational requirements. Therefore, for practical applications the upper frequency fmax can be chosen as

1.0 Hz, which is in agreement with the setting recommended in the SWAN user manual (SWAN team, 2006).

4.2.2 Directional resolution

The discrete directions are usually distributed over the full circle with a constant spacing of 10°. For specific applications it may be possible to limit these frequencies to a certain directional sector. Reducing the computational workload is often an argument to choose a limited directional sector. In such a case, one needs to ensure that all relevant directions of wave propagation are included in the schematization.

4.3 The effect of directional resolution on swell

propagation

4.3.1 Model setup

The sensitivity of swell penetration into the Waddensea on the directional resolution of the

SWAN model was studied on the curvi-linear grid AZG3A for the Amelander Zeegat. The characteristics of this grid are discussed in Section 6.

(28)

specified on the full circle. The water level was set at NAP +3m. As recommended in Section 3 the curvature-based convergence criterion was applied with the following settings: NUM STOPC 0.00 0.01 0.001 99.5 STAT mxitst=75

To simulate the penetration realistically, the source term for surf breaking was activated to allow breaking in shallow areas. Test computations indicated that the source term for quadruplet interactions could be turned off because they did not affect the results. This is not surprising since the magnitude of these interactions scales with f11. Therefore, quadruplets have a negligible effect on this long period swell waves compared to wind driven waves with typical periods of 5 s to 10 s (0.2 Hz to 0.1 Hz). An additional benefit is that the computations took less time.

For the sensitivity analysis the following combinations of incident directional spreading and directional resolution were used:

= 5°, 10°;

=2°, 3°, 4°, 5°, 6°, 8°, 10°.

For practical reasons only results for an incident directional spreading of 10° are shown.

4.3.2 Results

The penetration of swell waves is visualized in the form of contour plots of the significant wave height Hm0 and the spectral period Tm-1,0 for a number of combinations of directional

spreading and the directional resolution. Inclusion of all plots is impractical and not necessary for the purpose of this study.

The Figures 4.11 through 4.14 show the spatial variation of the significant wave height Hm0

and spectral period Tm-1,0 for the situation with a water level of NAP +3 m, an incident

directional spreading of 10°, and directional resolutions of 10°, 6°, 4° and 2°, respectively. The arrows indicate the mean wave directions. They are scaled with the square root of the significant wave height to create a suitable distribution of arrow lengths. White spots reflect dry areas. The figures show various interesting features of swell penetration into the tidal inlet. The swell waves shoal considerably on the shallow areas in the tidal inlet. The wave height is much lower in the Borndiep channel (i.e. the main channel in the tidal inlet, running just west of the island of Ameland) than on the shallow areas. This is probably due to the fact that the wave refract out of this channel. This is supported by the model finding that the wave periods are smaller in the Borndiep than on the shallow areas west of this channel. The arrows indicate that the swell waves bend more than 90° such that they propagate towards the southern shore of Terschelling. However, these strongly refracted waves have small heights.

(29)

dividing these differences by the results for the coarser resolution. The Figures 4.15, 4.16 and 4.17 show the differences for the significant wave height Hm0 and spectral period Tm-1,0.

The difference plots show practically no variation in parameter values on the North Sea. In the tidal inlet and in the Waddensea large variations exist between the results for different resolutions. Increasing the directional resolution, generally leads to an increase of the significant wave height in a narrow strip south of the tidal inlet of at most 15%. On both sides of this strip the wave height decreases with about 20%. The effect on the spectral period due to an increased resolution is much smaller. The wave period remains more or less constant in the narrow band south of the tidal inlet, and the differences are generally smaller than 1%.

The largest differences occur between the computations with directional resolutions of 2° and 10°. The differences are generally smaller between the results for the cases 6° and 2°, and even smaller between the results for the cases 4° and 2°. This behavior suggests that at some point a further increase in directional resolution does not lead to different results. This assumption was tested by plotting the convergence behavior of a number of integral wave parameter as a function of directional resolution in a number of output points.

For six output points south of the tidal inlet the convergence behavior of 5 integral wave parameters is plotted as a function of the directional resolution in the Figures 4.18 through 4.23. These results indicate that the wave parameters show a persistent trend with directional resolution. With increasing directional resolution, i.e. decreasing directional step , the significant wave height Hm0 increases in most output points, whereas the period

measures and the directional spreading decrease. It is also evident that the results are not yet converged with respect to the directional resolution . For the points under consideration, the significant wave height increases by 10% when the directional resolution is increased from 10° to 2°. The periods show a much weaker decrease of at most 1%. The directional spreading also shows a decrease of about 10%.

Of interest is also the spatial variation of the significant wave height Hm0, the spectral period

Tm-1,0, the mean direction and the directional spreading as a function of the directional

resolution . To that end, 11 cross sections are defined in grid AZG3A along which the variation of the above wave parameters as a function of the directional resolution are presented. The locations of the cross sections are shown in Figure 4.24. The most illustrative variations of the wave parameters are for the cross sections 6, 9 and 11. These variations are shown in the Figures 4.25 through 4.27.

The results for cross section 6 (Figure 4.25) show that the highest wave heights are found on the shallow area just west of the deepest channel (i.e. the Borndiep). This is possibly due to the effect of refraction which turns the waves to the sides of the channel. The results also indicate that the convergence behavior (w.r.t. the directional resolution) of the wave parameters is monotonous for large stretches. The parameters values either increase or decrease. The relative variation with directional resolution is larger for the significant wave height Hm0 than for the spectral period Tm-1,0. A strong variation with directional resolution is

(30)

Along cross section 9 (Figure 4.26) it is found that for the stretch around x=172 km the significant wave height Hm0 decreases with increasing directional resolution, whereas the

spectral period Tm-1,0 increases. For x=170 km it is found that the significant wave height

Hm0 increases with increasing directional resolution, whereas the spectral period Tm-1,0

hardly changes. For x=177 km the mean wave direction shows variations up to 40° and the directional spreading show a typical difference of 5°. However, for x=166 km and x=174 km the differences become as large as 30°.

Along cross section 11 (Figure 4.27) the variation in the significant wave height Hm0 with

directional resolution is generally smaller than for the rays 6 and 9. The variation of the spectral period Tm-1,0, however, shows a strong variation with directional resolution for x

larger than 174 km. The variation in the mean directional is small, whereas the directional spreading shows variations up to 10°.

The different results for different directional resolutions cannot be attributed to the poorer iteration behaviour of SWAN for finer directional resolutions. Only for a directional resolution of 2° convergence was not reached in some output point. See Section 4.3.3 for a discussion of the convergence behaviour ofSWAN.

The present computations were carried out without the effect of diffraction. It is possible that diffraction effects smooth the results in such a way that convergence with respect to the directional resolution is obtained. It is also possible that diffraction effects occur around the many shallow areas in the Waddensea. It is therefore recommended to assess the effect of diffraction using dedicated refraction-diffraction models.

4.3.3 Convergence behavior

Finally, the convergence behavior of theSWAN model was determined for all computations. The results for the cases with a directional resolution of 10°, 5°, 3° and 2° are shown in the Figures 4.28 through 4.31. The upper left panel shows the percentage of accepted points according to the curvature criterion. The upper right panel shows the locations of the test points. Each other panel shows the convergence behavior of one parameter for all test points. The location and results for one output point are highlighted in red.

The results of the convergence behavior indicate that with increasing resolution, more iterations are needed to reach acceptable convergence. For the case with a directional resolution of 2° convergence has not been reached in all output points after 33 iterations, despite the fact that the required number of accepted points was reached. Test computations (not shown here) with a 1° directional resolution showed that no convergence could be reached after 100 iterations in most points.

(31)

4.4 Summary

The results of this section indicate that the required spectral resolution depends on the conditions that need to be computed. The frequency range depends on the characteristics of the incoming waves and on the wind speed. The incident waves mainly affect the lowest frequency, whereas the wind speed determines the upper frequency.

The test carried out in this study indicate that choosing a higher upper frequency mainly affects the waves at the lee side of the islands and in the surf zone on the northern side of the Wadden islands. Further into the Waddensea the effect of choosing a higher upper frequency diminishes. These results are not surprising. At the lee side the waves are very small and it makes sense to choose an upper frequency higher than, say, 1 Hz. Note, that the value of 1 Hz is recommended in the SWAN manual (SWAN team, 2006). Choosing an even higher frequency does not increase the computational requirements considerably since the frequencies are distributed geometrically. The sensitivity of the results on the northern side of the island is probably related to the effects of the triads, which generate higher harmonics. A too small upper frequency will limit the generation of some higher harmonics, thereby affecting the mean periods.

The SWAN manual (SWAN team, 2006) recommends different directional resolutions for incident waves with different directional spreading. For most wind sea situations a directional resolution of 10° is sufficient. For narrow swell waves with a directional spreading of 10° or less a finer directional resolution is required.

The results presented in this section indicate that choosing different directional resolutions leads to different amounts of swell penetration into the Waddensea. From a numerical point of view, choosing a smaller directional resolution will theoretically lead to more accurate results of wave penetration. However, from an operational point of view, choosing a smaller directional resolution for mixed sea conditions leads probably to unfeasible computational requirements, both in memory requirements and in computing time per simulation.

The results also indicate that the computational results have not converged with respect to the directional resolution. Evidently, the results depend on the directional resolution. The results also indicate that the significant wave height is rather sensitive to the directional resolution, whereas the sensitivity of the spectral periods to the directional resolution is much smaller. Further, it is found that more iterations are needed with increasing directional resolution. Unfortunately, based on the present computational results no firm conclusion can be drawn about the recommended directional resolution for swell penetration into the Waddensea.

The present behavior (i.e. non-convergence w.r.t. directional resolution) might due to the complexity of the Amelander Zeegat, where the large variations in bottom level affect the refraction patterns and the resulting wave conditions. It is therefore recommended to analyze the sensitivity of the swell penetration on the directional resolution for an idealized and highly schematized tidal inlet.

(32)

peaked frequency spectra when a very fine frequency resolution of 3% was used. It is not yet clear if a similar effect occurs for low frequency swell waves when a very fine directional resolution is chosen. The present computations were carried out with long period swell waves, for which by definition, the quadruplet interactions are much smaller than for wind sea waves. The magnitude of these interactions scales with fc11, with fc a characteristic

frequency. It is therefore suggested to introduce a threshold frequency for activating the DIA to reduce the computational requirements. Further sensitivity analyses need to be carried out to investigate the effect of choosing a small directional resolution on wave evolution using a discrete spectral wave model that uses the DIA.

Sensitivity analysis of numerical aspects of SWAN: Activity 3.1 of SBW project Waddensea