The problem with refraction on coarse grids

Another issue is the accuracy with which the turning rate is computed on coarse grids. In SWAN this turning rate is computed as follows (see also Eq. (3.69))

  $$c_\theta = \frac{\sigma}{\sinh \,2kh} \left( \frac{\partial h}{\partial x}\sin \theta - \frac{\partial h}{\partial y}\cos \theta \right)$$ (3.80)

Here, wave refraction can only be caused by depth variation. As an alternative, the turning rate can be formulated in terms of phase velocity as follows (see Holthuijsen, 2007, pg. 210)

  $$c_\theta = -\frac{c_g}{c} \frac{\partial c}{\partial m}$$ (3.81)

or

  $$c_\theta = \frac{c_g}{c} \left( \frac{\partial c}{\partial x}\sin \theta - \frac{\partial c}{\partial y}\cos \theta \right)$$ (3.82)

An advantage of this formula is that refraction due to mud (non-rigid seafloor) can be included, which is implemented in SWAN version 41.01. Although these formulas are identical, they differ in result due to numerics. The first one, Eq. (3.80), seems to be rather inaccurate at relative coarse grids with steep bottom slopes. Experiments suggested that Eq. (3.82) with coarse resolution yields results that are similar to those using Eq. (3.82) or Eq. (3.80) with high resolution. By contrast, approximation based on Eq. (3.80) with coarse resolution yields considerably different $-$ often inaccurate $-$ result. Therefore, Eq. (3.80) is replaced by Eq. (3.82) since version 41.01AB.


Till version 41.01A, the derivative $\partial h/\partial x$ or $\partial c/\partial x$ has been approximated using a first order backward difference scheme, for both structured and unstructured grids,

  $$\frac{\partial c}{\partial x} \approx \frac{c_{i,j} - c_{i-1,j}}{\Delta x}$$ (3.83)

see Figure 3.4 for the used stencil. Note that grid point $(i,j)$ is the shallowest one. This approximation appeared to be rather inaccurate at coarse grids as well. Moreover, it can also lead to non-physical asymmetry in turning rate, and therefore wave energy. Therefore, since version 41.01AB, second order central differences are applied as follows

  $$\frac{\partial c}{\partial x} \approx \frac{c_{i+1,j} - c_{i-1,j}}{2\Delta x}$$ (3.84)

Hence, the refraction velocity, Eq. (3.82), is approximated in SWAN, using structured grids, as follows

  $${c_\theta}_{i,j} = \frac{{c_g}_{i,j}}{c_{i,j}} \left( \frac{c_{i+... ...lta x}\sin \theta - \frac{c_{i,j+1} - c_{i,j-1}}{2\Delta y}\cos \theta \right)$$ (3.85)

Note that grid point $(i,j)$ is the shallowest one and that the division by $c_{i,j}$ is not correct, i.e. it is not consistent with the Snel's law! It will overestimate the rate of turning. This error becomes rather large when bottom slopes are exceptionally large so that wave energy may turn over more than one directional bin. This justifies again the use of the refraction limiter, Eq. (3.73). An appropriate upper bound is obtained with $\alpha_\theta=0.9$ that particularly holds for relative long waves^3.8. For shorter waves, a smaller upper bound may be chosen (e.g. $\alpha_\theta = 0.5$). However, one may choose a larger CFL upper bound. For instance, referring to Figure 3.6, i.e. parallel depth contours within 90$^o$, waves can not turn more than 90$^o$ (in line with the Snel's law), which implies $\alpha_\theta = 9$, if $\Delta \theta = 10^o$.


In case of unstructured grids, first order approximations for the gradient of depth or phase velocity have been replaced by a more accurate formula based on the Green-Gauss formula like Eq. (8.36).