Energy transport along wave rays

In this section we focus on the following transport processes: shoaling and refraction. Calculating the wave shoaling and refraction effects is necessary to predict accurately shallow water wave conditions, either in the surf zone, across the main channels in estuaries, or across the seamounts. Throughout this section we assume the absence of the non-conservative source/sink terms, such as wind input, nonlinear wave-wave interactions and energy dissipation. The governing equation is the following wave energy balance (ambient current is not included)

  $$\frac{\partial E}{\partial t} + \nabla_{\vec{x}} \cdot ({\vec{c}}_g E) + \frac{\partial c_\theta E}{\partial \theta} = 0$$ (3.46)

Eq. (3.46) is linear and strict hyperbolic with nonlinear coefficients, whereas its characteristics (trajectories) in the $(\vec{x},\theta)-$space are generally not straight lines due to a varying seabed topography. Along the characteristics the wave energy fluxes ${\vec{c}}_g E$ and $c_\theta E$ are constant (Whitham, 1974).


We may rewrite Eq. (3.46) in a characteristic form, as follows

  $$\frac{dE}{dt} = -\left( \nabla_{\vec{x}} \cdot {\vec{c}}_g + \frac{\partial c_\theta}{\partial \theta} \right) E$$ (3.47)

with the total derivative of $E$ defined as

  $$\frac{dE}{dt} \equiv \frac{\partial E}{\partial t} + {\vec{c}}_g \cdot \nabla_{\vec{x}} \, E + c_\theta \frac{\partial E}{\partial \theta}$$ (3.48)

along a trajectory of energy propagation with slopes

  $$\frac{d \vec{x}}{dt} = {\vec{c}}_g = (c_{g,x},c_{g,y})\, , \quad \frac{d\theta}{dt} = c_\theta$$ (3.49)

Let us consider a stationary wave characteristic with slope

  $$\frac{d \vec{x}}{d\theta} = \frac{{\vec{c}}_g}{c_\theta}$$ (3.50)

If there is no change in the seabed along the wave ray, the group velocity is constant. In addition, if there is no depth variation along the wave crest, then the turning rate is zero. Hence, the characteristic (3.50) in the $(\vec{x},\theta)-$space is in this case the direction of propagation (with group velocity). By virtue of Eq. (3.47), the total energy $E$ is thus constant along the same characteristic. However, due to depth variations, wave energy will either increase or decrease along its curved characteristic.


The right hand side of Eq. (3.47) can be considered either as a source term or a sink term. This depends on the gradients of the group velocity and turning rate along the wave ray and wave crest, respectively. The rate at which the energy in- or decreases is related to the relaxation time $\tau$, which is the typical time scale for wave energy transport to reach steady state after being disturbed. It is given by

  $$\tau^{-1} = \left \vert \nabla_{\vec{x}} \cdot {\vec{c}}_g + \frac{\partial c_\theta}{\partial \theta} \right \vert$$ (3.51)

The importance of this relaxation time relates to the choice of discretization steps for the purpose of accurate integration of Eq. (3.47). In this respect, the numerical accuracy with respect to the change in wave energy from one step to the next along the characteristic, viz. Eq (3.48), is determined by the time step, grid size, and directional bin size (see Section 3.8.3). The associated step size is denoted by $\Delta T$. Causality requires that wave energy propagates in the right direction along its characteristic, and at the right speed. For example, causality problem can be present in an implicit scheme that propagates wave energy across a large distance using a large time step. In other words, if this time step is too large, some wave components getting ahead of themselves and leaving behind some other components ahead. To prevent this, temporal, spatial and directional changes in the numerical and exact solutions must go hand in hand. This implies that $\Delta T < \tau$. Hence, a sufficient condition for accurate integration reads

  $$\left \vert \nabla_{\vec{x}} \cdot {\vec{c}}_g + \frac{\partial c_\theta}{\partial \theta} \right \vert \, \Delta T < 1$$ (3.52)

This condition resembles the Lipschitz criterion^3.6and is generally less severe than the traditional CFL criterion for numerical stability. We may rewrite Eq. (3.52) in Cartesian coordinates as

  $$\left \vert \frac{\partial c_{g,x}}{\partial x} + \frac{\partial ... ...+ \frac{\partial k}{\partial y}\sin\theta \right) \right \vert \, \Delta T < 1$$ (3.53)

with $k = \vert\vec{k}\vert$ the length of the wave number vector (cf. Eq. (2.14)). This criterion implies that at locations with relatively large bottom slopes, the step size $\Delta T$ must reduce locally to prevent inaccuracies. This step size is determined by the temporal, spatial and directional resolutions. In this regard, the interpretation of the Lipschitz criterion (3.52) is that the maximum step size is determined by the numerical accuracy rather than by the numerical stability. This accuracy aspect is related to the curvature of the wave propagation field (due to change in wave direction over a certain distance) at the grid and time resolutions applied.


In shallower water, medium variations can be significant, i.e. changes in seabed and mean current can be large within a few wave lengths, whereas ocean waves feature coherent structures such as refraction over topography and currents. Therefore, one must choose a geographical grid size proportional to the resolution of the bathymetry capturing its local features. For instance, near the coast and in the surf zone, a typical grid size of 20 $-$ 50 m is not uncommon in SWAN applications. An important assumption made in this consideration is that the transport velocities $\vec{c}_g$ and $c_\theta$ do not much vary over a grid size, so that Lipschitz criterion (3.52) is most likely met. This is reasonable as the wave characteristics are more or less non-curved lines (during the elapsed time step), because the propagation in the $(\vec{x},\theta)-$space is slowly time varying. This is generally true in coastal applications.


In the open ocean, it is assumed that seabed topography is slowly varying in space so that the directionally spread wave field is rather spatially homogeneous. As such, refraction effects can be regarded as weak. The grid resolution is usually determined by the resolution for the wind field that generates the waves locally, which is often much coarser than the bottom resolution. For example, for oceanic waters, the employed grid size is typically 10 $-$ 50 km. However, there are exceptions. An example are the seamounts in the deeper parts of the open ocean. At such locations, the depth may vary from one grid point to the next by a factor of 10 or so or even more. In such a case, the value of the turning rate $c_\theta$ at the shallowest grid point is not representative anymore with respect to the distance between the two considered grid points. Hence, criterion (3.52) is violated, and computed spectral wave components will simply turn too much and jump multiple directional bins in one distance step. The result is that a number of wave components in one sweep within one distance step overtakes some other bins in another sweep ahead, which implies that causality is violated, resulting in a possible model instability. To circumvent this we must refine the computational grid locally to resolve the local bathymetric features, and the Lipschitz criterion (3.52) can be helpful in this. The consequences will be discussed in the following section.