Discretization in geographical space

Since the unknown $N$ and the propagation velocities are only given in points $(i,j,l,m)$, further approximation is needed. A first order upwind scheme in geographical space may be employed, since it is fully monotone, i.e. it can not to give rise to spurious oscillations. A disadvantage of this scheme is that it is numerically diffusive, which naturally degrades the accuracy of the model. This numerical diffusion is caused by gradients of wave action across geographic space, e.g. due to refraction by bathymetry or currents, which is often small in coastal areas. So the wave energy field can be considered as smooth. However, in the current SWAN version, two alternatives to this scheme are implemented, namely the second order SORDUP and Stelling/Leendertse schemes. These schemes produce far less numerical diffusion and are appropriate for ocean and shelf sea (regional) applications.


First order upwind scheme; BSBT


The fluxes $c_x N$ at $(i+1/2,j,l,m)$ and $c_y N$ at $(i,j+1/2,l,m)$ are approximated with an upwind scheme as follows

  $$c_x N\vert _{i+1/2,j,l,m} = \left\{ \begin{array}{l} c_x N\vert... ... c_x N\vert _{i+1,j,l,m}\,, \quad c_x\vert _{i,j,l,m} < 0 \end{array} \right.$$ (3.3)

and

  $$c_y N\vert _{i,j+1/2,l,m} = \left\{ \begin{array}{l} c_y N\vert... ... N\vert _{i,j+1,l,m}\,, \quad c_y\vert _{i,j,l,m} < 0 \end{array} \right.\, .$$ (3.4)

The fluxes at $(i-1/2,j,l,m)$ and $(i,j-1/2,l,m)$ are obtained from (3.3) and (3.4), respectively, by decreasing the indices by 1 in appropriate manner. According to the solution algorithm, to be explained later in Section 3.3, both fluxes at $(i+1/2,j,l,m)$ and $(i-1/2,j,l,m)$ acts together with the same sign of $c_x$ (either positive or negative). Note that the flux $c_x\,N$ is considered as a whole quantity, and the propgation velocity $c_x$ is taken in its points of definition^3.1. The same holds for the fluxes at $(i,j+1/2,l,m)$ and $(i,j-1/2,l,m)$ operating together with the same sign of $c_y$. These updates take place by ordering the grid points such that points solved later have no influence on the previous grid points. This ordering and the associated updates are carried out, for instance, during the first sweep if $c_x > 0$ and $c_y > 0$, and the resulting schemes read

  $$\left( \frac{(c_x N)_{i} - (c_x N)_{i-1}}{\Delta x} \right)^{n}_{j, l, m}$$ (3.5)

and

  $$\left( \frac{(c_y N)_{j} - (c_y N)_{j-1}}{\Delta y} \right)^{n}_{i, l, m}$$ (3.6)

In this way the solution at a grid point is determined solely by its upstream points, and causality is preserved. In addition, the above schemes represent the divergence of energy flux, $\nabla_{\vec{x}} \cdot ({\vec{c}}_g N)$, at the discrete level. If there is no net flux change in the corresponding grid point, the divergence must be zero, which is obviously the case. For further details on this solution update, see Section 3.3.


The combination of the time and geographic space discretizations in Eqs. (3.2), (3.3) and (3.4) is known as a first order, backward space, backward time (BSBT) scheme. This scheme can be applied both in stationary and nonstationary simulations. This dimensionally split scheme is unconditionally stable, monotone and compact, but not optimal in terms of numerical cross-diffusion. (An example of a genuine multidimensional BSBT scheme with a reduced amount of cross-diffusion is considered for unstructured grids as discussed in Section 8.3.) The BSBT scheme can also be derived from the method of characteristics and therefore, it can be considered as a semi-Lagrangian scheme. This scheme is therefore consistent with local wave characteristics. This view is the basic philosophy of the numerical approach used in SWAN and is related to the fact that the energy flux must remain constant along a wave characteristic (see, e.g., Whitham (1974), pg. 245 and Eq. (11.61), and Holthuijsen (2007), pg. 200). This makes sure that the scheme is flux conservative which is necessary for wave shoaling over rapidly varying sea beds.


SORDUP


For the SORDUP scheme, which is the default scheme for stationary computations, the second and third terms of Eq. (3.2) representing $x-$ and $y-$derivatives, respectively, are replaced by

  $$\left( \frac{3 (c_x N)_{i} - 4 (c_x N)_{i-1} + (c_x N)_{i-2}}{2 \Delta x} \right)^{n}_{j, l, m}$$ (3.7)

and

  $$\left( \frac{3 (c_y N)_{j} - 4 (c_y N)_{j-1} + (c_y N)_{j-2}}{2 \Delta y} \right)^{n}_{i, l, m}$$ (3.8)

in case of the first sweep, i.e. $c_x > 0$ and $c_y > 0$ (cf. Section 3.3). See also Rogers et al. (2002). This scheme, also known as a BDF scheme (Gear, 1971), is second accurate in space, but first order in time (not relevant), and is flux conservative (empirical evidence), but not monotone. In addition, this scheme preserves causality and produces less amount of numerical diffusion and is not significantly more expensive than the BSBT scheme.


In the neighboorhood of open boundaries, land boundaries and obstacles (i.e., the last two grids adjoining such grid points for the SORDUP scheme), SWAN will revert to the first order upwind BSBT scheme. This scheme has a larger numerical diffusion but that is usually acceptable over the small distances involved.


Stelling and Leendertse scheme


For the so-called cyclic scheme of Stelling and Leendertse (1992), which is the default scheme for nonstationary computations, the second and third terms of Eq. (3.2) representing $x-$ and $y-$derivatives, respectively, are replaced by

  $$\left( \frac{10(c_x N)_{i} - 15 (c_x N)_{i-1} + 6 (c_x N)_{i-2} -... ...left( \frac{(c_x N)_{i+1} - (c_x N)_{i-1}}{4 \Delta x} \right)^{n-1}_{j, l, m}$$ (3.9)

and

  $$\left( \frac{10 (c_y N)_{j} - 15 (c_y N)_{j-1} + 6 (c_y N)_{j-2} ... ...left( \frac{(c_y N)_{j+1} - (c_y N)_{j-1}}{4 \Delta y} \right)^{n-1}_{i, l, m}$$ (3.10)

in case of the first sweep, i.e. $c_x > 0$ and $c_y > 0$ (cf. Section 3.3). See also Rogers et al. (2002). This scheme is second accurate in time and space, unconditionally stable, preserves causality, and is flux conservative (empirical evidence). Moreover, the amount of numerical diffusion generated by this scheme is significantly much smaller than both BSBT and SORDUP schemes.


In the neighboorhood of open boundaries, land boundaries and obstacles (i.e., the last three grids adjoining such grid points for the Stelling and Leendertse scheme), SWAN will revert to the first order upwind BSBT scheme.


Although the Stelling and Leendertse scheme is unconditionally stable, i.e. there is no restriction on the chosen time step, there is, however, a practical limitation to this time step because of the spurious oscillations produced by this scheme. These so-called wiggles can appear much more strongly in the numerical solution if the associated Courant number is much larger than 1. Contrary to the BSBT scheme, this low-diffusion scheme is not able to suppress the wiggles sufficiently. The largest acceptable Courant number is 10 (based on the fastest wave) which is, however, a subjective criterion (Rogers et al., 2002).


Usually, the numerical diffusion of the Stelling and Leendertse scheme is so small that the so-called garden-sprinkler effect (GSE) may show up if propagation over very large distances is considered. This effect is due to the (coarse) spectral resolution (see Booij and Holthuijsen, 1987). It can be counteracted by anisotropic diffusion terms that have been explicitly added to the numerical scheme. Their values depend on the spectral resolution and the propagation time of the waves.


The diffusion applied in the propagation direction (locally along the wave ray) is

  $$D_{ss} = \frac{\Delta c_g^2T}{12}$$ (3.11)

where $\Delta c_g$ is the difference in the group velocities of successive frequencies, and $T$ is the so-called wave age, i.e. the time elapsed since the propagated swell was generated by the storm. The diffusion normal to the propagation direction (locally along the wave crest) is

  $$D_{nn} = \frac{c_g^2 \Delta \theta^2 T}{12}$$ (3.12)

From these, diffusion coefficients in Cartesian coordinates are calculated as

$$D_{xx} = D_{ss} \cos^2\theta + D_{nn} \sin^2 \theta \,,$$

$$D_{yy} = D_{ss} \sin^2\theta + D_{nn} \cos^2 \theta \,,$$

$$D_{xy} = (D_{ss}-D_{nn}) \cos \theta \sin \theta$$

The diffusion terms to be added to Eq. (3.1) are given by

  $$-D_{xx}\frac{\partial^2 N}{\partial x^2} -2D_{xy}\frac{\partial^2 N}{\partial x \partial y} -D_{yy}\frac{\partial^2 N}{\partial y^2}$$ (3.13)

Each of these second derivative terms is approximated in space using second order central differences, at time level $n-1$, as follows

  $$D_{xx} \left( \frac{N_{i+1} - 2N_{i} + N_{i-1}}{\Delta x^2} \right)^{n-1}_{j, l, m}$$ (3.14)

  $$D_{yy} \left( \frac{N_{j+1} - 2N_{j} + N_{j-1}}{\Delta y^2} \right)^{n-1}_{i, l, m}$$ (3.15)

  $$D_{xy} \left( \frac{N_{i,j} - N_{i-1,j} - N_{i,j-1} + N_{i-1,j-1}}{\Delta x \Delta y} \right)^{n-1}_{l, m}$$ (3.16)

This explicit scheme is fast (having little impact on computation time) but only conditionally stable. Through mathematical analysis (not shown) it can be shown that a likely stability condition is (Rogers et al., 2002)

  $$Q = \frac{\max(D_{xx},D_{yy},D_{xy})\Delta t}{\min(\Delta x,\Delta y)^2}\, \leq \, \frac{1}{2}$$ (3.17)

Thus, it is credible that Eq. (3.17) holds true for the two-dimensional Stelling and Leendertse scheme with this GSE correction. In short, by adding the GSE correction, the unconditionally stable advection scheme of SWAN becomes a conditionally stable advection-diffusion scheme. This restriction appears not to be severe as is commonly believed for a typical advection-diffusion equation. In experiments, it was found that with $Q \leq 0.48$, no instability was observed (Rogers et al., 2002). It is readily shown that for typical ocean applications $D_{nn}$ dominates the diffusion and $Q$ can be written as

  $$Q = \frac{{c_g}^2 T \Delta t \Delta \theta^2}{12 \Delta x^2}$$ (3.18)

The variable wave age $T$ could be computed during the computations of SWAN but it requires the same order of magnitude of computer memory as integrating the action balance equation. Instead a constant wave age $\overline{T}$ can be used as an approximation, so that Eq. (3.18) becomes

  $$Q = \frac{{\overline L} \mu \Delta \theta^2}{12 \Delta x}$$ (3.19)

where the characteristic travel distance of the waves is ${\overline L} = {c_g}{\overline T}$ (e.g., the dimension of the ocean basin) and $\mu = c_g \Delta t/\Delta x$ is the Courant number. For oceanic applications, the Courant number is typically $\mu \approx \frac{1}{2}$ so that $Q \leq \frac{1}{4}$ for typical values of $\Delta \theta \sim 10^o$ and ${\overline L}/\Delta x \sim 200$ (the number of grid points in one direction of the grid). This implies that the Stelling and Leendertse scheme with the GSE correction is stable for typical ocean cases. For shelf sea (regional) applications, the value of $\mu = {\cal O} (1)$ but the garden-sprinkler effect tends to be small on these scales and the diffusion can and should not be used to avoid the stability problem. For small-scale (local) applications, typically $\mu = {\cal O} (10-100)$. But such cases are usually treated as stationary and the SORDUP scheme (no GSE correction is included in this scheme), or preferably the BSBT scheme, should be used. See also Rogers et al. (2002) for further details.