The problem with refraction in non-stationary applications

In this section we discuss how refraction affects the accuracy of the discretization of the total derivative of wave energy along a non-curved wave ray. Initially, the energy is uniformly distributed over the wave directions. Under these conditions, we reconsider Eq. (3.47),

  $$\frac{dE}{dt} = -\frac{\partial c_\theta}{\partial \theta}\, E$$ (3.54)

and we examine the discretization of the following total derivative of the energy density

  $$\frac{dE}{dt} = \frac{\partial E}{\partial t} + c_x \frac{\partial E}{\partial x} + c_y\frac{\partial E}{\partial y}$$ (3.55)

with $(c_x,c_y) = {\vec{c}}_g$ the propagation velocity vector in the geographical space. This total derivative indicates that the time rate of change in energy is computed along the wave characteristics defined by the following ODEs

  $$\frac{dx}{dt} = c_x\, , \quad \frac{dy}{dt} = c_y\, .$$ (3.56)

For the purpose of illustration the spatial derivatives are replaced by first order upwind differences. If we assume $c_x$ and $c_y$ positive during the first sweep, see Figure 3.4,

Figure 3.4: Numerical scheme for wave propagation in geographic space with below the first quadrant for which the waves are propagated, and right the stencil.

then the approximation of Eq. (3.55) is as follows

  $$\frac{E^n_{i,j}-E^{n-1}_{i,j}}{\Delta t} + c_x \frac{E^n_{i,j} - E^n_{i-1,j}}{\Delta x} + c_y \frac{E^n_{i,j} - E^n_{i,j-1}}{\Delta y} \, .$$ (3.57)

Note that the time integration is based on the first order implicit Euler scheme. Furthermore, $\Delta t$ is the time step, and $\Delta x$ and $\Delta y$ are the mesh spaces. This approximation can be viewed as the well-known semi-Lagrangian approximation, which is rewritten as

  $$\left ( \frac{1}{\Delta t} + \frac{c_x}{\Delta x} + \frac{c_y}{\D... ...1}_{i,j} - \frac{c_x}{\Delta x} E^n_{i-1,j} - \frac{c_y}{\Delta y} E^n_{i,j-1}$$ (3.58)

which is, in turn, interpreted as an approximation of the total derivative, Eq. (3.55), as follows

  $$\frac{E^n_{i,j} - E^{n-1}_{i^*,j^*}}{\Delta T}$$ (3.59)

with

  $$i^* = i - p\, , \quad j^* = j - q\, , \quad p = \frac{c_x \Delta T}{\Delta x}\, , \quad q = \frac{c_y \Delta T}{\Delta y}\, .$$ (3.60)

Note that $p$ and $q$ are the Courant numbers. They are not integers, and therefore $(i^*,j^*)$ is not a grid point. This point, however, lies on the wave characteristic. The quantity $E^{n-1}_{i^*,j^*}$ can be interpreted as the value of $E$ at time $t^{n-1}$ in $(i^*,j^*)$ which is being convected in $(i,j)$ in a lapsed time $\Delta T$. This value is simply obtained from the interpolation of the surrounding values $E^{n-1}_{i,j}$, $E^{n}_{i,j}$, $E^{n}_{i-1,j}$ and $E^{n}_{i,j-1}$ in the $(t,x,y)-$plane. In this case, there is no restriction on time step $\Delta t$ as the characteristic lies inside the computational stencil. Note that the lapsed time $\Delta T$ is a function of time step $\Delta t$ and grid sizes $\Delta x$ and $\Delta y$, and is called the Lagrangian time step, which should not be confused with the Eulerian time step $\Delta t$. Generally, $\Delta T < \Delta t$.


Causality requires that wave energy being refracted from one directional bin to another further down must pass through all the bins along a wave crest between them. Hence, for an accurate time integration along the wave characteristic, viz. Eq. (3.59), we must employ the Lipschitz criterion (3.52). In the present context, this criterion reads

  $$\vert \frac{\partial c_\theta}{\partial \theta} \vert \, \Delta T < 1$$ (3.61)

which ensures that wave energy propagate from a bin to the adjacent one during one time step $\Delta T$. Thus, for physical consistency, a restriction on the time step must be imposed in order for the wave directions not to cross each other and the boundaries of a quadrant in the spectral space. Specifically, the Courant number based on $\Delta T$ and $\Delta \theta$ (i.e. directional bin) must be less than unity, that is,

  $$\mbox{Cr} \equiv \frac{\vert c_\theta\vert \Delta T}{\Delta \theta} < 1$$ (3.62)

with $c_\theta$ the turning rate^3.7. This condition is a sufficient one and implies the Lipschitz criterion.


A violation of this criterion implies that the energy can travel in one time or distance step over a number of directional bins or more than the length of one sweep (which in the absence of a current is 90$^o$, see Figure 3.4). For example, consider the first sweep, see Figure 3.4, the boundaries of the first quadrant are the lines 0$^o$ and 90$^o$. Next, we consider the directional sector in the spectral space associated with the considered sweep, see Figure 3.5.

Figure 3.5: Directional sector associated with the first sweep.

We shall show that condition, Eq. (3.62), is sufficient to assure that 1) during $\Delta T$ the distance travelled in $\theta-$direction is at most $\Delta \theta$ and 2) wave energy propagating in any bin in $\theta-$direction will not cross the boundaries of the directional sector, except for the first and last bins. Hence, this prevents wave rays from intersecting each other. Note that we implicitly assume that the net change in mean wave direction within the distance covered during $\Delta T$ is less than the directional resolution. Since, by definition,

  $$\frac{d\theta}{dt} = c_\theta\, ,$$ (3.63)

we may approximate this, as follows,

  $$\theta^n \approx \theta^{n-1} + c_\theta \Delta T\, .$$ (3.64)

For a given $n-1$, we choose an arbitrary point, $\theta^{n-1}$, inside the directional sector; see Figure 3.5. If $\theta^{n-1} > \Delta \theta$ then, Eqs. (3.64) and (3.62) imply $\vert\theta^n-\theta^{n-1}\vert < \Delta \theta$, and hence the chosen point at next time step, $\theta^n$, still lies inside the directional sector. If $\theta^{n-1} < \Delta \theta$, i.e. in the first bin, then one obtains

  $$\theta^n = \theta^{n-1} + \frac{\Delta T}{\Delta \theta} \left[ (... ...1}) c_\theta(\theta=0^o) + \theta^{n-1} c_\theta(\theta=\Delta \theta) \right]$$ (3.65)

or

  $$\theta^n = \theta^{n-1} \left ( 1 - \frac{\Delta T}{\Delta \theta... ...rac{\Delta T}{\Delta \theta} \theta^{n-1} c_\theta(\theta=\Delta \theta) \, .$$ (3.66)

If $c_\theta(\theta=0^o) > 0$ and $c_\theta(\theta=\Delta \theta) > 0$ then, the point at next step will be keep inside the directional sector, if Eq. (3.62) holds. In other cases the energy is leaving through boundary $\theta = 0^o$, though the change in direction is limited to the adjacent directional bin of the other quadrant. This holds also for the last bin and the corresponding right boundary $\theta = 90^o$.


Note that Eq. (3.62) is not required for the stability of the method but contributes to improve its physical accuracy, so that the principle of causality is obeyed. For instance, in a large-scale application, one often applies coarse (nested) grids with a desired time step. Both the mesh width and the time step might be too large to represent wave refraction appropriately. This is readily seen as follows. We consider stationary, long-crested waves in $(x,\theta)-$space, propagating at an angle with the positive $x-$axis. There is a bottom gradient along this axis so that refraction is present. There are no currents. The associated characteristic or wave ray is given by $d\theta/dx = c_\theta/c_g$. This is the spatial turning rate, i.e. the change in wave direction per unit forward distance. Hence, the directional turn of the wave crest over a distance $\Delta x$ is given by

  $$\frac{c_\theta}{c_g}\,\Delta x$$ (3.67)

and the Courant number becomes (cf. Eq. (3.62))

  $$\mbox{Cr} = \frac{\vert c_\theta\vert}{c_g} \frac{\Delta x}{\Delta \theta}$$ (3.68)

which must be smaller than unity to prevent crossing of wave rays. The turning rate of the wave direction per unit time is given by

  $$c_\theta = -\frac{1}{k} \frac{\partial \sigma}{\partial h}\frac{\partial h}{\partial m}\, ,$$ (3.69)

where $k$ is the wave number, $\sigma$ is the frequency, $h$ is the water depth and $m$ is the coordinate along the wave crest (i.e. orthogonal to the propagation direction). On a coarse grid, the depth difference in two adjacent grid points $\Delta h$, and thereby $c_\theta$, can be very large, especially for low-frequency components in very shallow water (see also Section 3.8.5). This implies that the Lipschitz condition may be violated, i.e. $\mbox{Cr} \geq 1$.


Complying with the Lipschitz criterion (3.62) simply prevents wave energy to jump over a number of directional bins, where this energy would go way beyond some other bins ahead. This was demonstrated to be a useful criterion with respect to the behavior of refraction when the bathymetry is poorly resolved by the model. As the refraction becomes excessive in a region with steep bottom gradients, it is possible that the wave rays falsely cross and that the wave energy focus toward a single grid point, creating unrealistically large wave heights and long periods; see Dietrich et al. (2013).


Yet the value of $c_\theta$ in the shallowest grid point can be simply too large due to a large difference in bottom levels over one mesh width, so that wave energy will change direction over more than some directional bins or even the directional sector, so that wave rays falsely intersect each other. To prevent this artefact a limitation on $c_\theta$ seems to be justified. Recalling Eq. (3.62), a limitation would be

  $$\vert c_\theta\vert < \frac{\Delta \theta}{ \Delta T} \, .$$ (3.70)

We estimate the reciprocal of the elapsed time $\Delta T$ as a fraction of

  $$\frac{1}{\Delta t} + \frac{c_x}{\Delta x} + \frac{c_y}{\Delta y}$$ (3.71)

and so,

  $$\vert c_\theta\vert \leq \alpha_\theta \Delta \theta \left ( \fra... ...} + \frac{\vert c_x\vert}{\Delta x} + \frac{\vert c_y\vert}{\Delta y} \right )$$ (3.72)

with $\alpha_\theta$ a user-defined maximum Courant number, which is generally smaller than 1. (In SWAN, the default value is $\alpha_\theta=0.9$.) Often the desired time step $\Delta t$ is such that the first term between brackets can be safely neglected. This implies a slightly more restriction on the turning rate and, in turn, a safety margin in the CFL condition, as follows

  $$\vert c_\theta\vert \leq \alpha_\theta \Delta \theta \left ( \frac{\vert c_x\vert}{\Delta x} + \frac{\vert c_y\vert}{\Delta y} \right ) \, .$$ (3.73)

It must be stressed that this limitation may affect the solution locally depending on $\alpha_\theta$. In fact, we need to find a good estimate for $c_\theta$ which solely depends on $\alpha_\theta$. However, it is unlikely that this measure affects the solution nearshore or on fine grids, since the turning rate will not too much vary over one (spatial) step. This is an effective survival measure in the sense that it prevents the excessive refraction without deteriorating the solution elsewhere.