[orcid=0000-0002-2164-2884]

[1]

url]https://www.researchgate.net/profile/Mark-Loveland-2

Methodology, Software, Validation, Writing - Original Draft

1]organization=Oden Institute for Computational Engineering and Sciences,
University of Texas at Austin, addressline=Oden Institute for Computational
Engineering and Sciences, The University of Texas at Austin. 201 E. 24th St.,
city=Austin, state=Texas, postcode=78712, country=USA

[orcid=0000-0001-6940-4191]

Methodology, Writing - Review & Editing

2]organization=The Department of Data Science, The Norwegian University of
Life Science, addressline=Drøbakveien 31, city=Ås, postcode=1433,
country=Norway 3]organization=Department of Scientific Computing and Numerical
Analysis, Simula Research Laboratory, addressline=Kristian Augusts gate 23,
city=Oslo, postcode=0164, country=Norway Methodology, Writing - Review &
Editing 4]organization=National Oceanic and Atmospheric Administration (NOAA),
NWS, NCEP, EMC, addressline=830 University Research Court, city=College Park,
state=MD, postcode=20740, country=USA

[orcid=0000-0001-7273-0684] Resources, Supervision, Project administration
[1]Corresponding author

# WAVEx: Stabilized Finite Elements for Spectral Wind Wave Models Using
FEniCSx

Mark Loveland<EMAIL_ADDRESS>[  Eirik Valseth  [  [  Jessica Meixner
[  Clint Dawson

###### Abstract

Several potential FEM discretizations of the Wave Action Balance Equation are
discussed. The methods, which include streamline upwind Petrov-Galerkin
(SUPG), least squares, and discontinuous Galerkin, are implemented using the
open source finite element library FEniCSx for simplified 2-D cases. Open
source finite element libraries, such as FEniCSx, typically only support
geometries up to dimension of 3. The Wave Action Balance Equation is 4
dimensions in space so this presents difficulties. A method to use a FEM
library, such as FEniCSx, to solve problems in domains with dimension larger
than 4 using the product basis is discussed. A new spectral wind wave model,
WAVEx, is formulated and implemented using the new finite element library
FEniCSx. WAVEx is designed to allow for construction of multiple FEM
discretizations with relatively small modifications in the Python code base.
An example implementation is then demonstrated with WAVEx using continuous
finite elements and SUPG stabilization in geographic/spectral space. For
propagation in time, a generalized one step implicit finite difference method
is used. When source terms are active, the second order operator splitting
scheme known as Strang splitting is used. In the splitting scheme, propagation
is solved using the aforementioned implicit method and the nonlinear source
terms are treated explicitly using second order Runge-Kutta. Several test
cases which are part of the Office for Naval Research Test Bed (ONR Test Bed)
are demonstrated both with and without 3rd generation source terms and results
are compared to analytic solutions, observations, and SWAN output.

###### keywords:

Spectral wind wave model Stabilized finite element methods FEniCSx Python

Several finite element methods are investigated for solving the Wave Action
Balance Equation (WAE) on simplified 2D problems

A finite element solver using SUPG stabilization, implicit time stepping, and
Strang splitting for the full 4D WAE is implemented using FEniCSx

A set of six test cases from the Office of Naval Research test bed are run
using FEniCSx on unstructured meshes and accuracy is compared with analytic
solutions and SWAN output

##  1 Introduction

Wind waves are those which typically range in period from 0.25-30 seconds, or
wave lengths of 0.1-1500 meters, and are the product of disturbances of the
ocean surface due to locally generated winds [1]. Wind wave models take into
consideration the local winds, bathymetry, and currents in order to
approximate how wind waves evolve over time.

This research focuses specifically on spectral wind wave models which are
defined through the Wave Action Balance Equation (WAE). These spectral wind
wave models work by estimating the transport of wave action distributed among
wind waves of various frequencies and directions as defined through the WAE
[2, 3, 4, 5, 6]. The WAE presents challenges for numerical discretization
techniques because it is a first order transient hyperbolic equation whose
domain lies in 4 dimensions (2 in space, frequency, direction) with a non-
divergent free velocity field and possibly nonlinear source/sink terms.
Numerical spectral wind wave models have been in operation since the 1970s and
are used daily as a tool in many fields including engineering, shipping, and
meteorology [7, 8, 2]. They are also important tools for both improving our
understanding of ocean processes and for predicting the impacts of severe
weather events. The wind wave models can produce results regarding the state
of the sea surface (sea state) such as average wave height, significant wave
height, mean weave direction, peak period, and wave radiation stresses.

Currently, there are many well-known spectral wind wave models used in
practice such as ECWAM [7, 9] which is supported by the European Centre for
Medium-Range Weather Forecasts (ECMWF), Simulating WAves Nearshore (SWAN) [10]
which is supported by Delft University of Technology, and WAVEWATCH III [8]
which is supported by the United States National Oceanic and Atmospheric
Administration (NOAA). A brief summary of the methods employed in contemporary
numerical spectral follows but we note that there are many great historical
reviews on numerical spectral wind wave models see e.g, [11, 12, 2, 13, 3, 5,
4, 14].

All numerical spectral wind wave models rely on numerical approximations to
the governing partial differential equation (PDE), i.e., the WAE.
Historically, spectral wind wave models used finite difference approximations
along with splitting schemes to separate the propagation and source terms of
the WAE [13]. For instance, the ECWAM model contains both a first order
explicit upwind scheme and a second order leapfrog method for the propagation
terms with a semi-implicit second order centered difference for the source
terms [7]. The original WAVEWATCH III model uses an explicit predictor-
corrector method which approximates Crank-Nicholson and an explicit Euler
scheme for the source terms [8]. The scheme of SWAN was unique at its time of
inception in that it used a fully implicit solver without splitting nonlinear
and linear terms using the ideas of a nonlinear solver by Patankar [15]. In
SWAN, a finite difference formulation based on an implicit upwind scheme is
used [10]. This idea of using an implicit time step scheme is relatively
uncommon in fluid problems such as the WAE but show great promise in
applications where the CFL condition is particularly restrictive [14].

It is also common when splitting the WAE into propagation and source to employ
a method of characteristics for the propagation portion [14]. This is an
intriguing method since without source terms the WAE is just an advection
problem and the characteristics should give the exact solution. However, it
can be difficult to use these methods for complex bathymetries and complicated
source terms and thus numerical approximations are necessary in practice. Two
popular models, STWAVE and TOMAWAC, were developed beginning in the 1990s use
this characteristic paradigm and in many cases with great success [16, 17].

As the reliability and usage of wave models increased, there was a push for
spectral wind wave models to be used in more shallow waters. Historically,
structured grids are common because they are easier to implement and may be
computationally less expensive than unstructured grids. However, in order for
structured grids to work for coastal applications, very high levels of
refinement are needed to capture often highly irregular changes in bathymetry
[14] as well as coastline. Thus, spectral wind wave models capable of
unstructured grids became more popular since they are often computationally
cheaper than refining a structured grid multiple levels.

The first known model capable of computaions on unstructured meshes is the
aforementioned TOMAWAC [17], which was released in 1997 and was able to
represent the geographic domain with unstructured triangles using what is
often known as a characteristic-Galerkin or a Lagrange-Galerkin approach. That
is, a method which uses finite element interpolation in combination with a
method of characteristics. In [18], Hsu _et al._ developed a finite element
method (FEM) version of SWAN which used a Taylor-Galerkin scheme in geographic
space while splitting the spectral frequency, direction, and source terms.
This model was further developed, focusing on improvements in eliminating
oscillations around steep gradients, into wind wave model (WWM) I and WWM II
which added an implicit and explicit fluctuation splitting schemes (also known
as residual distribution schemes) to the geographic space [19, 20, 21, 22].
This same approach was added into the WAVEWATCH III model, extending WAVEWATCH
III to allow for unstructured meshes and implicit time steps [23]. Around the
same time, SWAN became available on unstructured meshes made of triangular
elements by extending its finite difference upwinding scheme [24].

In addition to the model of Hsu, there have been a couple more finite element
methods applied to the WAE. In 2012, Yildirim _et al._ used a hybrid approach
with discontinuous Galerkin in geographic space and Fourier collocation in
spectral space to create a third generation model [25]. In 2014, Meixner used
a discontinuous Galerkin FEM in spectral and geographic space [26]. Both of
the models used Runge-Kutta type explicit time stepping schemes.

In addition to models which use finite difference methods, and FEMs, there are
models which use finite volume methods. FVCOM-SWAVE is an implementation of a
finite volume method which allows for unstructured grids [27]. A summary of
all of the aforementioned schemes are contained in Table 1.

Model Name |  Year |  Numerical Method |  Method Description |  Grid Type
---|---|---|---|---
WAM |  1988 |  Finite Difference |  Advection: Explicit 1st order upwind, leapfrog Source: Centered difference (semi-implicit) |  Structured
WAVEWATCH III |  1990 |  Finite Difference |  Advection: QUICKEST (Explicit 3rd order upwind) Source: Semi-implicit |  Structured
STWAVE |  1990s |  Method of Characterstics/ Finite Difference |  Advection: Method of characteristics Source: Finite difference |  Structured
TOMAWAC |  1997 |  Method of Charactersitcs/FEM |  Advection: Method of characteristics Source: semi-implicit (linear finite elements for interpolation) |  Unstructured
SWAN |  1999, 2008 |  Finite Difference |  Space: Upwinding Time: First order implicit |  Unstructured
WWMII/WW3 |  2009 |  Contour residual scheme |  Space: contour residual distribution Time: Implicit upwind |  Unstructured
FVCOM-SWAVE |  2009 |  Finite Volume |  FVCOM finite volume in space, centered difference in directional, flux corrected transport algorithm for frequency |  Unstructured
Yildirim |  2011 |  FEM/Spectral |  Geographic space: Discontinuous Galerkin Spectral space: Fourier collocation Time: TVD Runge-Kutta |  Unstructured
Meixner |  2014 |  FEM |  Space: Discontinuous Galerkin Time: SSP Runge-Kutta |  Unstructured
Table 1: Summary of numerical schemes in contemporary spectral wind wave
models.

This work focuses on the development and implementation of a new spectral wind
wave model, called WAVEx, that is capable of using unstructured meshes via the
FEM along with implicit time stepping using the open source finite element
library, FEnicSx [28, 29]. FEMs gained popularity starting in the mid
twentieth century as a means to numerically approximate differential equations
see e.g., [30, 31, 32, 33]. In general a FEM consists of recasting a
differential equation into an equivalent weak form, partitioning the
associated domain into a finite number of subdomains (i.e., elements),
approximating the solution on each element with a finite dimensional basis,
and subsequently solving the associated linear algebraic system of equations
to achieve an approximate solution. The FEM has many benefits in that it is
backed by rigorous mathematical theory that readily extends to high order
approximations, unstructured meshes, and adaptive meshing see e.g., [34, 35,
30] for thorough introductions and complete discussions on advantages and
weaknesses.

FEMs were originally developed for elliptic type problems which often arise in
structural mechanics [35, 34]. However, the standard (Bubnov-Galerkin) FEMs
have issues when applied to many common fluid problems, which are often
hyperbolic in nature, and has historically led to unstable solutions and
solutions that had spurious oscillations [36, 34, 33, 37]. These problems have
been thoroughly analyzed and many stabilized finite element techniques were
developed [38, 39, 40] to address the issues. Some of the methods developed
which became popular include artificial viscosity [34], Taylor-Galerkin [41],
Streamlined Upwind Petrov-Galerkin (SUPG) [42], least squares FEM [43],
discontinuous Galerkin (DG) [44, 45, 46, 47], and Lagrange-Galerkin or
characteristic-Galerkin [48, 49, 50, 51].

In the following study, a set of four finite element formulations of the WAE:
SUPG, DG, least squares, and standard Bubnov-Galerkin, will be considered and
implemented using the FEniCSx library in a simplified two-dimensional domain
(1 in space and 1 in frequency) along with implicit finite difference
approximation in time for all FEM schemes. Two test cases will be
demonstrated, one of a propagating sinusoidal signal and one with a developing
shock. For each test case, the $L^{2}$ and $l^{\infty}$ errors are calculated
and studied with respect to $h$ refinement. Then, employing the strategies of
Loveland _et al._ in [52], WAVEx is implemented for the full four dimensional
WAE. Thus, WAVEx is a finite element model that employs SUPG stabilization
with a second-order implicit time-stepping scheme (Crank-Nicolson) along with
a Strang splitting approach to handle source terms in the WAE implemented
using FEniCSx. The open-source framework FEniCSx was chosen, as its powerful
open-source framework allows for automation of coding finite element solvers
and in hopes that it will allow for users to more easily experiment with
different numerical techniques or with different source terms via its Python
interface.

Several test cases from the US Office of Naval Research (ONR) testbed [53] are
reproduced in order to demonstrate the validity and applicability of WAVEx.
The collection of tests include theoretical problems with analytic solutions
as well as one lab experiment that has physical measurements. To measure
accuracy, the solution generated by WAVEx is compared to the tabulated results
at the given points. From these comparison points, both the root mean square
error (RMSE), i.e., $l^{2}$ and the $l^{\infty}$ error are computed. For each
test case, the parameters include significant wave height and mean wave
direction and are computed in WAVEx. To the knowledge of the authors, no
finite element methods have been used for the WAE that use SUPG stabilization,
least squares FEM, or discontinuous Galerkin approach along with implicit time
stepping.

##  2 Methods

###  2.1 Problem Definition

The Wave Action Balance Equation (WAE) is derived as a conservation law of
wave action, $N$, which is a scalar-valued function of horizontal geographical
space $(x,y)$ and spectral space $(\sigma,\theta)$, i.e., frequency and
direction, and can be found in e.g., [1, 6, 2, 54]. Assuming source terms are
linear with respect to $N$, the WAE is a linear, scalar-valued, hyperbolic
equation in 4 dimensions with a varying (both in 4-D space and time), non-
divergence free velocity field, ($\mathbf{c}(x,y,\sigma,\theta)$), which can
be determined independently of the unknown $N(x,y,\sigma,\theta)$. This is
different than most conservation laws such as Navier-Stokes and related
transport equations since typically the propagation velocity is just a
function of the unknown i.e. $\mathbf{c}(N)$.

The proper boundary and initial conditions are as follows, assuming we have a
bounded domain in four-dimensional space $\Omega\subset\mathbb{R}^{4}$ that is
sufficiently regular. Similar to other advection problems, the boundary is
split up into 2 segments, inflow and outflow:

$\begin{split}\Gamma_{-}={x\in\partial\Omega:\mathbf{c}\cdot\mathbf{n}<0}=\textrm{inflow}\\\
\Gamma_{+}={x\in\partial\Omega:\mathbf{c}\cdot\mathbf{n}\geq
0}=\textrm{outflow},\end{split}$ (1)

where $\mathbf{n}$ is the outwards unit normal vector to the boundary. Then it
has been shown that the following problem possesses a unique solution [55, 41,
39]:

$\begin{split}N_{t}+\nabla\cdot(\mathbf{c}N)=\frac{S(N,x,y,\sigma,\theta,t)}{\sigma}\quad\textrm{on}\quad\Omega\times(0,T),\\\
N=N_{-}\quad\mathrm{on}\quad\Gamma_{-},\\\
N=N_{0}\quad\textrm{on}\quad\Omega\quad\textrm{at}\quad t=0,\end{split}$ (2)

where $N_{t}$ denotes the (partial) time derivative of $N$, $S$ the
source/sink terms, $\Omega$ the computational domain, and $N_{-}$ the
specified essential boundary condition on the outflow boundary. In the case
when $S$ is either $0$ or is independent of $N$, a fairly straightforward
analytic solution to the above problem can be obtained via the method of
characteristics [41]. However, in practice $S$ is quite complex, often
nonlinear, and realistic problems become too complex to analyze with pencil
and paper and thus the need for numerical methods to approximately solve (2)
arises.

The propagation velocities, $\mathbf{c}$, are all determined using the
constitutive relation called the dispersion relation which is derived under
assumptions of Airy Wave Theory. The dispersion relation relates relative
radial frequency, $\sigma$, to wavenumber magnitude, $k$:

$\sigma^{2}=gk\tanh(kd).$ (3)

Here $g$ is the constant of gravitational acceleration and $d$ is total water
depth. The velocity $\mathbf{c}$ is a non-constant, non-divergence free vector
quantity defined as:

$\mathbf{c}=(c_{g}cos\theta+u,c_{g}sin\theta+v,c_{\sigma},c_{\theta}),$ (4)

where $u,v$ are the mean water velocities in $x,y$ directions respectively,
$c_{g}$ the relative group velocity (defined as
$\frac{\partial\sigma}{\partial k}$, whereas the absolute group velocity is
$\frac{\partial\omega}{\partial k}=c_{g}+\|(u,v)\|$), and $c_{\sigma}$,
$c_{\theta}$ define the advection of the spectra with respect to direction and
relative frequency, respectively. The group velocity, $c_{g}$, can be directly
obtained by differentiating the dispersion relation in (3) for relative
frequency, $\sigma$, with respect to wavenumber magnitude $k$. The explicit
expression is:

$c_{g}=\frac{1}{2}\left(1+\frac{2kd}{\sinh{2kd}}\right)\sqrt{\frac{g}{k}\tanh{(kd)}}.$
(5)

The propagation velocity $c_{\sigma}$, represents frequency shifting due to
changes in depths and currents. By applying the chain rule, we can find an
expression of $c_{\sigma}$:

$\begin{split}c_{\sigma}=\frac{d\sigma}{dt}=&\frac{k\sigma}{\sinh{(2kd)}}\left(\frac{\partial
d}{\partial t}+u\frac{\partial d}{\partial x}+v\frac{\partial d}{\partial
y}\right)-\\\ &c_{g}k\left(\frac{\partial u}{\partial
x}\text{cos}^{2}(\theta)+\frac{\partial u}{\partial
y}\text{cos}(\theta)\text{sin}(\theta)+\frac{\partial v}{\partial
x}\text{sin}(\theta)\text{cos}(\theta)+\frac{\partial v}{\partial
y}\text{cos}^{2}(\theta)\right).\end{split}$ (6)

The full derivation is quite lengthy but can be found in Appendix D in [1].
Now the propagation velocity $c_{\theta}$, with respect to $\theta$,
represents the shift in the spectrum due to refraction and diffraction.
Similarly to $c_{\sigma}$, the derivation is quite lengthy but can be obtained
by applying the chain rule:

$\begin{split}c_{\theta}=\frac{d\theta}{dt}=&\frac{\sigma}{\sinh{(2kd)}}\left(\frac{\partial
d}{\partial x}\text{sin}(\theta)-\frac{\partial d}{\partial
y}\text{cos}(\theta)\right)+\\\ &\frac{\partial u}{\partial
x}\text{cos}(\theta)\text{sin}(\theta)-\frac{\partial u}{\partial
y}\text{cos}^{2}(\theta)+\frac{\partial v}{\partial
x}\text{sin}^{2}(\theta)-\frac{\partial v}{\partial
y}\text{cos}(\theta)\text{sin}(\theta).\end{split}$ (7)

###  2.2 Source Terms

In the WAE (2), the source term $S(N,x,y,\sigma,\theta,t)$ was left arbitrary.
However, in practice, it can take many forms but in general it can be thought
of as a sum of three key sources/sinks:

$S=S_{in}+S_{diss}+S_{nl}.$ (8)

where $S_{in}$ represents any contribution to the spectrum due to wind input,
$S_{diss}$ represents any change in the energy spectrum due to dissipation
which can include things like whitecapping and surf zone breaking, and
$S_{nl}$ represents changes in the spectrum due to nonlinear wave
interactions. When source terms are active, the WAVEx model in this study will
use the default 3rd generation source term package used in SWAN cycle III
Version 41.41, see details in [56]. The default wind input term is that from
WAM Cycle III and first outlined in [57], the dissipation from whitecapping,
bottom friction, and depth-induced breaking are included and the DIA
approximation is used to account for nonlinear interactions. More details on
these source terms can be found in Chapter 3 of [58].

###  2.3 Weak Formulations For the Wave Action Balance Equation

In this section, several weak forms of the WAE will be presented. These weak
forms will eventually lead to to fully discrete systems that are then solved
with FEniCSx. First, for comparison purposes, two forms of the standard
Bubnov-Galerkin approach will be derived for the WAE and the problems with the
standard Bubnov-Galerkin discretizations will be discussed. Three alternatives
to the standard Bubnov-Galerkin approach: the least squares FEM, the SUPG
method, and the DG method will be derived and discussed. All of these weak
forms will then be run on two simplified test cases using FEniCSx in Section
3.1 in order to investigate convergence properties of all the stabilized
finite element methods with respect to $h$ refinement. The full discretization
of an SUPG method will be described for the full four-dimensional setting in
the subsequent section. Note that we utilize inner product notation for
integrals, i.e., $\left(u,v\right)_{\Omega}=\int_{\Omega}uv\,\rm{d}\Omega$.

####  2.3.1 Bubnov-Galerkin

Taking the strong form of the WAE as in (2), we obtain an equivalent problem
but in the weak form by multiplying both sides by a test function $w\in
L^{2}(\Omega)$ and integrating over a bounded and Lipschitz domain
$\Omega\subset\mathbb{R}^{4}$. Note that we do not integrate over time since
we will eventually employ finite difference approximation in this dimension.
Thus we have the following weak form which we will call CG Strong since there
are more regularity requirements on the trial space. Find $N\in U$ such that:

$\left(\frac{\partial N}{\partial
t},w\right)_{\Omega}+\left(\nabla\cdot(\textbf{c}N),w\right)_{\Omega}=\left(\frac{S}{\sigma},w\right)_{\Omega}\quad\forall
w\in L^{2}(\Omega).$ (9)

Thus, it is required that $U=\\{N\in
H^{1}({\Omega},t);N(\Gamma_{-},t)=N_{-};\frac{\partial N}{\partial t}\in
L^{2}(\Omega,t)\\}$, recalling that $\Gamma_{-}$ is the inflow boundary and
$N_{-}$ is the Dirichlet condition on that part of the boundary. We will also
assume that velocity vector $\mathbf{c}$ and source term $S$ is sufficiently
well behaved for the weak form to have meaning. It can be shown that (9) is
well-posed (in the sense that a solution exists, it is unique, and
continuously depends on the data) see e.g., the book of Ern and Guermond [34].
However, it is also known that if (9) is discretized using the standard
Bubnov-Galerkin approach then the discrete inf-sup condition depends on mesh
size, $h$, and therefore doesn’t provide optimal error estimates [34]. By
optimal, we mean for polynomial approximations of order $k$ the convergence of
the error in the $L^{2}$ sense to of the same order of convergence as the
underlying FE interpolants, usually $h^{k+1}$ (this classical result can be
found in e.g., [59, 35, 33, 37]). In practice this leads to oscillatory
solutions since only $L^{2}$ stability can be shown [38, 30]. Also, the
standard Bubnov-Galerkin approximation in this setting often leads to ill-
conditioned systems of equations and can even lead to divergence in the finite
element solution especially if the exact solution is non-smooth [30, 36].

We can also derive a second weak form if we integrate by parts, which we will
call the CG form, find $N\in U$:

$\left(\frac{\partial N}{\partial
t},w\right)_{\Omega}-\left(\textbf{c}N,\nabla
w\right)_{\Omega}+\left(\mathbf{c}N\cdot\mathbf{n},w\right)_{\Gamma_{+}}=\left(\frac{S}{\sigma},w\right)_{\Omega}\quad\forall
w\in H^{1}(\Omega).$ (10)

In order for the term on the boundary to have meaning we must require
additional regularity than just $N\in L^{2}(\Omega)$ and so the analysis is
not so straightforward but has been shown to be well posed for cases of
homogeneous boundary conditions and in cases with lower dimensionality [60].
When discretized, similar problems can be expected to that of the strong CG
from in (9). In order to mitigate the problems inherent in these standard
Bubnov-Galerkin discretizations, we will investigate a least squares method,
an SUPG method, and a DG method.

####  2.3.2 Least Squares

The least squares approach is fundamentally different to the Galerkin approach
in that the least squares method is posed as a minimization problem of a
linear differential operator. A good reference on the properties of least
squares can be found in the text of Bochev and Gunzburger [61]. To derive the
least squares form for the WAE, we start with the strong form as in (2), we
define the differential operator $\mathcal{L}$ and right hand side as in [61]
so that the above equation $\mathcal{L}$ would be defined:

$\mathcal{L}\bullet=\frac{\partial\bullet}{\partial
t}+\nabla\cdot(\textbf{c}\bullet).$ (11)

Then the least squares variational form then is defined as follows. Find $N\in
U$:

$(\mathcal{L}N,\mathcal{L}w)_{\Omega}=(\frac{S}{\sigma},\mathcal{L}w)_{\Omega}\quad\forall
w\in U.$ (12)

This problem is equivalent to minimizing the functional
$J(V)=\frac{1}{2}(\mathcal{L}(V),\mathcal{L}(V))_{\Omega}-(\mathcal{L}V,S/\sigma)_{\Omega}$.
Explicitly writing out the definition of $\mathcal{L}$ gives the formulation
for the WAE will be (Note: the time derivative on the test space will be
omitted since the time derivative will be discretized via finite difference
and not finite elements) , find $N\in U$:

$(\frac{\partial N}{\partial
t}+\nabla\cdot(\textbf{c}N),\nabla\cdot(\textbf{c}w))_{\Omega}=(\frac{S}{\sigma},\nabla\cdot(\textbf{c}w))_{\Omega}\quad\forall
w\in U,$ (13)

where $U=\\{N\in H^{1}(\Omega,t);N(\Gamma_{-},t)=N_{-};\frac{\partial
N}{\partial t}\in L^{2}(\Omega)\\}$, and again $\mathbf{c},S$ are assumed
sufficiently regular to give the integrals meaning. Now the least squares
problem generates a symmetric problem (trial and test spaces are the same) and
thus when discretized exhibits different properties than the Bubnov-Galerkin
approximations. It has been shown that for similar problems, such as the
advection-reaction equation, the least squares method can reduce oscillations
and provide stability [34]. It can be shown that least squares produces
optimal error convergence in the streamline directions however it has yet to
be shown that the least squares method produces optimal errors with respect to
$h$ refinement in the $L^{2}$ norm [34].

####  2.3.3 SUPG

The next stabilized form that will be investigated is the SUPG method. Since
the WAE is an advection dominated equation, the method is a natural choice to
employ. The method originally from [42] was designed for fluid problems and
has had much success in stabilizing advection dominated problems for fluids by
adding diffusion in the upwind direction in the weak formulation. The SUPG
scheme is formulated by adding a term with a vanishing residual to the
standard Bubnov-Galerkin form in the discrete section. This in turn modifies
the discrete test space which results in the change from a standard Bubnov-
Galerkin form to a so called Petrov-Galerkin form. We can add this residual
term to either CG strong as in (9) or the CG form as in (10). We will
demonstrate both possible SUPG formulations now. We start with CG strong form
and taking it into the discrete setting where the test and trial functions are
piecewise smooth and globally continuous. We will call the resulting form SUPG
strong which is, find $N_{h}\in U_{h}(\Omega)$:

$\left(\frac{\partial N_{h}}{\partial
t}+\nabla\cdot(\textbf{c}N_{h}),w_{h}\right)_{\Omega}+\left(\frac{\partial
N_{h}}{\partial
t}+\nabla\cdot(\textbf{c}N_{h})-\frac{S}{\sigma},\tau\mathbf{c}\cdot\nabla(w_{h})\right)_{\Omega_{e}}=\left(\frac{S}{\sigma},w_{h}\right)_{\Omega}\quad\forall
w_{h}\in V_{h}(\Omega).$ (14)

Here $\Omega_{e}$ is the interior over each finite element. Since in the
analytic setting the residual term is identically 0, then this problem
inherits the well-posedness from (9). The weak form can can be simplified by
adding the integrals together:

$\left(\frac{\partial N_{h}}{\partial
t}+\nabla\cdot(\textbf{c}N_{h}),w_{h}+\tau\mathbf{c}\cdot\nabla(w_{h})\right)_{\Omega_{e}}=\left(\frac{S}{\sigma},w_{h}+\tau\mathbf{c}\cdot\nabla(w_{h})\right)_{\Omega_{e}}\quad\forall
w_{h}\in V_{h}(\Omega).$ (15)

This will be called the SUPG strong form that will be discretized with
continuous finite elements in FEniCSx. It has been shown that adding this
diffusion in the streamline direction improves the stability of the discrete
problem by improving the discrete inf-sup condition (see [38, 42]) and
allowing for discrete stability in the $H^{1}$ sense. Furthermore, error
estimates for very similar problems to the WAE, such as the advection-reaction
equation show nearly optimal rates at $O(h^{k+0.5})$ [38, 62, 63].

We can also arrive at an SUPG form by taking (10) into the discrete setting
and adding the residual to get, find $N_{h}\in U_{h}(\Omega)$:

$\begin{split}\left(\frac{\partial N_{h}}{\partial
t},w_{h}\right)_{\Omega}-\left(\textbf{c}N_{h},\nabla
w_{h}\right)_{\Omega}+\left(\mathbf{c}N_{h}\cdot\mathbf{n},w_{h}\right)_{\Gamma_{+}}+\\\
\left(\frac{\partial N_{h}}{\partial
t}+\nabla\cdot(\textbf{c}N_{h})-\frac{S}{\sigma},\tau\mathbf{c}\cdot\nabla(w_{h})\right)_{\Omega_{e}}=\left(\frac{S}{\sigma},w_{h}\right)_{\Omega}\quad\forall
w_{h}\in V_{h}(\Omega).\end{split}$ (16)

After descritization of both the SUPG strong and SUPG forms, $\tau$ is taken
from [64, 62] and is chosen to be:

$\tau=\frac{h}{\|\mathbf{c}\|},$ (17)

where $h$ is the characteristic length of an element.

####  2.3.4 Discontinuous Galerkin

The next stabilized form that will be examined is the DG method as in the work
of Meixner [26]. The DG method is different to any of the other methods
discussed up to this point in that it allows for discontinuities between
element interfaces. The DG formulation is obtained by taking the strong form
of the WAE, multiplying by a test function, and integrating by parts over each
element. Find $N_{h}\in U_{h}$:

$\left(\frac{\partial N_{h}}{\partial
t},w_{h}\right)_{\Omega_{e}}-\left(\textbf{c}N_{h},\nabla
w_{h}\right)_{\Omega_{e}}+\left(\mathbf{c}\hat{N}_{h}\cdot\mathbf{n},w_{h}\right)_{\partial\Omega_{e}}=\left(\frac{S}{\sigma},w_{h}\right)_{\Omega_{e}}\quad\forall
w_{h}\in W_{h}.$ (18)

The test and trial spaces are virtually the same with $W_{h}=\\{v\in
L^{2}(\Omega):v|_{\Omega_{e}}\in P^{k}(\Omega_{e})\,\forall\Omega_{e}\in
T_{h}\\}$ and $U_{h}$ being $W_{h}$ but with a finite energy lift on the
inflow boundary. $P^{k}$ are polynomials of degree $k$, and $\Omega_{e}$
corresponds to the interior of an element while $T_{h}$ is s finite element
partition of the domain $\Omega$. $\hat{N}_{h}$ is a uniquely defined flux on
the boundary which is required to form a closed system of equations (see [26,
45, 44]). In this implementation, the fluxes on the element boundaries are
uniquely defined in this case using the Lax-Friedrichs flux which is described
in detail in the work of Meixner _et al._ 2014 [26]. The DG method has well
documented analysis on stability and error estimation for hyperbolic problems
such as the advection reaction equation, which is essentially what the Wave
Action Balance Equation is [45, 47, 46]. It has been shown that for a DG
scheme like the one above, stability in the $H^{1}$ sense exists and that
almost optimal convergence, $O(h^{k+\frac{1}{2}})$ can be shown [47].

###  2.4 WAVEx: Implementation of Wave Action Balance Equation Solver with
FEniCSx

We now introduce WAVEx which is capable of solving the WAE in its full
4-dimensional spatial domains using the techniques described in [52], and is
implemented with the finite element library FEniCSx. WAVEx in its current
state allows for SUPG stabilization along with Strang splitting to deal with
the nonlinear source terms. The codebase for WAVEx is publically available and
on Github at https://github.com/Markloveland/WAVEx. This section will explain
in detail the numerics currently available in WAVEx including the finite
element schemes as well as the time stepping schemes.

####  2.4.1 Semi-Discrete Form

This particular implementation in WAVEx will use a continuous FEM approach
along with SUPG stabilization. We first limit ourselves to cases where the
global domain is a Cartesian product between two subdomains (geographic space,
and spectral space). This means we can construct the product basis:

$N(x,y,\sigma,\theta)\approx\sum_{i=1}^{N}\sum_{j=1}^{M}N_{i,j}\phi_{i}(x,y)\psi_{j}(\sigma,\theta),$
(19)

where $N_{i,j}$ are the scalar weights of the interpolating functions.

As we know, the standard Bubnov-Galerkin approximation of advective equations
leads to ill conditioned discrete problems and as in Section 2.3, we will use
SUPG stabilization. SUPG was chosen since it performed well in our preliminary
2D studies (see Section 3), and the fact that it is also computationally
cheapest. We begin by first recalling the two possible weak forms as described
in Section 2.3, SUPG strong and SUPG. The latter being when integration by
parts is conducted and the former is when integrating by parts is omitted.
Both options are included in the WAVEx code repository. For the form without
integration by parts, SUPG strong, we begin with (15) and substitute in the
product basis to get the semi-discrete (i.e., discrete in space, infinite
dimensional in time) system:

$(\frac{\partial(\phi_{i}\psi_{j}N_{i,j})}{\partial
t}+\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}N_{i,j})-\phi_{i}\psi_{j}S_{i,j},\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=0.$
(20)

Since the source terms, $S$ in general are quite complex, we will only seek
the $L^{2}$ projections of their exact form as shown in (20). Analogously, the
SUPG weak form of the WAE (16) in the semidiscrete setting is defined by
substituting in our product basis into (16) to get:

$\begin{split}(\frac{\partial(\phi_{i}\psi_{j}N_{i,j})}{\partial
t},\gamma_{k}\beta_{l})_{\Omega}-(\mathbf{c}\phi_{i}\psi_{j}N_{i,j},\nabla(\gamma_{k}\beta_{l}))_{\Omega}+(\mathbf{c}\phi_{i}\psi_{j}\cdot\mathbf{n}N_{i,j},\gamma_{k}\beta_{l})_{\partial\Omega}+\\\
(\frac{\partial(\phi_{i}\psi_{j}N_{i,j})}{\partial
t}+\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}N_{i,j})-\phi_{i}\psi_{j}S_{i,j},\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=(\phi_{i}\psi_{j}S_{i,j},\gamma_{k}\beta_{l})_{\Omega}.\end{split}$
(21)

This takes care of the terms not related to the time step, the choice of
finite difference approximation in time will dictate the full discretization.

####  2.4.2 Time Step Choice and Fully Discrete Form

WAVEx uses a generalized one step finite difference approximation for the time
derivative in (20) and (21):

$\begin{split}(\frac{\partial(\phi_{i}\psi_{j}N_{i,j})}{\partial
t},\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=F(N,t)\rightarrow\\\
(\frac{\phi_{i}^{n+1}\psi_{j}^{n+1}N_{i,j}^{n+1}-\phi_{i}^{n}\psi_{j}^{n}N_{i,j}^{n}}{\Delta
t},\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=\alpha_{t}F(N,t^{n+1})+(1-\alpha_{t})F(N,t^{n}).\end{split}$
(22)

where term $F(N,t)$ is shorthand from (20) and equals:

$F(N,t)=-(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}N_{i,j})-\phi_{i}\psi_{j}S_{i,j},\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega},$
(23)

and for (21):

$\begin{split}F(N,t)=(\mathbf{c}\phi_{i}\psi_{j}N_{i,j},\nabla(\gamma_{k}\beta_{l}))_{\Omega}-(\mathbf{c}\phi_{i}\psi_{j}\cdot\mathbf{n}N_{i,j},\gamma_{k}\beta_{l})_{\partial\Omega}+(\phi_{i}\psi_{j}S_{i,j},\gamma_{k}\beta_{l})_{\Omega}\\\
-(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}N_{i,j})-\phi_{i}\psi_{j}S_{i,j},\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega},\end{split}$
(24)

the superscripts represent the discrete level in time, $n$ will be the
previous time step and $n+1$ is the point in time to be solved for, and
$\alpha_{t}$ is a time step parameter between $0$ to $1$. The scheme is
implicit unless $\alpha_{t}=0$, as $\alpha_{t}=0$ would be the first order
explicit Euler approximation, $\alpha_{t}=1$ would be the first order implicit
Euler approximation, and $\theta=0.5$ would be the second order Crank-Nicolson
scheme. The terms with superscript $n$ are known and move to the RHS while the
terms $n+1$ will be part of the global stiffness matrix at each time step.

The definition given for the finite difference approximation to the time
derivative leads to a fully discrete system described by (22) depending on the
choice of $F(N,t)$. Choosing (23) will yield the fully discrete system for
SUPG strong (no integration by parts) and choosing (24) will give the fully
discrete system for SUPG weak. Once a choice of finite basis is made for trial
and test spaces, the forms above will yield a system of equations that in turn
can be solved for the action balance at the next time step, $N^{n+1}_{i,j}$.
However, there are two difficulties that need to be addressed in order to
reach a final implementation: $i)$, FEniCS(x) can not evaluate these integrals
on the entire 4-dimensional domain $\Omega$. In the next section, the discrete
equations in (22) will be broken down into products of the geographic and
spectral subdomains using the product rule as in [52] into a discrete form
that is implementable in FEniCSx. $ii)$, in general the source term $S$ is
nonlinear with respect to $N$. This means that the discrete form of (22) will
not yield a system of linear equations as long as we are using an implicit
time step, $\alpha_{t}>0$. To avoid this, an operator splitting method will
described in Section 2.4.4 in order to avoid having to solve a nonlinear
system. Although in the future a fully implicit solver for the nonlinear
problem is of great interest.

####  2.4.3 FEniCSx Compatible Format

Now since we want to only evaluate functions in the 2-dimensional subdomains
with FEniCSx we need to take advantage of the product rule. As an example,
starting with the second term from (21) we get:

$-(\mathbf{c}\phi_{i}\psi_{j},\nabla(\gamma_{k}\beta_{l}))_{\Omega}=-(\mathbf{c}\phi_{i}\psi_{j},\gamma_{k}\nabla(\beta_{l})+\beta_{l}\nabla(\gamma_{k}))_{\Omega}.$
(25)

Similar to the previous examples, using the property that the basis functions
are only non-constant in their respective subdomains we get the following
simplifications ($\nabla_{1}$ denotes gradient operator in the geographic
subdomain, and $\nabla_{2}$ denotes the gradient operator in the spectral
subdomain, and same for velocity):

$-(\mathbf{c}\phi_{i}\psi_{j},\gamma_{k}\nabla(\beta_{l})+\beta_{l}\nabla(\gamma_{k}))_{\Omega}=-(\mathbf{c}_{1}\phi_{i}\psi_{j},\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega}-(\mathbf{c}_{2}\phi_{i}\psi_{j},\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega}.$
(26)

This importantly allows us to pull integrands apart in the following fashion
so that we can use FEniCSx:

$\begin{split}-(\mathbf{c}_{1}\phi_{i}\psi_{j},\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega}-(\mathbf{c}_{2}\phi_{i}\psi_{j},\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega}=-\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\phi_{i}\mathbf{c}_{1}\cdot\nabla_{1}\gamma_{k}dxdy\\\
-\int_{\Omega_{2}}\psi_{j}\nabla_{2}\beta_{l}\cdot\int_{\Omega_{1}}\mathbf{c}_{2}\phi_{i}\gamma_{k}dxdy.\end{split}$
(27)

Where (27) represents a term that can be evaluated using FEniCSx. Integrands
inside $\Omega_{1}$ each create a $N_{dof1}\times N_{dof1}$ sparse matrix and
to form the full discrete system we will need a matrix at each quadrature
point in the second subdomain. Each matrix will be different because of the
changes in $\mathbf{c}$ in both subdomains. The remainder of the terms in (22)
can be rewritten in this way via the product rule. The full process is the
same can be seen in detail in Appendix A.

####  2.4.4 Operator Splitting

Since the source terms are in general non-linear, in the interest of saving
computational complexity we will employ a splitting scheme to handle the
source terms of the WAE. More extensive notes on operator splitting can be
found in [65]. We will be employing the second order Strang splitting [66]
scheme which can be defined on a system of ODEs as follows:

$\begin{split}\frac{dX}{dt}+A(X,t)=0,\in(0,T),\\\ X(0)=X_{0}.\\\ \end{split}$
(28)

In our case $A$ will be the discrete finite element operator in 4D space and
$X$ will be the action balance at all of the degrees of freedom at a
particular point in time. The Strang splitting scheme splits up $A$ as a sum
of two operators $A_{1}$ and $A_{2}$ for $n\geq 0,X^{n}\rightarrow
X^{n+\frac{1}{2}}\rightarrow\hat{X}^{n+\frac{1}{2}}\rightarrow X^{n+1}$.
Solve:

$\begin{split}\frac{dX_{1}}{dt}+A_{1}(X_{1},t)=0,t\in(t^{n},t^{n}+\frac{1}{2}),\\\
X_{1}(t^{n})=X^{n},\\\ \end{split}$ (29)

then set $X^{n+\frac{1}{2}}=X_{1}(t^{n+\frac{1}{2}})$ and solve:

$\begin{split}\frac{dX_{2}}{dt}+A_{2}(X_{2},t)=0,t\in(t^{n},t^{n}+\Delta
t),\\\ X_{2}(t^{n})=X^{n+\frac{1}{2}},\\\ \end{split}$ (30)

then set $\hat{X}^{n+\frac{1}{2}}=X_{2}(\Delta t)$ and solve:

$\begin{split}\frac{dX_{1}}{dt}+A_{1}(X_{1},t)=0,t\in(t^{n+\frac{1}{2}},t^{n+1}),\\\
X_{1}(t^{n+\frac{1}{2}})=\hat{X}^{n+\frac{1}{2}}.\\\ \end{split}$ (31)

Which gives $X^{n+1}=X_{1}(t^{n+1})$. We will apply this second order Strang
splitting to the WAE by setting $A_{1}$ to be the source terms and $A_{2}$ as
all remaining terms. In particular, for the case of the semi-discretization in
(20), the operators will be:

$\begin{split}A_{1}=-(\phi_{i}\psi_{j}S_{i,j},\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega},\\\
A_{2}=(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}N_{i,j}),\gamma_{k}\beta_{l}+\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}.\end{split}$
(32)

We will implement the Strang splitting with an implicit-explicit scheme. The
time step associated with operator $A_{2}$ will be the scheme from Section
2.4.2, whereas the two time steps associated with operator $A_{1}$ will be an
explicit second order Runge-Kutta (RK2) scheme [67]. The RK2 used in WAVEx is
defined as:

$\begin{split}N^{(1)}=N^{n}+\Delta tA_{1}(N^{n},t^{n}),\\\
N^{n+1}=\frac{1}{2}N^{n}+\frac{1}{2}(N^{(1)}+\Delta
tA_{1}(N^{(1)},t^{n+1})).\end{split}$ (33)

To summarize, WAVEx utilizes SUPG stabilization in 4 dimensions and employs
Strang splitting to separate advection from the source terms. The source terms
are integrated in time with an explicit RK2 scheme while the advection part is
integrated in time with an implicit generalized one step finite difference.
Since WAVEx uses the FEniCSx framework, it is automatically parallelized via
MPI and the linear algebra backend that is used is the Portable, Extensible
Toolkit for Scientific Computation (PETSc) [68]. The linear solver that is
used in all the following cases is the default GMRES solver with block Jacobi
preconditioning. Also, WAVEx is designed to allow for the addition of
different discretizations in both time and space besides just the specific
choices made in this case.

##  3 Results

###  3.1 Convergence Test Cases

To investigate the asymptotic convergence properties of the stabilized methods
introduced in Section 2.3, two test cases will be considered. These test cases
will be run with reduced dimension (2 instead of 4 in space) and the velocity
fields will be simplified. For each case the error in the finite element
solution, $N_{h}$, relative to the exact solution, $N_{exact}$, will be
computed both in the $L^{2}(\Omega)$ and the $l^{\infty}(\Omega)$ norm which
are defined as:

$\begin{split}\|N_{h}-N_{exact}\|_{L^{2}}=\sqrt{\int_{\Omega}(N_{h}-N_{exact})^{2}dx},\\\
\|N_{h}-N_{exact}\|_{l^{\infty}}=\max_{i\in
n_{dof}}(|N_{h}(x_{i})-N_{exact}(x_{i})|).\end{split}$ (34)

####  3.1.1 Propagation of a Sinusoidal Wave

The first test case is a harmonic wave propagating through a square domain
over time. The goal of this test case is to evaluate the capabilities of the
different FEM schemes in a "pure advection" scenario by setting source, $S$,
to 0. Harmonic waves are propagating across 2D domain in the direction of the
velocity vector, $\mathbf{c}=(1,1)$, which is constant space and time. The
domain is a square, $\Omega=(0,10)\times(0,10)$. The exact solution is:

$N_{analytic}=sin(x-c_{x}t)+cos(y-c_{y}t).$ (35)

Which dictates the initial and boundary conditions:

$\begin{split}N_{0}(x,y)=sin(x)+cos(y),\\\
N(x,y,t)=N_{analytic}\quad\textrm{for }(x,y)\in\Gamma_{-},\end{split}$ (36)

and as a consequence of the velocity field the Dirichlet boundary,
$\partial\Omega_{d}$, is the bottom and left sides of the square domain. The
test case has a total simulation time of 5 units which is split up into 1000
time steps. The time step is set small to minimize error from the
approximation in time and focus on FEM error convergence with respect to $h$
refinement.

Each scheme introduced in Section 2.3 is discretized with uniform triangular
elements in space and order 1 Lagrange finite element functions for both trial
and test space. Each scheme is fully implicit in time with order 1 implicit
Euler finite difference approximation. The computational grid is refined in
space with a uniform spacings of $L/4$, $L/8$, $L/12$, $L/16$, $L/20$, $L/24$.
For each finite element scheme, the global $L^{2}$ and $l^{\infty}$ errors are
computed. The errors at the final time $T=5$ are plotted for both norms and
for all schemes in Figure 1.

Figure 1: $L^{2}$ and $l^{\infty}$ errors at T = 5 for all schemes for case 1.

For the smooth solution present in case 1, nearly optimal convergence rates of
$O(h^{k+1})$ were shown for the SUPG cases. The CG methods without
stabilization converged at a rate of $h^{2}$ until machine precision was met
at $L/20$ which brings down the average shown in Table 2 to 1.559.

Table 2:  Error convergence rates for case 1. Method | Avg. L2 rate | Avg. Linf rate
---|---|---
CG Strong | 1.559 | 1.445
CG | 1.559 | 1.445
SUPG strong | 1.988 | 2.01
SUPG | 1.988 | 2.01
Least Squares | 1.790 | 1.730
DG | 1.611 | 1.981

The errors in the SUPG converged at a better rate than it’s theoretical rate
of $O(h^{1.5})$ most likely due to the minimal amount of residual present in
the smooth solution. Least squares had the highest in magnitude of $L^{2}$
error and converged at a rate between $h^{1.5}$ and $h^{2}$. Lastly, we
observe that the DG method had the lowest magnitude in errors for stabilized
schemes and converged near the optimal theoretical rate of $h^{1.5}$.

####  3.1.2 Formation of a Shock Near the Boundary

This test case is an adaptation of a common numerical test for advection
dominated PDEs and is adapted from the work of Egger and Schöberl [69]. This
case has an analytic solution that begins as a smooth bump then shifts into a
shock towards the right boundary. The aim is to see how the stabilized schemes
handle ill-conditioned problems and formation of discontinuities. The domain
is the unit square $\Omega=(0,1)\times(0,1)$ with a uniform velocity field
$\textbf{c}=\\{1,0\\}$. The analytic solution is:

$N_{exact}=(-4(y-0.5)^{2}+1)\left[x+\frac{e^{\gamma\mathbf{c}_{x}x}-1}{1-e^{\gamma\mathbf{c}_{x}}}\right],$
(37)

where $\gamma=100t$, $\mathbf{c}_{x}=1$, and $t\in[0.01,1]$. The number of
time steps for each run is set to 1000. A consequence of this analytic
solution is that S should be:

$(-4(y-0.5)^{2}+1)\left[\frac{100e^{100t}(e^{100tx}-1)}{(1-e^{100t})^{2}}+\frac{100xe^{100tx}}{1-e^{100t}}\right]+(-4(y-0.5)^{2}+1)\left[1+\frac{100te^{\gamma\cdot\mathbf{c}_{x}\cdot
x}}{1-e^{\gamma\cdot\mathbf{c}_{x}}}\right].$ (38)

The domains are the same structure as in case 1, and the levels of h
refinement for this case are $L/4$, $L/8$, $L/12$, $L/16$, $L/20$, $L/24$,
$L/28$, $L/32$, $L/36$. The errors at final time $T=1$ are plotted for both
$L^{2}$ and $l^{\infty}$ norms are displayed in Figure 2.

Figure 2: $L^{2}$ and $l^{\infty}$ errors at T = 1 for all schemes for case 2.

Even with a small time step that is guaranteed stable, the solutions from the
CG weak forms both produced the infamous oscillations in the final solution. A
comparison in the solution at $T=1$ of CG and SUPG implementations are shown
in Figure 3.

Figure 3: Solution for case 2 on a 32x32 node grid, SUPG right CG left. Table 3:  Error convergence rates for case 2. Method | Avg. L2 rate | Avg. Linf rate
---|---|---
CG Strong | 1.287 | 1.270
CG | 1.287 | 1.270
SUPG strong | 1.948 | 1.461
SUPG | 1.795 | 1.461
Least Squares | 1.907 | 1.440
DG | 1.999 | 1.460

Consequently, the resulting error convergence rates in the CG implementations
were not optimal as in the former case case. Interestingly, the SUPG and SUPG
strong results did have differences as presented in Figure 2 with the SUPG
strong version performing best. All the stabilized schemes did not produce
visible oscillations at any level of refinement and all schemes converged at
equal to or higher than the best theoretical rates as presented in Table 3.
The DG method had the lowest error of all the stabilized schemes which isn’t
surprising since for the same polynomial order the DG method has nearly double
the degrees of freedom as the continuous implementations.

###  3.2 Full Scale Test Cases

In this section, 6 test cases from the ONR testbed will be implemented. For
each of the case the RMSE and $l^{\infty}$ errors will be computed with
respect to the WAVEx output and analytic solution or lab data depending on the
case. The RMSE is computed using the definition:

$RMSE(x,y)=\sqrt{\frac{\sum_{i=1}^{N_{d}}(x_{i}-y_{i})^{2}}{N_{d}}},$ (39)

where $x$ represents the set of parameters generated by the wave model, $y$ is
the set of parameters that come from either the analytic solution or
observation, and $N_{d}$ is the number of points the ONR test bed has data
for. Using the same nomenclature as in (39), the $l^{\infty}$ measure is
defined as:

$l^{\infty}(x,y)=\max_{1\leq i\leq N_{d}}\lvert x_{i}-y_{i}\rvert.$ (40)

In all cases, the data that is compared to the WAVEx model is either
significant wave height ($H_{s}$), or both $H_{s}$ and mean wave direction
($\theta_{mean}$). Significant wave height is defined as the highest $1/3$ of
recorded waves at a point, and mean wave direction is as it sounds. In WAVEx,
the aforementioned statistical values must be estimated. The significant wave
height is estimated by first calculating the zeroth moment of the spectrum,
$m_{0}$. In general the $i^{th}$ moment of the spectrum is defined as:

$m_{i}(x,y)=\int_{-\pi}^{\pi}\int^{\sigma_{max}}_{\sigma_{min}}\sigma^{i}E(x,y,\sigma,\theta)d\sigma
d\theta,$ (41)

and recall $E=N\sigma$. The significant wave height is subsequently estimated
by:

$H_{s}(x,y)\approx 4\sqrt{m_{0}}.$ (42)

The reasoning behind this being a decent approximation to the highest one
third of waves is not trivial but its assumptions and derivation can be found
in several books such as the [1, 54]. The mean wave direction in degrees is
computed as:

$\theta_{mean}(x,y)=\frac{180}{\pi}\frac{\int_{-\pi}^{\pi}\int^{\sigma_{max}}_{\sigma_{min}}\cos{(\theta)}E(x,y,\sigma,\theta)d\sigma
d\theta}{\int_{-\pi}^{\pi}\int^{\sigma_{max}}_{\sigma_{min}}sin(\theta)E(x,y,\sigma,\theta)d\sigma
d\theta}.$ (43)

In all the presented test cases, WAVEx is run using the SUPG strong scheme as
in (20) and the implicit time stepping uses the Crank-Nicolson scheme which is
when $\alpha_{t}=0.5$ from (22) for second order accuracy in time.

####  3.2.1 A21: Shoaling

The first test case is a situation where some waves are propagating from the
deep ocean directly into the coastline. The water depth in the domain is
uniformly decreasing from 20 meters to 0 meters with a slope of 0.005. The
waves are propagating perpendicular to the coastline and are approximately
monochromatic. The waves at the inflow boundary have a significant wave height
of 1 meter with a mean period of 10 seconds. All source terms are zero and
there are no currents in this case. The analytic solution for the significant
wave height can be computed applying the analytic solution to the advection
equation with nonuniform velocity combined with the definition of the velocity
field for the wave action balance equation and the definition of significant
wave height. The result is that the significant wave height should be:

$H_{s}(x)=H_{s}(0)\sqrt{\frac{c_{g}(0)}{c_{g}(x)}}.$ (44)

An unstructured mesh is used in the WAVEx simulations for cases A11 and A21
illustrated in Figure 4. The mesh has 7458 nodes comprising 13841 triangular
elements and has variable resolution with element size ranging from 800 meters
near on the bottom to 20 meters near the top.

Figure 4: The unstructured geographic mesh used for test cases A11, A21.

In WAVEx, the spectral grid is capable of being unstructured but for all of
the ONR cases it is structured with logarithmic spacing in frequency and
uniform in direction. An example is shown in Figure 5.

Figure 5: Logarithmic spaced spectral grid as used in cases A11, A21.

The time step is set to 50 seconds for all the cases, and WAVEx solves until
the iterative solver reaches a steady state. The significant wave height is
computed along the propagation direction and results are plotted in Figure 6.
The results in this figure also clearly demonstrate the effect of instability
in the Bubnov-Galerkin solution as its oscillations makes the solution very
poor. The computed RMSE error vs the analytic solution is 0.000695 m for WAVEx
with SUPG stabilization and 0.313 m without, while $l^{\infty}$ error is
0.0018 m with SUPG and .9916 m without.

Figure 6: Significant wave heights from case A21.

####  3.2.2 A11: Refraction

The second test case is similar to case A21 in that the bathymetry is the same
and the domain is identical. The only difference is that the waves on the
inflow boundary now approach the shoreline at a direction 30 degrees above the
direction perpendicular to the shoreline. This too has an analytic solution.
Again, no source terms are present and the currents are all 0. In addition to
significant wave height, the mean wave direction is measured and the analytic
solutions are:

$\begin{split}\theta(x)=arcsin(c(x)\frac{sin(30)}{c(0)}),\\\
H_{s}(x)=H_{s}(0)\sqrt{\frac{c_{g}(0)cos(30)}{c_{g}(x)cos(\theta(x))}}.\end{split}$
(45)

The significant wave height is plotted against the analytic solution in Figure
7 while the mean wave direction is plotted in Figure 8. In both figures, the
WAVEx and analytic results agree very well. The observed RMSE error in HS is
0.00261 m and RMSE error in degrees is 0.119 degrees while the $l^{\infty}$
error is 0.00865 m and 0.1946 degrees.

Figure 7: Significant wave height for case A11. Figure 8: Mean wave direction
for case A11.

####  3.2.3 A31-34: Currents

Cases A31 through A34 are situations in deep water where there are interacting
currents but no source terms. All simulations are performed on the same domain
which is a rectangle shape of constant arbitrarily deep water (set to 10000m)
and is 4000 meters long in the principle direction. To avoid any noise from
the side boundaries, the width of the domain is set to 10000 meters. Each case
has either different currents or different mean wave directions on the inflow
boundary. Case A31 has an inflow boundary where the mean wave direction is
parallel with the principle axis with a current that is also parallel with the
principle axis and increases from 0 m/s at the inflow boundary to 2 m/s at the
outflow boundary. Case A32 has the same boundary conditions as A31 but the
currents in this case decrease from 0 m/s at the inflow boundary to -2 m/s at
the outflow boundary. Case A33 has the same current profile as case A31 but
the inflow boundary condition now has a mean wave direction of 120 degrees
relative to the positive direction of the principle axis. Finally, case A34 is
the same as A33 but instead of 120 degrees, it is set to 60 degrees. All waves
at the inflow boundary have a significant wave height of 1 meter and peak
period of 10 seconds. Similar to cases A11 and A21, an analytic expression for
the directions and significant wave heights may be derived. These are more
complex due to the nonzero currents and are different for cases A31, A32 and
the A33, A34 pair. The solutions are different since the former have incoming
spectra parallel to the currents while the latter are at an angle. The
analytic solution for significant wave height in cases A31 and A32 is:

$\begin{split}H_{s}(x)=H_{s}(0)\sqrt{\frac{c(0)^{2}}{c(x)(c(x)+2u(x))}},\\\
c(x)=c(0)(0.5+0.5\sqrt{1+4\frac{u(x)}{c(0)}}),\end{split}$ (46)

while the solution for cases A33 and A34 along with the analytic solution for
mean wave direction is:

$\begin{split}H_{s}(x)=H_{s}(0)\sqrt{\frac{sin(2\theta(0)}{sin(2\theta(x))}},\\\
\theta(x)=arccos(\frac{gk(0)cos(\theta(0))}{(\frac{2\pi}{Pkper}-u(x)k(0)cos(\theta(0)))^{2}}).\end{split}$
(47)

For these cases, the mesh is uniformly spaced with triangular elements. The
mesh contains 10,201 nodes which comprise 20,000 elements. This mesh is 10000
meters in the x direction and 4000 meters in the y direction. Each element has
a size of around 110 meters. The mesh is 100 elements wide and 100 elements
tall. The mesh is illustrated in Figure 9.

Figure 9: Geographic mesh for cases A31-A34.

The resulting significant wave heights from case A31 and A32 are shown in
Figures 10,11, and the computed RMSE are 0.000255 m, and 0.00109 m
respectively while $l^{\infty}$ error is 0.0003737 m and 0.00257 m.
Significant wave height for cases A33 and A34 are plotted in Figure 12 and
Figure 14. The RMSE for significant wave height for A33 is 0.000315 m, while
A34 is 0.000899 m. The $l^{\infty}$ for case A33 is 0.000659 m while A34 is
0.00128 m. The resulting mean wave direction for cases A33 and A34 are shown
in Figure 13 and Figure 15. The RMSE for mean wave direction in case A33 is
0.1793 degrees while for case A34 it is 0.436. The $l^{\infty}$ for case A33
is 0.028 degrees while case A34 is 0.0436 degrees. We observe close agreement
of WAVEx and the analytic solutions in all cases. This agreement is further
witnessed in the low errors in both mean wave direction and significant wave
height.

Figure 10: Significant wave height from case A31. Figure 11: Significant wave
height from case A32. Figure 12: Significant wave height from case A33. Figure
13: Mean wave direction from case A33. Figure 14: Significant wave height from
case A34. Figure 15: Mean wave direction from case A34.

####  3.2.4 L11: Wave Breaking on a Beach

The final test case is the L11 case, a laboratory experiment intended to
simulate groups of waves breaking on a beach. The bathymetry varies in slope
as it approaches the beach but is uniform in the along-shore direction. The
input spectrum is a JONSWAP spectrum with peakedness parameters specified in
the manual [53]. The mesh that is used is in this case and is set as the
rectangle $[0,113]\times[7.4,30]$, it is a structured, uniform mesh with
triangular elements. Each element is approximately 1 meter wide and there are
4681 total nodes comprising 9000 elements. The mesh is shown in Figure 16.

Figure 16: Geographic mesh for case L11.

In addition to comparison to the measured results from the ONR testbed [53],
we also compare the WAVEx results to results from SWAN. The SWAN numerical
model is set up with the third generation source terms are turned on, however
there is no wind in this case. The significant wave heights are recorded both
from SWAN and WAVEx and plotted against the observations in Figure 17 (note
that SWAN was run in 1-D mode). The RMSE in significant wave height with
respect to the observations were 0.0064m for WAVEx and 0.00556m for SWAN while
the $l^{\infty}$ errors were 0.0179m and 0.015m respectively.

Figure 17: Significant wave height for case L11.

##  4 Conclusions

The purpose of this work was to develop a finite element spectral wind wave
model using the open source FEniCSx library, which we call WAVEx. The
motivation behind this was to create a platform so that different numerical
schemes and source terms can easily be implemented and tested in the future.
First, some stabilized FEM methods were introduced (least squares, SUPG, DG)
and convergence rates with respect to mesh refinement were verified in simple
2D test cases. All of the stabilized schemes demonstrated convergence in the
$L^{2}$ and $l^{\infty}$ norms with respect to $h$ refinement while also
reducing oscillations when compared to the standard Bubnov-Galerkin CG
approach which produced oscillations when the solution became less smooth.

Next, a spectral wave model with SUPG stabilization in 4 dimensions was
implemented. The WAE was split into propagation and source terms via Strang
splitting for computational efficiency. The propagation portion was
implemented with an implicit time stepping scheme while the (nonlinear) source
term portion was implemented with an explicit time stepping scheme in order to
keep the discrete problem linear. The new wave model showed accuracy and
stability in the analytic test cases and a laboratory test case which included
source terms.

In the future, a first step to improving WAVEx would be to introduce a custom
FEniCSx assembler in order to speed up run times. One drawback of the way
WAVEx is currently implemented is that the assembly routine is relatively
slow, but this can be circumvented by developing new tools in the FEniCSx.
Other future investigations include how WAVEx scales with more complex test
cases in addition to fully implicit time stepping by extending WAVEx to solve
nonlinear problems. Additionally, implementing the DG and least squares
schemes in the full (4D) operational setting is a natural extension that we
have not done at the present time. While we have considered several FEMs,
there is a plethora of methods of interest that show promise in handling
problems like the WAE such as Lagrange-Galerkin or discontinuous Petrov-
Galerkin (DPG) [70] methods.

##  5 Acknowledgements

Author Loveland has been supported by the CSEM Fellowship from the Oden
Institute at the University of Texas at Austin. Authors Loveland, Valseth, and
Dawson have been supported by the United States National Science Foundation -
NSF PREEVENTS Track 2 Program, under NSF Grant Number 1855047 and the
Department of Homeland Security Coastal Resilience Center research project
"Accurate and Fast Wave Modeling and Coupling with ADCIRC". Author Meixner has
been supported by the National Oceanic and Atmospheric Administration. The
authors would also like to thank the reviewers of this manuscript for their
time and thoughtful suggestions.

##  Appendix A Appendix A

This appendix illustrates how to rewrite both (20) and (21) into a product of
2-dimensional integrals so that FEniCSx can interpret it through the unified
form language (UFL) framework. We evaluate the third term in (21):

$\begin{split}(\mathbf{c}\phi_{i}\psi_{j}\cdot\mathbf{n},\gamma_{k}\beta_{l})_{\partial\Omega}=(\mathbf{c}\phi_{i}\psi_{j}\cdot\mathbf{n},\gamma_{k}\beta_{l})_{\partial\Omega_{1}\times\Omega_{2}\cup\Omega_{1}\times\partial\Omega_{2}}=\\\
\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\partial\Omega_{1}}\phi_{i}\gamma_{k}\mathbf{c_{1}}\cdot\mathbf{n}_{1}dxdy+\int_{\partial\Omega_{2}}\psi_{j}\beta_{l}\mathbf{n}_{2}\cdot\int_{\Omega_{1}}\mathbf{c_{2}}\phi_{i}\beta_{k}dxdy,\end{split}$
(48)

and for the continuous Galerkin components of the form in (20) we have:

$(\frac{\partial N}{\partial t}+\nabla\cdot(\mathbf{c}N)-S,v)_{\Omega}.$ (49)

Replacing the functions with the finite element basis and omitting the
transient and source term we get:

$\begin{split}(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}),\gamma_{k}\beta_{l})_{\Omega}=(\mathbf{c}\cdot\nabla(\phi_{i}\psi_{j})+(\phi_{i}\psi_{j})\nabla\cdot\mathbf{c},\gamma_{k}\beta_{l})_{\Omega}=,\\\
(\mathbf{c_{1}}\cdot\psi_{j}\nabla_{1}\phi_{i}+\mathbf{c_{2}}\cdot\phi_{i}\nabla_{2}\psi_{j}+\phi_{i}\psi_{j}\nabla_{1}\cdot\mathbf{c_{1}}+\phi_{i}\psi_{j}\nabla_{2}\cdot\mathbf{c_{2}},\gamma_{k}\beta_{l})_{\Omega},\end{split}$
(50)

which decomposes into the following terms:

$\begin{split}(\mathbf{c_{1}}\cdot\psi_{j}\nabla_{1}\phi_{i},\gamma_{k}\beta_{l})_{\Omega}=\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\mathbf{c_{1}}\cdot\nabla_{1}(\phi_{i})\gamma_{k}dxdy,\
\\\
(\mathbf{c_{2}}\cdot\phi_{i}\nabla_{2}\psi_{j},\gamma_{k}\beta_{l})_{\Omega}=\int_{\Omega_{2}}\beta_{l}\nabla_{2}\psi_{j}\cdot\int_{\Omega_{1}}\mathbf{c_{2}}\phi_{i}\gamma_{k}dxdy,\\\
(\phi_{i}\psi_{j}\nabla_{1}\cdot\mathbf{c_{1}},\gamma_{k}\beta_{l})_{\Omega}=\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\nabla_{1}\cdot(\mathbf{c_{1}})\phi_{i}\gamma_{k}dxdy,\\\
(\phi_{i}\psi_{j}\nabla_{2}\cdot\mathbf{c_{2}},\gamma_{k}\beta_{l})_{\Omega}=\int_{\Omega_{1}}\phi_{i}\gamma_{k}\int_{\Omega_{2}}\nabla_{2}\cdot(\mathbf{c_{2}})\psi_{j}\beta_{l}dydx.\end{split}$
(51)

Now for SUPG stabilization terms in both forms we need to add the upwind
component:

$(\frac{\partial N}{\partial
t}+\nabla\cdot(\mathbf{c}N)-S,\tau\mathbf{c}\cdot\nabla(v))_{\Omega}.$ (52)

The source term will be absorbed into the RHS and we will first focus our
attention on the streamline diffusion term:

$(\nabla\cdot(\mathbf{c}N),\tau\mathbf{c}\cdot\nabla v)_{\Omega}.$ (53)

Using our product basis expansion, we consider this term now in the discrete
setting and expand using the product rule:

$\begin{split}(\nabla\cdot(\mathbf{c}N),\tau\mathbf{c}\cdot\nabla
v)_{\Omega}\rightarrow(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}),\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=,\\\
(\nabla\cdot(\mathbf{c}\phi_{i}\psi_{j}),\tau\beta_{l}\mathbf{c}_{1}\cdot\nabla_{1}\gamma_{k}+\tau\gamma_{k}\mathbf{c}_{2}\cdot\nabla_{2}\beta_{l})_{\Omega}=,\\\
(\mathbf{c}\cdot\nabla(\phi_{i}\psi_{j})+(\phi_{i}\psi_{j})\nabla\cdot\mathbf{c},\tau\beta_{l}\mathbf{c}_{1}\cdot\nabla_{1}\gamma_{k}+\tau\gamma_{k}\mathbf{c}_{2}\cdot\nabla_{2}\beta_{l})_{\Omega}=,\\\
(\mathbf{c_{1}}\cdot\psi_{j}\nabla_{1}(\phi_{i})+\mathbf{c_{2}}\cdot\phi_{i}\nabla_{2}(\psi_{j})+(\phi_{i}\psi_{j})\nabla_{1}\cdot\mathbf{c_{1}}+(\phi_{i}\psi_{j})\nabla_{2}\cdot\mathbf{c_{2}},\tau\beta_{l}\mathbf{c}_{1}\cdot\nabla_{1}\gamma_{k}+\tau\gamma_{k}\mathbf{c}_{2}\cdot\nabla_{2}\beta_{l})_{\Omega}.\end{split}$
(54)

Since we have 4 scalar terms on the left hand part of the inner product and 2
on the right, if we expand all of these out it will lead to a sum of 8 total
terms:

$\begin{split}(\mathbf{c_{1}}\cdot\psi_{j}\nabla_{1}(\phi_{i}),\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega},\\\
(\mathbf{c_{1}}\cdot\psi_{j}\nabla_{1}(\phi_{i}),\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega,}\\\
(\mathbf{c_{2}}\cdot\phi_{i}\nabla_{2}(\psi_{j}),\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega},\\\
(\mathbf{c_{2}}\cdot\phi_{i}\nabla_{2}(\psi_{j}),\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega},\\\
((\phi_{i}\psi_{j})\nabla_{1}\cdot\mathbf{c_{1}},\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega},\\\
((\phi_{i}\psi_{j})\nabla_{1}\cdot\mathbf{c_{1}},\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega},\\\
((\phi_{i}\psi_{j})\nabla_{2}\cdot\mathbf{c_{2}},\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega},\\\
((\phi_{i}\psi_{j})\nabla_{2}\cdot\mathbf{c_{2}},\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega}.\\\
\end{split}$ (55)

Now looking term by term we have first:

$(\mathbf{c_{1}}\psi_{j}\cdot\nabla_{1}\phi_{i},\tau\beta_{l}\mathbf{c_{1}}\cdot\nabla_{1}\gamma_{k})_{\Omega}=\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\tau\mathbf{c_{1}}\cdot\nabla_{1}\phi_{i}\mathbf{c_{1}}\cdot\nabla_{1}\gamma_{k}dxdy.$
(56)

The second term we get:

$(\mathbf{c}_{1}\psi_{j}\cdot\nabla_{1}\phi_{i},\tau\gamma_{k}\mathbf{c}_{2}\cdot\nabla_{2}\beta_{l})_{\Omega}=\int_{\Omega_{2}}\psi_{j}\nabla_{2}\beta_{l}\cdot\int_{\Omega_{1}}\mathbf{c_{2}}\tau\gamma_{k}\mathbf{c_{1}}\cdot\nabla_{1}\phi_{i}dxdy.$
(57)

The third term will be:

$(\mathbf{c}_{2}\phi_{i}\cdot\nabla_{2}\psi_{j},\tau\beta_{l}\mathbf{c_{1}}\cdot\nabla_{1}\gamma_{k})_{\Omega}=\int_{\Omega_{2}}\nabla_{2}\psi_{j}\beta_{l}\cdot\int_{\Omega_{1}}\tau\mathbf{c}_{2}\phi_{i}\mathbf{c_{1}}\cdot\nabla_{1}\gamma_{k}dxdy.$
(58)

Fourth term:

$(\mathbf{c_{2}}\phi_{i}\cdot\nabla_{2}\psi_{j},\tau\gamma_{k}\mathbf{c_{2}}\cdot\nabla_{2}\beta_{l})_{\Omega}=\int_{\Omega_{1}}\phi_{i}\gamma_{k}\int_{\Omega_{2}}\tau\mathbf{c_{2}}\cdot\nabla_{2}\psi_{j}\mathbf{c_{2}}\cdot\nabla_{2}\beta_{l}dydx.$
(59)

Fifth term:

$((\phi_{i}\psi_{j})\nabla_{1}\cdot\mathbf{c_{1}},\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega}=\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\tau\phi_{i}\nabla_{1}\cdot\mathbf{c_{1}}\mathbf{c_{1}}\cdot\nabla_{1}(\gamma_{k})dxdy.$
(60)

Sixth:

$((\phi_{i}\psi_{j})\nabla_{1}\cdot\mathbf{c_{1}},\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega}=\int_{\Omega_{2}}\psi_{j}\nabla_{2}(\beta_{l})\cdot\int_{\Omega_{1}}\mathbf{c_{2}}\tau\phi_{i}\gamma_{k}\nabla_{1}\cdot\mathbf{c_{1}}dxdy.$
(61)

Seventh:

$((\phi_{i}\psi_{j})\nabla_{2}\cdot\mathbf{c_{2}},\tau\mathbf{c_{1}}\cdot\beta_{l}\nabla_{1}(\gamma_{k}))_{\Omega}=\int_{\Omega_{1}}\phi_{i}\nabla_{1}(\gamma_{k})\cdot\int_{\Omega_{2}}\mathbf{c_{1}}\tau\psi_{j}\beta_{l}\nabla_{2}\cdot\mathbf{c_{2}}dydx.$
(62)

Eighth:

$((\phi_{i}\psi_{j})\nabla_{2}\cdot\mathbf{c_{2}},\tau\mathbf{c_{2}}\cdot\gamma_{k}\nabla_{2}(\beta_{l}))_{\Omega}=\int_{\Omega_{1}}\phi_{i}\gamma_{k}\int_{\Omega_{2}}\tau\mathbf{c_{2}}\cdot\nabla_{2}(\beta_{l})\psi_{j}\nabla_{2}\cdot\mathbf{c_{2}}dydx.$
(63)

Now the last part that remains is the upwinding portion based on the choice of
time step. Choosing a simple implicit time step we get:

$(\frac{\partial N}{\partial t},\tau\mathbf{c}\cdot\nabla
v)_{\Omega}\approx(\frac{\phi_{i}^{n+1}\psi_{j}^{n+1}-\phi_{i}^{n}\psi_{j}^{n}}{\Delta
t},\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}.$ (64)

The component with the previous time step $n$ will go to RHS and the part that
is in current time step $n+1$ will be in stiffness matrix. Both
descritizations however will be the same so I will just do the $n+1$ time step
after multiplying by $\Delta t$:

$(\phi_{i}\psi_{j},\tau\mathbf{c}\cdot\nabla(\gamma_{k}\beta_{l}))_{\Omega}=(\phi_{i}\psi_{j},\tau(\mathbf{c}_{x}\cdot\beta_{l}\nabla_{x}(\gamma_{k})+\mathbf{c}_{y}\cdot\gamma_{k}\nabla_{y}(\beta_{l})))_{\Omega}.$
(65)

Focusing on left term we can evaluate using this decomposition:

$(\phi_{i}\psi_{j},\tau\mathbf{c}_{x}\cdot\beta_{l}\nabla_{x}(\gamma_{k}))_{\Omega}=\int_{\Omega_{2}}\psi_{j}\beta_{l}\int_{\Omega_{1}}\tau\phi_{i}\mathbf{c_{x}}\cdot\nabla_{x}\gamma_{k}dxdy.$
(66)

Then on the right we get:

$(\phi_{i}\psi_{j},\tau\mathbf{c}_{y}\cdot\gamma_{k}\nabla_{y}(\beta_{l}))_{\Omega}=\int_{\Omega_{2}}\psi_{j}\nabla_{y}\beta_{l}\cdot\int_{\Omega_{1}}\tau\mathbf{c}_{y}\phi_{i}\gamma_{k}dxdy.$
(67)

## References

  * [1] L. H. Holthuijsen, Waves in Oceanic and Coastal Waters, Cambridge University Press, 2007.  doi:10.1017/CBO9780511618536.
  * [2] G. J. Komen, L. Cavaleri, M. Donelan, K. Hasselmann, S. Hasselmann, P. A. E. M. Janssen, Dynamics and Modelling of Ocean Waves, Cambridge University Press, 1994\.  doi:10.1017/CBO9780511628955.
  * [3] J. Monbaliu, Chapter 5 Spectral wave models in coastal areas, in: V. Lakhan (Ed.), Advances in Coastal Modeling, Vol. 67 of Elsevier Oceanography Series, Elsevier, 2003, pp. 133–158.  doi:https://doi.org/10.1016/S0422-9894(03)80122-8.
URL https://www.sciencedirect.com/science/article/pii/S0422989403801228

  * [4] M. Khandekar, Operational Analysis and Prediction of Ocean Wind Waves, Vol. 33, 1989\.  doi:10.1007/978-1-4613-8952-1.
  * [5] P. Janssen, Progress in ocean wave forecasting, Journal of Computational Physics 227 (2008) 3572–3594.  doi:10.1016/j.jcp.2007.04.029.
  * [6] I. R. Young, Wind generated ocean waves, Elsevier, 1999.
  * [7] T. W. Group, The WAM model—a third generation ocean wave prediction model, Journal of Physical Oceanography 18 (12) (1988) 1775–1810.
  * [8] H. L. Tolman, A Third-Generation Model for Wind Waves on Slowly Varying, Unsteady, and Inhomogeneous Depths and Currents, Journal of Physical Oceanography 21 (6) (1991) 782 – 797.  doi:10.1175/1520-0485(1991)021<0782:ATGMFW>2.0.CO;2.
URL
https://journals.ametsoc.org/view/journals/phoc/21/6/1520-0485_1991_021_0782_atgmfw_2_0_co_2.xml

  * [9] P. Janssen, The interaction of ocean waves and wind, Cambridge University Press, 2004.
  * [10] N. Booij, R. C. Ris, L. H. Holthuijsen, A third-generation wave model for coastal regions: 1. Model description and validation, Journal of Geophysical Research: Oceans 104 (C4) (1999) 7649–7666.  arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/98JC02622, doi:https://doi.org/10.1029/98JC02622.
URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/98JC02622

  * [11] A. Roland, Development of WWM II: Spectral wave modeling on unstructured meshes, Ph.D. thesis (01 2008).
  * [12] J. Meixner, Discontinuous Galerkin Methods for Spectral Wave/Circulation Modeling, Ph.D. thesis (07 2013).
  * [13] J.-H. Alves, F. Ardhuin, A. Babanin, M. Banner, K. Belibassakis, M. Benoit, M. Donelan, J. Groeneweg, T. Herbers, P. Hwang, P. Janssen, T. Janssen, I. Lavrenov, R. Magne, J. Monbaliu, M. Onorato, V. Polnikov, D. Resio, W. Rogers, L. Cavaleri, Wave modelling - the state of the art, Progress In Oceanography 75 (2007) 603–674.  doi:10.1016/j.pocean.2007.05.005.
  * [14] A. Roland, F. Ardhuin, On the developments of spectral wave models: numerics and parameterizations for the coastal ocean, Ocean Dynamics 64 (2014) 833–846.
  * [15] S. V. Patankar, Numerical heat transfer and fluid flow, CRC press, 1980.
  * [16] J. M. Smith, A. R. Sherlock, D. T. Resio, STWAVE: Steady-state spectral wave model user’s manual for STWAVE, Version 3.0, Tech. rep., Engineer Research and Development Center Vicksburg MS Coastal and Hydraulics Lab (2001).
  * [17] M. Benoit, F. Marcos, F. Becq, TOMAWAC. A prediction model for offshore and nearshore storm waves (1997).
  * [18] T.-W. Hsu, S.-H. Ou, J.-M. Liau, Hindcasting nearshore wind waves using a FEM code for SWAN, Coastal Engineering 52 (2005) 177–195.  doi:10.1016/j.coastaleng.2004.11.005.
  * [19] A. Roland, U. Zanke, T.-W. Hsu, S.-H. Ou, J.-M. Liau, Spectral wave modelling on unstructured grids with the WWM (Wind Wave Model) I: The deep water case, in: Third Chinese–German Joint Symposium on Coastal and Ocean Engineering, Tainan, Taiwan, 2006, p. 17.
  * [20] A. Roland, A. Cucco, C. Ferrarin, T.-W. Hsu, J.-M. Liau, S.-H. Ou, G. Umgiesser, U. Zanke, On the development and verification of a 2-D coupled wave-current model on unstructured meshes, Journal of Marine Systems 78 (2009) S244–S254.
  * [21] A. Roland, Application of Residual Distribution (RD) schemes to the geographical part of the Wave Action Equation, in: ECMWF workshop on ocean waves, 2012, pp. 25–27.
  * [22] R. Abgrall, Residual distribution schemes: current status and future trends, Computers & Fluids 35 (7) (2006) 641–669.
  * [23] A. Abdolali, A. Roland, A. Van Der Westhuysen, J. Meixner, A. Chawla, T. J. Hesser, J. M. Smith, M. D. Sikiric, Large-scale hurricane modeling using domain decomposition parallelization and implicit scheme implemented in WAVEWATCH III wave model, Coastal Engineering 157 (2020) 103656.
  * [24] M. Zijlema, Computation of wind-wave spectra in coastal waters with SWAN on unstructured grids, Coastal Engineering 57 (3) (2010) 267–277.
  * [25] B. Yildirim, G. E. Karniadakis, A hybrid spectral/DG method for solving the phase-averaged ocean wave equation: Algorithm and validation, Journal of Computational Physics 231 (14) (2012) 4921–4953.  doi:https://doi.org/10.1016/j.jcp.2012.04.013.
URL https://www.sciencedirect.com/science/article/pii/S0021999112001854

  * [26] J. Meixner, J. C. Dietrich, C. Dawson, M. Zijlema, L. H. Holthuijsen, A Discontinuous Galerkin Coupled Wave Propagation/Circulation Model, Journal of Scientific Computing 59 (2) (2014) 334–370.  doi:10.1007/s10915-013-9761-5.
URL https://doi.org/10.1007/s10915-013-9761-5

  * [27] J. Qi, C. Chen, R. C. Beardsley, W. Perrie, G. W. Cowles, Z. Lai, An unstructured-grid finite-volume surface wave model (FVCOM-SWAVE): Implementation, validations and applications, Ocean Modelling 28 (1) (2009) 153–166, the Sixth International Workshop on Unstructured Mesh Numerical Modelling of Coastal, Shelf and Ocean Flows.  doi:https://doi.org/10.1016/j.ocemod.2009.01.007.
URL https://www.sciencedirect.com/science/article/pii/S1463500309000067

  * [28] A. T. McRae, G.-T. Bercea, L. Mitchell, D. A. Ham, C. J. Cotter, Automated generation and symbolic manipulation of tensor product finite elements, SIAM Journal on Scientific Computing 38 (5) (2016) S25–S47.
  * [29] M. Habera, J. S. Hale, C. Richardson, J. Ring, M. Rognes, N. Sime, G. N. Wells, FEniCSX: A sustainable future for the FEniCS Project (2 2020).  doi:10.6084/m9.figshare.11866101.v1.
URL
https://figshare.com/articles/presentation/FEniCSX_A_sustainable_future_for_the_FEniCS_Project/11866101

  * [30] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The finite element method: its basis and fundamentals, Elsevier, 2005.
  * [31] E. B. Becker, G. F. Carey, J. T. Oden, Finite Elements: An Introduction, Vol. 1, Prentice Hall, 1981.
  * [32] G. F. Carey, J. T. Oden, Finite Elements: A Second Course, Prentice-hall, 1983.
  * [33] J. N. Reddy, Introduction to the finite element method, McGraw-Hill Education, 2019\.
  * [34] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Vol. 159, Springer Science & Business Media, 2013.
  * [35] S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Vol. 3, Springer, 2008.
  * [36] C. Johnson, Numerical solution of partial differential equations by the finite element method, Courier Corporation, 2012.
  * [37] J. T. Oden, L. F. Demkowicz, Applied functional analysis, Chapman and Hall/CRC, 2017\.
  * [38] C. Johnson, U. Navert, J. Pitkaranta, Finite element methods for linear hyperbolic problems, Computer methods in applied mechanics and engineering 45 (1984) 285–312.
  * [39] O. Pironneau, O. Pironneau, Finite element methods for fluids, Wiley Chichester, 1989.
  * [40] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, The finite element method for fluid dynamics, Butterworth-Heinemann, 2013.
  * [41] J. Donea, L. Quartapelle, An introduction to finite element methods for transient advection problems, Computer Methods in Applied Mechanics and Engineering 95 (2) (1992) 169–203.  doi:https://doi.org/10.1016/0045-7825(92)90139-B.
URL https://www.sciencedirect.com/science/article/pii/004578259290139B

  * [42] A. N. Brooks, T. J. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 32 (1) (1982) 199–259.  doi:https://doi.org/10.1016/0045-7825(82)90071-8.
URL https://www.sciencedirect.com/science/article/pii/0045782582900718

  * [43] P. B. Bochev, M. D. Gunzburger, Finite element methods of least-squares type, SIAM review 40 (4) (1998) 789–837.
  * [44] C. Dawson, J. Proft, A priori error estimates for interior penalty versions of the local discontinuous Galerkin method applied to transport equations, Numerical Methods for Partial Differential Equations: An International Journal 17 (6) (2001) 545–564.
  * [45] B. Cockburn, G. E. Karniadakis, C.-W. Shu, The development of discontinuous Galerkin methods, Springer, 2000.
  * [46] F. Brezzi, B. Cockburn, L. D. Marini, E. Süli, Stabilization mechanisms in discontinuous Galerkin finite element methods, Computer Methods in Applied Mechanics and Engineering 195 (25-28) (2006) 3293–3310.
  * [47] F. Brezzi, L. D. Marini, E. Süli, Discontinuous Galerkin methods for first-order hyperbolic problems, Mathematical models and methods in applied sciences 14 (12) (2004) 1893–1903.
  * [48] J. Douglas, Jr., T. F. Russell, Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM Journal on Numerical Analysis 19 (5) (1982) 871–885.  arXiv:https://doi.org/10.1137/0719063, doi:10.1137/0719063.
URL https://doi.org/10.1137/0719063

  * [49] R. Bermejo, J. Carpio, L. Saavedra, New error estimates of lagrange-galerkin methods for the advection equation (2023).  doi:10.48550/ARXIV.2301.03438.
URL https://arxiv.org/abs/2301.03438

  * [50] K. Morton, A. Priestley, E. Suli, Stability of the lagrange-galerkin method with non-exact integration, ESAIM: Mathematical Modelling and Numerical Analysis 22 (4) (1988) 625–653.
  * [51] C. Dawson, C. Van Duijn, M. Wheeler, Characteristic-Galerkin methods for contaminant transport with nonequilibrium adsorption kinetics, SIAM Journal on Numerical Analysis 31 (4) (1994) 982–999.
  * [52] M. Loveland, E. Valseth, M. Lukac, C. Dawson, Extending fenics to work in higher dimensions using tensor product finite elements, Journal of Computational Science 64 (2022) 101831.
  * [53] R. Ris, L. Holthuijsen, J. Smith, N. Booij, A. van Dongeren, The onr test bed for coastal and oceanic wave models, 2003, pp. 380–391.  doi:10.1142/9789812791306_0033.
  * [54] P. LeBlond, L. Mysak, Waves in the Ocean, ISSN, Elsevier Science, 1981.
URL https://books.google.com/books?id=TYGIwzz_k8kC

  * [55] R. J. LeVeque, R. J. Leveque, Numerical methods for conservation laws, Vol. 214, Springer, 1992.
  * [56] S. Team, et al., SWAN, Scientific and Technical Documentation, SWAN Cycle III version 41.41, http://www. swan. tudelft. nl, Delft University of Technology (2022).
URL https://swanmodel.sourceforge.io/download/zip/swantech.pdf

  * [57] G. J. Komen, S. Hasselmann, K. Hasselmann, On the existence of a fully developed wind-sea spectrum, Journal of Physical Oceanography 14 (8) (1984) 1271 – 1285.  doi:10.1175/1520-0485(1984)014<1271:OTEOAF>2.0.CO;2.
URL
https://journals.ametsoc.org/view/journals/phoc/14/8/1520-0485_1984_014_1271_oteoaf_2_0_co_2.xml

  * [58] M. Loveland, Wavex: A new sepctral wind wave model using stabilized finite elements, Ph.D. thesis, The University of Texas at Austin (2023).
  * [59] I. Babuska, J. Whiteman, T. Strouboulis, Finite elements: an introduction to the method and error estimation, Oxford University Press, 2010.
  * [60] A. Ern, M. Vohralík, M. Zakerzadeh, Guaranteed and robust L2-norm a posteriori error estimates for 1D linear advection problems, ESAIM: Mathematical Modelling and Numerical Analysis 55 (2021) S447–S474.
  * [61] P. Bochev, M. Gunzburger, Least-Squares Finite Element Methods, Vol. 166, 2006.  doi:10.1007/b13382.
  * [62] E. Burman, Consistent SUPG-method for transient transport problems: Stability and convergence, Computer Methods in Applied Mechanics and Engineering 199 (17-20) (2010) 1114–1123.
  * [63] E. Burman, G. Smith, Analysis of the space semi-discretized SUPG method for transient convection–diffusion equations, Mathematical Models and Methods in Applied Sciences 21 (10) (2011) 2049–2068.
  * [64] R. Codina, On stabilized finite element methods for linear systems of convection–diffusion-reaction equations, Computer Methods in Applied Mechanics and Engineering 188 (1) (2000) 61–82.  doi:https://doi.org/10.1016/S0045-7825(00)00177-8.
URL https://www.sciencedirect.com/science/article/pii/S0045782500001778

  * [65] R. Glowinski, S. Osher, W. Yin, Splitting Methods in Communication, Imaging, Science and Engineering, 2016.  doi:10.1007/978-3-319-41589-5.
  * [66] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis 5 (3) (1968) 506–517.  arXiv:https://doi.org/10.1137/0705041, doi:10.1137/0705041.
URL https://doi.org/10.1137/0705041

  * [67] S. Gottlieb, Z. J. Grant, L. Isherwood, Strong stability preserving integrating factor runge-kutta methods.  doi:10.48550/ARXIV.1708.02595.
URL https://arxiv.org/abs/1708.02595

  * [68] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al., Petsc users manual (2019).
  * [69] H. Egger, J. Schöberl, A hybrid mixed discontinuous Galerkin finite-element method for convection–diffusion problems, IMA Journal of Numerical Analysis 30 (4) (2010) 1206–1234.
  * [70] L. Demkowicz, J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. Part I: The transport equation, Computer Methods in Applied Mechanics and Engineering 199 (23) (2010) 1558–1572.  doi:https://doi.org/10.1016/j.cma.2010.01.003.
URL https://www.sciencedirect.com/science/article/pii/S0045782510000125