Introduction

The accuracy with which physical processes for wave growth are approximated numerically is of crucial importance in assessing the predictive realism of spectral wave models. There is a need to separate these numerical errors from errors due to physical modelling. Third-generation wave models pose a numerical difficulty caused by the presence of multiple time scales. This is a reflection of the physical nature of wind waves, which consist of a wide range of frequencies. The ratio of the largest to the smallest time scale of spectral components is often substantially larger than one. When this is the case, the action balance equation is called stiff (Press et al., 1993)^3.4. Taking proper account of these time scales is a necessary condition for numerical accuracy. This would require the use of a very small time step in a numerical algorithm, which may be impractical. Moreover, the action balance equation is usually so stiff that its numerical implementation combined with economically large time steps often prevent a stable solution. In this respect, nonlinear four-wave interaction usually poses the biggest problem, since this process is associated with high sensitivity to spectral change.


In a number of papers concerning spectral wave computation, numerical measures are proposed to achieve stable model results economically. WAMDI Group (1988) suggest to use a semi-implicit time integration scheme with a time step that matches the time scale of low-frequency waves. However, numerically stable solution of the resulting sytem of equations can not be guaranteed (Hargreaves and Annan, 2001). The ratio of the largest eigenvalue to the smallest eigenvalue of the stiff system of equations, called the condition number, can be so large that even a fully-implicit method combined with large time steps precludes a stable solution. For counterexamples, see Hargreaves and Annan (2001). The only remedy is time step reduction or under-relaxation so that the modified system of equations has a spectrum of eigenvalues with a more favourable condition number.


To guarantee numerical stability at relatively large time steps, the so-called action density limiter has been introduced in WAM in the early 1980's (Hersbach and Janssen, 1999). This limiter restricts the rate of change of the energy spectrum at each time step. Because low-frequency waves carry the most energy, it is desirable to solve the balance equation in this part of the spectrum accurately without intervention by the limiter, whereas for high-frequency waves using an equilibrium level is sufficient. Although this approach lacks a rigorous foundation and is not generally applicable or valid, it appears to guarantee numerical stability at relatively large time steps even when these do not match the time scales of wave growth. Moreover, it is believed that the limiter will not affect the stationary solution when convergence is reached. This assumption is widely employed as a justification for the use of limiters. For an overview, we refer to Hersbach and Janssen (1999) and Tolman (2002) and the references quoted therein. Tolman (1992) proposes an alternative to the action density limiter in which the time step is dynamically adjusted where necessary to ensure accurate wave evolution. The calculation of this optimal time step is related to the action density limiter. Further details can be found in Tolman (1992, 2002). .


The steady-state solution in the SWAN model is obtained in an iterative manner, which can be regarded as a time marching method with a pseudo time step. This pseudo time step generally does not match the relatively small time scale in frequency space and consequently, divergence will occur. Therefore, SWAN makes use of the action density limiter to stabilize the iteration process (Booij et al., 1999). However, experience with SWAN has revealed that the limiter acts not only in the equilibrium space, but also in the energy-containing part of the wave spectrum. This finding is also confirmed by Tolman (2002). Furthermore, the limiter appears to be active over almost all spectra in the geographical domain and during the entire iteration process. This activity has been associated with poor convergence behaviour, such as small-amplitude oscillation in frequency space. Ris (1999) demonstrated that stationary SWAN results are influenced by the settings of the action limiter while De Waal (2001) suspects that the limiter acts as a hidden sink in the source term balance under equilibrium conditions. The question to what extent this limiter adversely affects the stationary solution of SWAN has not been addressed previously, and is considered here.


An alternative way to restrict the high rate of change at higher frequencies is under-relaxation, i.e. making smaller updates by means of a much smaller (pseudo) time step (Ferziger and Perić, 1999). Consequently, a limiter may no longer be needed. Although this approach may be suitable to SWAN, it slows down convergence significantly. Here, we propose a new method that finds a compromise between fast convergence on the one hand and minimizing the role of the limiter in the energetic part of the spectrum on the other. The key idea to achieve this is to link the extent of updating to the wave frequency - the larger the frequency, the smaller the update. This approach is therefore called frequency-dependent under-relaxation.