Discretization procedure

For the sake of clarity of the algorithm description below, we put all the terms but the time derivative and propagation term in the geographical space of Eq. (2.16) in one term $F(\vec{x},\sigma,\theta)$:

  $$\frac{\partial N}{\partial t} + \nabla_{\vec{x}} \cdot [\vec{c}_{\vec{x}} N] = F$$ (8.5)

with $\vec{c}_{\vec{x}} = {\vec{c}}_g + \vec{u}$ the geographic velocity vector.


For the time being, we restrict ourselves to triangular meshes. However, other type of meshes can be employed as well, e.g. hybrid grids (consisting of both triangles and quadrilaterals). We consider a triangulation of a geographical domain in which Eq. (8.5) is solved; see Fig 8.1.

Figure 8.1: An example of triangulation.

Every vertex and all the triangles around this vertex are taken into account. Observe that the number of cells around a vertex can be different for all vertices. A vertex-based scheme is used in which the wave action $N$ is stored at the vertices and Eq. (8.5) is solved in each vertex. We note that the values at boundary vertices are fixed during the computation.


For the time integration, we adopt the first order implicit Euler scheme, as follows

  $$\frac{N^n - N^{n-1}}{\Delta t} + \nabla_{\vec{x}} \cdot [\vec{c}_{\vec{x}} N^n] = F^n$$ (8.6)

where $\Delta t$ is the time step and $n$ is the time step counter. The main property of this approximation is that it does not suffer from the stability restriction imposed by the CFL condition inherent in the explicit methods as employed in most spectral models. In principle, the time step is limited only by the desired temporal accuracy. This procedure, however, involves the solution of a large system of equations.


A point-by-point multi-directional Gauss-Seidel iteration technique is employed for updating all grid vertices. A key feature of this technique is that it takes advantage of the newly acquired vertex values during an iteration. It is locally implicit but globally explicit. In other words, it circumvents the need to build or store large matrices and remains stable at any time steps. This means that this numerical procedure can converge to steady state much more rapidly than explicit methods without requiring too much computational work and memory as do implicit methods.


We consider the update of a vertex as labeled 1 in Figure 8.2. This involves looping over each cell of this vertex.

Figure 8.2: Update of the wave action at vertex 1 in a triangle $\triangle$123 and the shaded directional sector in spectral space for which the waves are propagated.

We want to find an approximation for the propagation term of Eq. (8.6). To this end, we employ some vector calculus. We consider a triangular cell as depicted in Figure 8.3.

Figure 8.3: A triangular cell with geometrical quantities used for discretization in geographical space. Definitions of these quantities are provided in the text.

In vertex 1, we apply a mapping from a local coordinate system $\vec{\xi}$ = ($\xi$,$\eta$) to the Cartesian one $\vec{x}$ = ($x$,$y$). Based on this transformation $\vec{x}(\vec{\xi})$, we have the following base vectors that are tangential to the coordinate lines $\xi$ and $\eta$, respectively,

  $${\vec{e}}_{(1)} = \frac{\partial \vec{x}}{\partial \xi} \, , \quad {\vec{e}}_{(2)} = \frac{\partial \vec{x}}{\partial \eta} \, .$$ (8.7)

The vectors

  $${\vec{e}}^{(1)} = \mbox{grad} \, \xi \, , \quad {\vec{e}}^{(2)} = \mbox{grad} \, \eta$$ (8.8)

are normal to the coordinate surface of constant $\xi$ and $\eta$, respectively (see Figure 8.3). Moreover, they are reciprocal to the base vectors, i.e.

  $${\vec{e}}_{(\alpha)} \cdot {\vec{e}}^{(\beta)} = \delta_{\alpha}^{\beta} \, , \quad \alpha, \beta = \{1,2\} \, ,$$ (8.9)

where $\delta_{\alpha}^{\beta}$ is Kronecker delta (which is unity if $\alpha$ = $\beta$, and zero otherwise). Using Cramer's rule, one can find

  $$\vec{e}^{(1)} = \frac{1}{D} ( e^2_{(2)},-e^1_{(2)} )^{\top}\, ,\,... ...1_{(1)} )^{\top}\, , \, \, D = e^2_{(2)} e^1_{(1)} - e^2_{(1)} e^1_{(2)} \, .$$ (8.10)

Next, we expand the propagation term of Eq. (8.6):

  $$\nabla_{\vec{x}} \cdot [\vec{c}_{\vec{x}} N] = \frac{\partial c_x N}{\partial x} + \frac{\partial c_y N}{\partial y} \, ,$$ (8.11)

where $c_x$ and $c_y$ are the $x-$ and $y-$components of the wave propagation vector $\vec{c}_{\vec{x}}$, respectively. Using the chain rule, we obtain

  $$\nabla_{\vec{x}} \cdot [\vec{c}_{\vec{x}} N] = e^{(1)}_1 \frac{\p... ...ial c_y N}{\partial \xi} + e^{(2)}_2 \frac{\partial c_y N}{\partial \eta} \, .$$ (8.12)

Further, we approximate the derivatives in Eq. (8.12). The most simplest one is a one-sided first order difference scheme, as follows

$$\frac{\partial c_x N}{\partial \xi} \approx \frac{c_x N_1 - c_x N... ...rtial c_x N}{\partial \eta} \approx \frac{c_x N_1 - c_x N_3 }{\Delta \eta} \, ,$$

$$\frac{\partial c_y N}{\partial \xi} \approx \frac{c_y N_1 - c_y N... ...rtial c_y N}{\partial \eta} \approx \frac{c_y N_1 - c_y N_3 }{\Delta \eta} \, ,$$ (8.13)

where the action densities at vertices 1, 2 and 3 are denoted by $N_1$, $N_2$ and $N_3$, respectively. Here, we choose the mapping $\vec{x}(\vec{\xi})$ such that $\Delta \xi$ = $\Delta \eta$ = 1. The approximation is completed by substituting (8.13) in (8.12):

  $$\nabla_{\vec{x}} \cdot [\vec{c}_{\vec{x}} N] \approx c_x N\vert _... ...rt _3^1 e^{(2)}_1 + c_y N\vert _2^1 e^{(1)}_2 + c_y N\vert _3^1 e^{(2)}_2 \, .$$ (8.14)

Note that the components of the vectors ${\vec{e}}^{(1)}$ and ${\vec{e}}^{(2)}$ in Eq. (8.14) are given by Eqs. (8.10), while the base vectors are calculated according to

  $${\vec{e}}_{(1)} = {\vec{x}}_1 - {\vec{x}}_2 \, , \quad {\vec{e}}_{(2)} = {\vec{x}}_1 - {\vec{x}}_3$$ (8.15)

with ${\vec{x}}_i = (x_i, y_i)$ the position vector of vertex $i$ in a Cartesian coordinate system. This space discretization is of lowest order accurate and conserves action (see Section 8.7).


The upwind difference scheme BSBT (8.14) is employed for three reasons. First, it is compact, i.e. operating on one triangle only. Second, it enforces the propagation of wave action to follow the characteristics. Hence, this scheme produces numerical diffusion along these characteristics, and as such, it is multidimensional with a minimum amount of cross-diffusion. Third, it is monotone, i.e. guaranteeing $N > 0$ everywhere. This finite difference BSBT scheme appears to be identical to the well-known N(arrow) scheme, based on the fluctuation splitting method, on regular, triangular grids obtained from triangulation of rectangles with diagonals aligning with the wave characteristic as much as possible (see Struijs, 1994, pp. 63-64). Hence, the BSBT scheme share with the N scheme in many respects: it is multidimensional, first order accurate, optimal in terms of cross-diffusion, monotone, conservative, narrow stencil and consistent with local wave characteristics (see Section 3.2.1), but they are not identical in general.


Given the action densities $N^n_2$ and $N^n_3$ at vertices 2 and 3 of triangle $\triangle$123, the wave action in vertex 1 is readily determined according to

$$\left[ \frac{1}{\Delta t} + c_{x,1} \left( e^{(1)}_1 + e^{(2)}_1 \right) + c_{y,1} \left( e^{(1)}_2 + e^{(2)}_2 \right) \right] N_1^n =$$

$$\frac{N_1^{n-1}}{\Delta t}+\left( c_{x,2} e^{(1)}_1 + c_{y,2} e^{... ...) N^n_2 + \left( c_{x,3} e^{(2)}_1 + c_{y,3} e^{(2)}_2 \right) N^n_3 + F^n \, .$$ (8.16)

The wave directions between faces $\vec{e}_{(1)}$ and $\vec{e}_{(2)}$ enclose all wave energy propagation in between the corresponding directions $\theta_1$ and $\theta_2$ as indicated as a shaded sector in Figure 8.2. This sector is the domain of dependence of Eq. (8.16) in vertex 1. Since, the wave characteristics lie within this directional sector, this ensures that the CFL number used will properly capture the propagation of wave action towards vertex 1. So, propagation is in line with the causality principle and is not subjected to a CFL stability criterion. Next, the term $F^n$ in Eq. (8.16) is discretized implicitly in the sector considered. Since the approximation in the spectral space and the linearization of the source terms are well explained in Section 3.3, we shall not pursue them any further. Eq. (8.16) constitute a coupled set of linear, algebraic equations for all spectral bins within the sector considered at vertex 1. The solution is found by means of an iterative solver; see Section 3.3 for details.


The update of vertex 1 is completed when all surrounding cells have been treated. This allows waves to transmit from all directions. Due to refraction and nonlinear interactions, wave energy shifts in the spectral space from one directional sector to another. This is taken into account properly by repeating the whole procedure with converging results.