SQGE.Tests.tex

The main goal of this section is twofold. First, we show that the FE
discretization of the streamfunction formulation of the SQGE \eqref{eqn:SQGE_Psi}
with the Argyris element produces accurate numerical approximations. To this
end, we benchmark our numerical results against those in the published
literature \cite{Vallis06, Cascon, Myers}. The second goal is to show that the
numerical results follow the theoretical error estimates in
\autoref{thm:EnergyNorm} and \autoref{thm:Errors}, i.e., we compare the observed
rates of convergence to the theoretical rates of convergence developed in
\autoref{sec:SQGEErrors}.

Although the pure streamfunction formulation of the steady SQGE
\eqref{eqn:SQGE_Psi} is our main concern, we also test our Argyris FE
discretization on two simplified settings: \begin{inparaenum}[(i)] \item the
\emph{linear Stommel} model; and \item \emph{linear Stommel-Munk} model.
\end{inparaenum} The reason for using these two additional numerical tests is
that they are standard test problems in the geophysical fluid dynamics
literature (see, e.g., Chapter 13 in Vallis \cite{Vallis06} as well as the
reports of Myers \emph{et al.} \cite{Myers}, and Cascon \emph{et al.}
\cite{Cascon}). This allows us to benchmark our numerical results against
those in the published literature. Since both the linear Stommel and the
linear Stommel-Munk models lack the nonlinearity present in the SQGE
\eqref{eqn:SQGE_Psi}, they represent good stepping stones to verifying our FE
discretization.

The first simplified model used in our numerical investigation is the
\emph{linear Stommel} model (equation (14.22) in \cite{Vallis06} and equation
(11) in \cite{Myers})
\begin{equation}
  \epsilon_s \Delta \psi + \frac{\partial \psi}{\partial x} = f,
  \label{eqn:Stommel}
\end{equation}
where $\epsilon_s$ is the Stommel number. The Stommel number is given by
(equation (10) in \cite{Myers})
\begin{equation*}
  \epsilon_s = \frac{\gamma}{\beta L},
\end{equation*}
where $\gamma,\, \beta, \text{ and } L$ are the bottom fricition decay, the
measure of the $\beta$-plane effect, and characteristic length scale,
respectively. We note that the linear Stommel Model \eqref{eqn:Stommel} is just
the linear Stommel-Munk model \eqref{eqn:Stommel-Munk} in which the biharmonic
term is dropped (i.e. $\epsilon_m=0$). Thus it is a good starting point, since
is less complex than both the SQGE and the linear Stommel-Munk model.

The \emph{linear Stommel-Munk} model (see p. 587 in \cite{Vallis06} and problem
2 in \cite{Cascon}) is given by
\begin{equation}
  \epsilon_s \Delta \psi - \epsilon_m \Delta^2 \psi + \frac{\partial \psi}{\partial x} = f,
  \label{eqn:Stommel-Munk}
\end{equation}
where $\epsilon_s \text{ and } \epsilon_m$ are the Stommel number and Munk
scale, respectively. The Stommel number and Munk scale are given by (equation
(10) in \cite{Myers})
\begin{equation*}
  \epsilon_m = \frac{A}{\beta L^3},
\end{equation*}
where $A$ is the eddy viscosity.  The model is supplemented with appropriate
boundary conditions, which will be described for each of the subsequent
numerical tests.

We note that the linear Stommel-Munk model \eqref{eqn:Stommel-Munk} is similar
in form to the SQGE \eqref{eqn:SQGE_Psi}. Indeed, both models contain the
biharmonic operator $\Delta^2 \psi$, the rotation term $\dfrac{\partial
\psi}{\partial x}$, and the forcing terms $f$ and $F$, respectively. The two
main differences between the two models are the following: First, the SQGE are
nonlinear, since they contain the Jacobian term $J(\psi,q)$; the Stommel-Munk
model is linear, since it doesn't contain the Jacobian term. The second
difference is that the linear Stommel-Munk model contains a Laplacian term
$\Delta \psi$, whereas the SQGE does not.

We also note that the two models use different parameters: the Reynolds number
$Re$ and the Rossby number $Ro$ in \eqref{eqn:SQGE_Psi} and the Stommel number
$\epsilon_s$ and the Munk scale $\epsilon_m$ in the linear Stommel-Munk model.
These two sets of parameters, however, are related by the following relations:
\begin{align}
  \epsilon_m &= Ro\, Re^{-1}, \label{eqn:Munkscale}\\
  \epsilon_s &= Ro \frac{\gamma\, L}{U}. \label{eqn:Stommelnumber}
\end{align}
%\textcolor{red}{If $\gamma \sim T$ then $\epsilon = Ro$. However, I am unsure this is a valid
%assumption and so I am having trouble relating the Stommel number to Reynolds and Rosby.}

The rest of the section is organized as follows: In \autoref{sse:LSM} we
present results for the linear Stommel model \eqref{eqn:Stommel}. In
\autoref{sse:SMM} we present results for the linear Stommel-Munk model
\eqref{eqn:Stommel-Munk}. Finally, in \autoref{sse:SQGE} we present results for
the nonlinear SQGE \eqref{eqn:SQGE_Psi}.

\subsection{Linear Stommel Model} \label{sse:LSM} This subsection presents the
results for the FE discretization of the linear Stommel model
\eqref{eqn:Stommel} by using the Argyris element. The computational domain is
$\Omega = [0,1]\times [0,1]$. For completeness, we present results for two
numerical tests. The first test, denoted by Test 1, corresponds to the exact
solution used by Vallis (equation (14.38) in \cite{Vallis06}), while the second
test, denoted by Test 2, corresponds to the exact solution used by Myers
\emph{et al.} (Equations 15 and 16 in \cite{Myers}).

\tbf{Test 1a:} In this test, we choose the same setting as that used by Vallis
(equation (14.38) in \cite{Vallis06}).  In particular, the forcing term and the
non-homogeneous Dirichlet boundary conditions are chosen to match those given by
the exact solution
\begin{equation}
  \psi(x,y) = (1-x-e^{-\nicefrac{x}{\epsilon_s}})\, \pi\, \sin \pi y.
  \label{eqn:MyerExact}
\end{equation}
We choose the same Stommel number as that used by Vallis, i.e.
$\epsilon_s=0.04$. The exact solution \eqref{eqn:MyerExact} considered by Vallis
satisfies $\psi \to 0$ as $x \to 0$, but does not vanish at $x=1$. In our
numerical tests, we used a standard lifting procedure to treat these
non-homogeneous boundary conditions, i.e. for a problem of the form
\begin{align*}
  L \psi&=f \text{ on } \Omega\\
  \psi &=g \text{ on } \partial \Omega,
\end{align*}
we reformulate it to be
\begin{align*}
  L\tilde{\psi} &= \tilde{f} \text{ on } \Omega \\
  \tilde{\psi} &= 0 \text{ on } \partial \Omega
\end{align*}
where $L\tilde{\psi} = Lu - LS = f - g = \tilde{f}$. The solution $\psi$ from
the original problem can then be found by $\psi(x,y) =\tilde{\psi}+S$. The
function $S(x,y)$ is assumed to have the form
\begin{equation*}
  S(x,y) = A(y) (1-x) + B(y) x
\end{equation*}
and satisfies the boundary conditions given by \eqref{eqn:MyersExact}. After
some simple algebra we see that
\begin{equation*}
  S(x,y) = -x e^{-x/\epsilon_s}\pi \sin(\pi y).
\end{equation*}
The function $\tilde{f}$ can be determined by applying the operator $L$
corresponding to the \emph{linear Stommel} problem to $u - S$.

Applying the finite element method to the \emph{linear Stommel} problem with the
new modified $\tilde{f}$, corresponding to the exact solution given by Vallis,
and homogeneous boundary conditions using Argyris Finite Elements we get a
solution that matches the solution presented by Vallis, as can be seen in
\autoref{fig:StommelVallis}. Additionally, the table of errors
\autoref{tab:StommelErrorsVallis} shows the order of convergence appears to be
approaching the expected rates for $L^2, H^1, \text{ and } H^2$ norms.
\autoref{fig:StommelVallis} presents the streamlines of the approximate solution
obtained by using the Argyris Finite Element on a mesh with $h=\frac{1}{32}$ and
$9670$ DoFs. Comparing \autoref{fig:StommelVallis} with Figure $14.5$ in
\cite{Vallis06}, we notice that our approximation is close to his. Since the
exact solution is available, we can compute the errors in various norms.
\autoref{tab:StommelErrorsVallis} presents the errors $e_0,\, e_1, \text{ and }
e_2$ (i.e., the $L^2,\, H^1, \text{ and } H^2$ errors, respectively) for various
values of the mesh sizes, $h$ (the DoFs are also included).

\begin{table}%[H]
\begin{center}
%{\footnotesize
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $70$ & $1.15\times 10^{-1}$ & $-$ & $1.81\times 10^0$ & $-$ & $8.37\times 10^1$ & $-$ \\[0.2em] % & $6.81\times 10^{-3}$ & $-$ \\
  $\nicefrac{1}{4}$ & $206$ & $1.02\times 10^{-2}$ & $3.50$ & $3.12\times 10^{-1}$ & $2.54$ & $2.55\times 10^1$ & $1.72$ \\[0.2em] % & $7.801\times 10^{-3}$ & $-1.961\times 10^{-1}$ \\
  $\nicefrac{1}{8}$ & $694$ & $4.46\times 10^{-4}$ & $4.51$ & $2.59\times 10^{-2}$ & $3.59$ & $3.90\times 10^0$ & $2.71$ \\[0.2em] % & $0.4.681\times 10^{-2}$ & $4.059$ \\
  $\nicefrac{1}{16}$ & $2534$ & $1.09\times 10^{-5}$ & $5.36$ & $1.22\times 10^{-3}$ & $4.41$ & $3.49\times 10^{-1}$ & $3.48$ \\[0.2em] % & $1.732\times 19^{-5}$ & $4.756$ \\
  $\nicefrac{1}{32}$ & $9670$ & $1.97\times 10^{-7}$ & $5.79$ & $4.35\times 10^{-5}$ & $4.80$ & $2.34\times 10^{-2}$ & $3.90$ \\[0.2em] % & $5.396\times 19^{-7}$ & $5.005$ \\
  \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the linear Stommel Model
  \eqref{eqn:Stommel}, Test 1 \cite{Vallis06}.}
\label{tab:StommelErrorsVallis}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{figures/StommelAConvergence.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 1a \cite{Vallis06}:
      Observed errors versus expected rate of convergence.}
    \label{fig:StommelErrorsVallis}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/linearStommelVallis.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 1a \cite{Vallis06}: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $9670$ DoFs.}
    \label{fig:StommelVallis}
  \end{center}
\end{figure}
We note that the errors in \autoref{tab:StommelErrorsVallis} follow the
theoretical rates of convergence predicted by the estimates \eqref{eqn:H2Error}
- \eqref{eqn:L2Error} in \autoref{thm:Errors}. The orders of convergence in
\autoref{tab:StommelErrorsVallis} are close to the theoretical ones for the fine
meshes, but not as close for the coarse meshes. We think that the inaccuracies
on the coarse meshes are due to their inability to capture the thin boundary
layer on the left-hand side (i.e., at $x=0$). The finer the mesh gets, the
better this boundary layer is captured and the better the numerical accuracy
becomes.

\tbf{Test 1b:}
In the Second part of Test 1, we verify the hypothesis above, that is, whether
the degrading accuracy of the approximation is indeed due to the thin western
boundary layer. To this end, we change the Stommel number in Test 1a to be
$\epsilon_s=1$, which will result in a much thicker western boundary layer. We
then run the same numerical test as before, but with the new Stommel number. As
can be seen in \autoref{tab:StommelErrorsVallise1}, the rates of convergence are
the expected theoretical orders of convergence. This shows that the reason for
the inaccuracies in \autoref{tab:StommelErrorsVallis} were indeed due to the
thin western boundary layer.

\begin{table}%[H]
\begin{center}
%{\scriptsize
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $70$ & $1.69\times 10^{-5}$ & $-$ & $3.43\times 10^{-4}$ & $-$ & $8.72\times 10^{-3}$ & $-$ \\[0.2em] % & $4.306\times 10^{-6}$ & $-$ \\
  $\nicefrac{1}{4}$ & $206$ & $3.72\times 10^{-7}$ & $5.50$ & $1.34\times 10^{-5}$ & $4.68$ & $5.62\times 10^{-4}$ & $3.96$ \\[0.2em] % & $5.542\times 10^{-7}$ & $2.958$ \\
  $\nicefrac{1}{8}$ & $694$ & $4.89\times 10^{-9}$ & $6.25$ & $3.76\times 10^{-7}$ & $5.16$ & $3.25\times 10^{-5}$ & $4.11$ \\[0.2em] % & $9.043\times 10^{-9}$ & $5.937$ \\
  $\nicefrac{1}{16}$ & $2534$ & $7.08\times 10^{-11}$ & $6.11$ & $1.12\times 10^{-8}$ & $5.07$ & $1.96\times 10^{-6}$ & $4.05$ \\[0.2em] % & $1.355\times 10^{-10}$ & $6.06$ \\
  $\nicefrac{1}{32}$ & $9670$ & $1.08\times 10^{-12}$ & $6.04$ & $3.44\times 10^{-10}$ & $5.02$ & $1.21\times 10^{-7}$ & $4.02$ \\[0.2em] % & $2.169\times 10^{-12}$ & $5.965$ \\
  \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the linear Stommel Model
  \eqref{eqn:Stommel}, Test 1b \cite{Vallis06}.}
\label{tab:StommelErrorsVallise1}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{figures/StommelBConvergence.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 1b \cite{Vallis06}:
      Observed errors versus expected rate of convergence.}
    \label{fig:StommelErrorsVallis1}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/StommelVallise1.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 1b \cite{Vallis06}: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $9670$ DoFs with $\epsilon_s=1$.}
    \label{fig:StommelVallise1}
  \end{center}
\end{figure}

\tbf{Test 2:}
For our second test we use the exact solution given by Myers (Equations 15 and
16 in \cite{Myers}), i.e.
{\footnotesize
\begin{equation}
  \psi(x,y) =\frac{\sin(\pi y)}{\pi(1+4\pi^2\epsilon_s^2)}\left\{2\pi\epsilon_s\sin(\pi x)+cos(\pi x)+\frac{1}{e^{R_1}-e^{R_2}}\left[(1+e^{R_2})e^{R_1x}-(1+e^{R_1})e^{R_2x}\right]\right\},
  \label{eqn:MyersExact}
\end{equation}
}
where $R_1\text{ and } R_2$ are the positive and negative roots, respectively,
of
\begin{equation*}
  R = \frac{-1\pm\sqrt{1+4\pi^2 \epsilon_s^2}}{2\epsilon_s}.
\end{equation*}
The forcing term and the homogeneous Dirichlet boundary conditions are chosen to
match those given by the exact solution \eqref{eqn:MyersExact}. We choose the
same Stommel number as that used by Myers, i.e. $\epsilon_s=0.05$.

\autoref{fig:StommelMyers} presents the streamlines of the approximate solution
obtained by using the Argyris Finite Element on a mesh with $h=\frac{1}{32}$ and
$9670$ DoFs. Comparing \autoref{fig:StommelMyers} with Figure $2$ in
\cite{Myers}, we notice that our approximation is close to that in \cite{Myers}.
Since the exact solution is available, we can compute the errors in various
norms. \autoref{tab:StommelErrorsMyers} presents the errors $e_0,\, e_1, \text{
and } e_2$ (i.e., the $L^2,\, H^1, \text{ and } H^2$ errors, respectively) for
various values of the mesh sizes, $h$.

\begin{table}%[H]
\begin{center}
%{\footnotesize
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $70$ & $5.65\times 10^{-3}$ & $-$ & $1.45\times 10^{-1}$ & $-$ & $6.60\times 10^0$ & $-$ \\[0.2em] % & $8.815\times 10^{-3}$ & $-$ \\
  $\nicefrac{1}{4}$ & $206$ & $4.28\times 10^{-4}$ & $3.72$ & $2.08\times 10^{-2}$ & $2.80$ & $1.63\times 10^0$ & $2.02$ \\[0.2em] % & $0.3.139\times 10^{-2}$ & $1.49$ \\
  $\nicefrac{1}{8}$ & $694$ & $1.46\times 10^{-5}$ & $4.87$ & $1.41\times 10^{-3}$ & $3.89$ & $2.07\times 10^{-1}$ & $2.98$ \\[0.2em] % & $2.15\times 10^{-5}$ & $3.868$ \\
  $\nicefrac{1}{16}$ & $2534$ & $2.95\times 10^{-7}$ & $5.63$ & $5.83\times 10^{-5}$ & $4.59$ & $1.65\times 10^{-2}$ & $3.65$ \\[0.2em] % & $7.097\times 10^{-7}$ & $4.921$ \\
  $\nicefrac{1}{32}$ & $9670$ & $4.97\times 10^{-9}$ & $5.89$ & $2.00\times 10^{-6}$ & $4.87$ & $1.07\times 10^{-3}$ & $3.95$ \\[0.2em] % & $2.054\times 10^{-8}$ & $5.111$ \\
  \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the linear Stommel Model
  \eqref{eqn:Stommel}, Test 2 \cite{Myers}.}
\label{tab:StommelErrorsMyers}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{Figures/StommelMyersConvergence.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 2 \cite{Myers}:
      Observed errors versus expected rate of convergence.}
    \label{fig:StommelErrorsMeyers}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/linearStommelMyers.pdf}
    \caption{Linear Stommel model \eqref{eqn:Stommel}, Test 2 \cite{Myers}: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $9670$ DoFs.}
    \label{fig:StommelMyers}
  \end{center}
\end{figure}

We note that the errors in \autoref{tab:StommelErrorsMyers} follow the
theoretical rates of convergence predicted by the estimates \eqref{eqn:H2Error}
- \eqref{eqn:L2Error} in \autoref{thm:Errors}. Again, we see that the orders of
convergence in \autoref{tab:StommelErrorsMyers} are close to the theoretical
ones for the fine meshes, but not as close for the coarse meshes. We attribute
this to the inaccuracies at the thin boundary layer on the left-hand side (i.e.,
at $x=0$). The finer the mesh gets, the better this boundary layer is captured
and the better the numerical accuracy becomes.

\subsection{Linear Stommel-Munk Model}\label{sse:SMM}
This subsection presents results for the FE discretization of the linear
Stommel-Munk model \eqref{eqn:Stommel-Munk} by using the Argyris element. Our
computational setting is the same as that used by Cascon \emph{et al.}
\cite{Cascon}: The computational domain is $\Omega = [0,3]\times[0,1]$, the Munk
scale is $\epsilon_m=6\times 10^{-5}$, the Stommel number is $\epsilon_s=0.05$,
and the boundary conditions are
\begin{equation} \label{eqn:SMProb}
  \psi = \frac{\partial \psi}{\partial \mathbf{n}}=0 \quad \text{ on } \partial\Omega
\end{equation}
For completeness, we present results for two numerical tests, denoted by Test 3
and Test 4, both corresponding to Test 1 and Test 2 in \cite{Cascon},
respectively.

\tbf{Test 3:}
For our third test we use the exact solution given by Test 1 in \cite{Cascon},
i.e.
\begin{equation}
  \psi(x,y) = \sin^2 \frac{\pi x}{3} \sin^2 \pi y.
  \label{eqn:CasconExact1}
\end{equation}
The forcing term is chosen to match that given by the exact solution
\eqref{eqn:CasconExact1}.

For this third test we take $F$ corresponding to applying the linear operator
$L$ associated with the \emph{linear Stommel-Munk} model to the exact solution
\eqref{eqn:CasconExact1}.

\autoref{fig:StommelMunkSin} presents the streamlines of the approximate
solution obtained by using the Argyris Finite Element on a mesh with
$h=\frac{1}{32}$ and $28550$ DoFs. Comparing \autoref{fig:StommelMunkSin} with
Figure $7$ in \cite{Cascon}, we notice that our approximation is close to that in
\cite{Myers}. Since the exact solution is available, we can compute the errors
in various norms. \autoref{tab:SMsinErrors} presents the errors $e_0,\, e_1,
\text{ and } e_2$ (i.e., the $L^2,\, H^1, \text{ and } H^2$ errors,
respectively) for various values of the mesh sizes, $h$.

\begin{table}%[H]
\begin{center}
%{\scriptsize
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $170$ & $2.99\times 10^{-3}$ & $-$ & $4.08\times 10^{-2}$ & $-$ & $7.62\times 10^{-1}$ & $-$ \\[0.2em] % & $3.423\times 10^{-3}$ & $-$ \\
  $\nicefrac{1}{4}$ & $550$ & $3.22\times 10^{-5}$ & $6.54$ & $1.03\times 10^{-3}$ & $5.31$ & $4.08\times 10^{-2}$ & $4.23$ \\[0.2em] % & $1.885\times 10^{-5}$ & $7.505$ \\
  $\nicefrac{1}{8}$ & $1958$ & $3.44\times 10^{-7}$ & $6.55$ & $2.49\times 10^{-5}$ & $5.37$ & $2.25\times 10^{-3}$ & $4.18$ \\[0.2em] % & $2.371\times 10^{-7}$ & $6.313$ \\
  $\nicefrac{1}{16}$ & $7366$ & $4.57\times 10^{-9}$ & $6.23$ & $7.03\times 10^{-7}$ & $5.15$ & $1.34\times 10^{-4}$ & $4.07$ \\[0.2em] % & $4.296\times 10^{-9}$ & $5.786$ \\
  $\nicefrac{1}{32}$ & $28550$ & $6.70\times 10^{-11}$ & $6.09$ & $2.11\times 10^{-8}$ & $5.06$ & $8.26\times 10^{-6}$ & $4.02$ \\[0.2em] % & $6.86\times 10^{-11}$ & $5.969$ \\
 \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the linear Stommel-Munk Model
  \eqref{eqn:SMProb}, Test 3 \cite{Cascon}.}
\label{tab:SMsinErrors}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SMsinConvergence.pdf}
    \caption{Linear Stommel-Munk model \eqref{eqn:SMProb}, Test 3 \cite{Cascon}:
      Observed errors versus expected rate of convergence.}
    \label{fig:SMSinErrors}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/StommelMunk1.pdf}
    \caption{Linear Stommel-Munk model \eqref{eqn:SMProb}, Test 3 \cite{Cascon}:
    Streamlines of the approximation, $\psi^h$, on a mesh size,
    $h=\frac{1}{32}$, and $28550$ DoFs.}
    \label{fig:StommelMunkSin}
  \end{center}
\end{figure}

We note that the errors in \autoref{tab:SMsinErrors} follow the theoretical
rates of convergence predicted by the estimates \eqref{eqn:H2Error} -
\eqref{eqn:L2Error} in \autoref{thm:Errors}.  This time, we see that the orders
of convergence in \autoref{tab:SMsinErrors} are close to the theoretical ones
for the fine meshes, but are higher than expected for the for the coarse meshes.
We attribute this to the fact that the exact solution \eqref{eqn:CasconExact1}
does not display any boundary layers that could be challenging to capture by the
Argyris element on a coarse mesh.

\tbf{Test 4:}
For our fourth test we use the exact solution given by Test 2 in \cite{Cascon},
i.e.
{\small
\begin{equation}
  \psi(x,y) = \left[\left(1 - \frac{x}{3}\right)\left(1-e^{-20x}\right) \sin \pi y\right]^2.
  \label{eqn:CasconExact2}
\end{equation}
}
Again we take the forcing term $F$ corresponding the exact solution
\eqref{eqn:CasconExact2}.

\autoref{fig:SMe} presents the streamlines of the approximate solution obtained
by using the Argyris Finite Element on a mesh with $h=\frac{1}{32}$ and $28550$
DoFs. Comparing \autoref{fig:SMe} with Figure $10$ in \cite{Myers}, we notice
that our approximation is close to \cite{Myers}. Since the exact solution is
available, we can compute the errors in various norms. \autoref{tab:SMeErrors}
presents the errors $e_0,\, e_1, \text{ and } e_2$ (i.e., the $L^2,\, H^1,
\text{ and } H^2$ errors, respectively) for various values of the mesh sizes,
$h$.

\begin{table}%[H]
\begin{center}
%{\small
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $170$ & $6.04\times 10^{-2}$ & $-$ & $1.16\times 10^0$ & $-$ & $3.90\times 10^1$ & $-$ \\[0.2em] % & $2.907\times 10^{-2}$ & $-$ \\
  $\nicefrac{1}{4}$ & $550$ & $1.13\times 10^{-2}$ & $2.41$ & $4.00\times 10^{-1}$ & $1.54$ & $2.14\times 10^1$ & $0.866$ \\[0.2em] % & $5.678\times 10^{-3}$ & $2.356$ \\
  $\nicefrac{1}{8}$ & $1958$ & $8.40\times 10^{-4}$ & $3.75$ & $5.91\times 10^{-2}$ & $2.76$ & $5.66\times 10^0$ & $1.92$ \\[0.2em] % & $0.5.973\times 10^{-2}$ & $3.249$ \\
  $\nicefrac{1}{16}$ & $7366$ & $2.82\times 10^{-5}$ & $4.90$ & $4.01\times 10^{-3}$ & $3.88$ & $7.38\times 10^{-1}$ & $2.94$ \\[0.2em] % & $2.979\times 10^{-5}$ & $4.326$ \\
  $\nicefrac{1}{32}$ & $28550$ & $5.59\times 10^{-7}$ & $5.66$ & $1.61\times 10^{-4}$ & $4.64$ & $5.97\times 10^{-2}$ & $3.63$ \\[0.2em] % & $8.632\times 10^{-7}$ & $5.109$ \\
 \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the linear Stommel-Munk Model
  \eqref{eqn:SMProb}, Test 4 \cite{Cascon}.}
\label{tab:SMeErrors}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SMeConvergence.pdf}
    \caption{Linear Stommel-Munk model \eqref{eqn:SMProb}, Test 4 \cite{Cascon}:
      Observed errors versus expected rate of convergence.}
    \label{fig:SMeErrors}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/StommelMunk2.pdf}
    \caption{Linear Stommel-Munk model \eqref{eqn:SMProb}, Test 4 \cite{Cascon}: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $28550$ DoFs.}
    \label{fig:SMe}
  \end{center}
\end{figure}

We note that the errors in \autoref{tab:SMeErrors} follow the theoretical rates
of convergence predicted by the estimates \eqref{eqn:H2Error} -
\eqref{eqn:L2Error} in \autoref{thm:Errors}. Again, we see that the orders of
convergence in \autoref{tab:SMeErrors} are close to the theoretical ones for the
fine meshes, but not as close for the coarse meshes. As stated previously, we
attribute this to the inaccuracies at the thin boundary layer on the left-hand
side (i.e., at $x=0$). The finer the mesh gets, the better this boundary layer
is captured and the better the numerical accuracy becomes.

\subsection{SQGE}\label{sse:SQGE}
This subsection presents results for the FE discretization of the streamfunction
formulation of the SQGE \eqref{eqn:SQGE_Psi} by using the Argyris element. Our
computational domain is $\Omega=[0,3]\times[0,1]$, the Reynolds number is
$Re=1.667$, and the Rossby number is $Ro=10^{-4}$.  For completeness, we present
results for two numerical tests, denoted by Test 5 and Test 6, both
corresponding to the exact solutions given in Test 1 and Test 2 of
\cite{Cascon}, respectively.

\tbf{Test 5:}
In this test, we take the same exact solution presented in \emph{Test 3}, i.e.
\begin{equation}
  \psi(x,y) = \sin^2 \frac{\pi x}{3} \, \sin^2 \pi y.
  \label{eqn:StreamfunctionExact1}
\end{equation}
Again, the forcing term $F$ and homogeneous boundary conditions, $\psi = \frac{\partial
\psi}{\partial \mathbf{n}} = 0$, correspond to the exact solution \eqref{eqn:StreamfunctionExact1}.

\autoref{fig:SQGEsin} presents the streamlines of the approximate solution
obtained by using the Argyris Finite Element on a mesh with $h=\frac{1}{32}$ and
$36150$ DoFs. We note that the streamlines look as we expect and are similar to
those given by Figure $7$ in \cite{Myers}, which uses the same exact solution.
Since the exact solution is available, we can compute the errors in various
norms. \autoref{tab:SQGEsinErrors} presents the errors $e_0,\, e_1, \text{ and }
e_2$ (i.e., the $L^2,\, H^1, \text{ and } H^2$ errors, respectively) for various
values of the mesh sizes, $h$.

\begin{table}%[H]
\begin{center}
%{\scriptsize
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $138$ & $4.10\times 10^{-3}$ & $-$ & $4.91\times 10^{-2}$ & $-$ & $8.61\times 10^{-1}$ & $-$\\
  $\nicefrac{1}{4}$ & $536$ & $2.23\times 10^{-5}$ & $7.52$ & $7.66\times 10^{-4}$ & $6.00$ & $2.96\times 10^{-2}$ & $4.86$\\
  $\nicefrac{1}{8}$ & $2349$ & $2.35\times 10^{-7}$ & $6.57$ & $1.60\times 10^{-5}$ & $5.58$ & $1.31\times 10^{-3}$ & $4.50$\\
  $\nicefrac{1}{16}$ & $9152$ & $3.03\times 10^{-9}$ & $6.28$ & $4.33\times 10^{-7}$ & $5.21$ & $7.53\times 10^{-5}$ & $4.12$\\
  $\nicefrac{1}{32}$ & $36150$ & $4.64\times 10^{-11}$ & $6.03$ & $1.37\times 10^{-8}$ & $4.99$ & $4.84\times 10^{-6}$ & $3.96$\\
  $\nicefrac{1}{64}$ & $146090$ & $9.86\times 10^{-13}$ & $5.56$ & $4.10\times
  10^{-10}$ & $5.06$ & $2.92\times 10^{-7}$ & $4.05$ \\[0.2em]
 \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the SQGE \eqref{eqn:SQGE_Psi}, Test 5.}
\label{tab:SQGEsinErrors}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SQGEsin2Convergence.pdf}
    \caption{SQGE \eqref{eqn:SQGE_Psi}, Test 5: Observed errors versus expected
      rate of convergence.}
    \label{fig:SQGEsinErrors}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SQGEsin.pdf}
    \caption{SQGE \eqref{eqn:SQGE_Psi}, Test 5: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $36150$ DoFs.}
    \label{fig:SQGEsin}
  \end{center}
\end{figure}

We note that the errors in \autoref{tab:SQGEsinErrors} follow the theoretical
rates of convergence predicted by the estimates \eqref{eqn:H2Error} -
\eqref{eqn:L2Error} in \autoref{thm:Errors}. Again, since the exact solution
\eqref{eqn:StreamfunctionExact1} does not display any boundary layers, we see
that the orders of convergence in \autoref{tab:SQGEsinErrors} are close to the
theoretical ones for the fine meshes, but are higher than expected for the for
the coarse meshes. We also note the drop in the rate of convergence in the
$L^2$-norm for $\nicefrac{1}{64}$ is due to reaching the machine epsilon and
thus any rate of convergence for mesh size smaller than $\nicefrac{1}{64}$ is
meaningless.

\tbf{Test 6:}
In this test, we take the same exact solution as used in \emph{Test 4}, i.e.
\begin{equation}
  \psi(x,y) = \left[\left(1 - \frac{x}{3}\right)\left(1-e^{-20x}\right) \sin \pi y\right]^2.
  \label{eqn:StreamfunctionExact2}
\end{equation}
The forcing term $F$ and the homogeneous boundary conditions correspond to the exact solution
\eqref{eqn:StreamfunctionExact2}.

\autoref{fig:SQGEe} presents the streamlines of the approximate solution
obtained by using the Argyris Finite Element on a mesh with $h=\frac{1}{32}$ and
$36150$ DoFs. We note that the streamlines look as we expect and are similar to
those given by Figure $10$ in \cite{Myers}, which uses the same exact solution.
Since the exact solution is available, we can compute the errors in various
norms. \autoref{tab:SQGEeErrors} presents the errors $e_0,\, e_1, \text{ and }
e_2$ (i.e., the $L^2,\, H^1, \text{ and } H^2$ errors, respectively) for various
values of the mesh sizes, $h$.

\begin{table}%[H]
\begin{center}
%{\small
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em] % & $e_{\infty}$ & $L_{\infty}$ order \\
  \hline
  $\nicefrac{1}{2}$ & $138$ & $4.36\times 10^{-1}$ & $-$ & $2.22\times 10^0$ & $-$ & $4.58\times 10^1$ & $-$\\
  $\nicefrac{1}{4}$ & $536$ & $7.21\times 10^{-3}$ & $5.92$ & $2.86\times 10^{-1}$ & $2.96$ & $1.65\times 10^1$ & $1.48$\\
  $\nicefrac{1}{8}$ & $2349$ & $1.92\times 10^{-3}$ & $1.91$ & $1.07\times 10^{-1}$ & $1.42$ & $8.06\times 10^{0}$ & $1.03$\\
  $\nicefrac{1}{16}$ & $9152$ & $2.89\times 10^{-5}$ & $6.05$ & $3.94\times 10^{-3}$ & $4.76$ & $6.84\times 10^{-1}$ & $3.56$\\
  $\nicefrac{1}{32}$ & $36150$ & $6.63\times 10^{-7}$ & $5.45$ & $1.75\times 10^{-4}$ & $4.50$ & $5.77\times 10^{-2}$ & $3.57$\\
  $\nicefrac{1}{64}$ & $146090$ & $1.37\times 10^{-8}$ & $5.60$ & $6.91\times 10^{-6}$ & $4.66$ & $4.36\times 10^{-3}$ & $3.73$\\
  $\nicefrac{1}{128}$ & $585048$ & $1.47\times 10^{-10}$ & $6.55$ & $1.63\times 10^{-7}$ & $5.40$ & $2.28\times 10^{-4}$ & $4.26$ \\
 \hline
\end{tabular}
%}
\end{center}
\caption{Errors and rate of convergence for the SQGE \eqref{eqn:SQGE_Psi}, Test 6.}
\label{tab:SQGEeErrors}
\end{table}

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SQGEeConvergence.pdf}
    \caption{SQGE \eqref{eqn:SQGE_Psi}, Test 6: Observed errors versus expected
      rate of convergence.}
    \label{fig:SQGEeErrors}
  \end{center}
\end{figure}

\begin{figure}%[H]
  \begin{center}
    \includegraphics[scale=0.5]{Figures/SQGEe.pdf}
    \caption{SQGE \eqref{eqn:SQGE_Psi}, Test 6: Streamlines of the approximation,
    $\psi^h$, $h=\frac{1}{32}$, and $36150$ DoFs.}
    \label{fig:SQGEe}
  \end{center}
\end{figure}

We note that the errors in \autoref{tab:SQGEeErrors} follow the theoretical
rates of convergence predicted by the estimates \eqref{eqn:H2Error} -
\eqref{eqn:L2Error} in \autoref{thm:Errors}. Again for an exact solution which
has a western boundary layer, we see that the orders of convergence in
\autoref{tab:SQGEeErrors} are close to the theoretical ones for the fine meshes,
but not as close for the coarse meshes. We attribute this to the inaccuracies at
the thin boundary layer on the left-hand side (i.e., at $x=0$). The finer the
mesh gets, the better this boundary layer is captured and the better the
numerical accuracy becomes.

\subsection*{Mediterranean Sea}
We have created a FE mesh of the Mediterranean Sea using GMSH \cite{GMSH}.  The
coastline data was obtained from GSHHS \cite{GSHHS}. Major islands such as
Corsica, Sardinia, and Sicily were connected to the nearest land mass in order
to ensure a unique streamfunction (see the discussions in
\cite{Gunzburger89,van-Gijzen1998}).  Additionally, the Atlantic Ocean was
closed off from the Mediterranean Sea at the Straits of Gibraltar, from the Red
Sea at the Suez Canal, and the Sea of Marmara at the Dardanelles Strait, while
the Gulf of Corinth was treated as land. The resultant FE mesh is displayed in
\autoref{fig:MedMesh}.

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.4]{Figures/MediterraneanMesh.png}
    \caption{Mesh of the Mediterranean Sea created using GMSH \cite{GMSH}. The
    mesh size corresponds to $h=\frac{1}{320}$ with $DoF=240,342$.}
    \label{fig:MedMesh}
  \end{center}
\end{figure}

The following experiment used a forcing function given by $F =
\sin\left(\frac{\pi}{4} y\right)$ which is the same forcing function given by Bryan
\cite{Bryan1963}. Note: Bryan, in fact, used $\sin(\frac{\pi}{2} y)$, but on a
domain with a vertical extent which was of length two. In our case the vertical
extent of the Northern Hemisphere is of length one. Thus, the appropriate forcing
function becomes $\sin(\frac{\pi}{4} y)$. We choose similar parameters as those
used in \cite{delSastre04} and are summarized in \autoref{tab:OceanParameters}.

\begin{table}
  \begin{center}
  \begin{tabular}{|l|l|}
    \hline
    $A$ & $2000\,m^2s^{-1}$\\
    \hline
    $\theta_0$ & $40^\circ$ \\
    \hline
    $\omega$ & $7.2526\times 10^{-5}\,s^{-1}$ \\
    \hline
    $H$ & $1000\,m$ \\
    \hline
    $L$ & $1000\,km$ \\
    \hline
    $r_e$ & $6378.1\,km$ \\
    \hline
    $\rho$ & $1024\, \nicefrac{kg}{m^3}$ \\
    \hline
  \end{tabular}
  \end{center}
  \caption{Table of parameter values used for the simulations of the
    Mediterranean Sea \cite{delSastre04}, where $A,\theta_0,\omega,H,L,r_e,\rho$ are the
    eddy viscosity, reference angle for the $\beta$-plane approximation, angular
    velocity of the Earth, domain height, domain length, radius of the Earth, and
    density of seawater, respectively.}
  \label{tab:OceanParameters}
\end{table}

Taking the derivative of relation \eqref{eqn:beta_plane} with respect to $y$
gives the following relation for $\beta$ to $f$ and using the equation (2.80)
from \cite{Vallis06} gives
\begin{equation}
  \beta = \frac{2\omega}{r_e}\cos \theta_0.
  \label{eqn:Beta}
\end{equation}
From the parameters given in \autoref{tab:OceanParameters} we see that
\begin{equation*}
  \beta \approx 1.742\times 10^{-11}\, m^{-1}\,s^{-1}.
\end{equation*}
Using this approximation for $\beta$ and \eqref{eqn:velocity_scale} with
$\tau_0 = 0.6\, dyne\, cm^{-2}$ \cite{Hellerman} gives the following
approximation for the characteristic velocity:
\begin{equation*}
%  \begin{split}
%    U &= \frac{3.1415 \, 0.6 dyne\, cm^{-2}}{1027 \nicefrac{kg}{m^3} \,
%      1000\, m \, 1.742 \times 10^{-11}\,m^{-1}\,s^{-1} \, 1000\, km} \\
    U \approx 1.054\times 10^{-2} \nicefrac{m}{s}.
%  \end{split}
\end{equation*}
Therefore by \eqref{eqn:rossby_number}, the Rossby number is
\begin{equation*}
%  Ro = \frac{1.054\times 10^{-2} \nicefrac{m}{s}}{1.742\times 10^{-11} m^{-1}
%    s^{-1} (1000 km)^2}
  Ro = 6.051\times 10^{-4}
\end{equation*}
and by \eqref{eqn:reynolds_number} the Reynolds number is
\begin{equation*}
  Re = 5.27.
\end{equation*}

The solution obtained by applying the SQGE to the Mediterranean with parameters
given in \autoref{tab:OceanParameters} and forcing function $F =
\sin(\frac{\pi}{2} y)$ can be seen in \autoref{fig:SQGEMed}. The shape of the
streamfunction is very much in agreement with the solution obtained by Galan del
Sastre (Figure 2.22 in \cite{Galan-del-Sastre2004}). Additionally,  we note that
the resultant large scale structures in \autoref{fig:SQGEMed} are very similar
to those in \autoref{fig:GyresMed}. It should be noted that some of the
differences between our solution and the solution obtained by Galan del Sastre's
and our solution and the observed large scale structures seen in
\autoref{fig:GyresMed} are likely a result of his inclusion of more islands and
not connected those islands to Europe or Africa. The rates of convergence are
much lower than expected, however this likely is due to the fact that the domain
is non-smooth which is likely causing the solution to be in a space which is
less regular than the $H^6$ regularity \cite{Soane2010} assumed by the error
analysis in \autoref{sec:SQGEErrors}.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|}%c|c|}
  \hline
  $h$ & $DoFs$ & $e_0$ & $L_2$ order & $e_1$ & $H^1$ order & $e_2$ & $H^2$ order \\[0.2em]
  \hline
  $\nicefrac{1}{20}$ & $1122$ & $2.08\times 10^{-6}$ & $-$ & $1.95\times 10^{-4}$ & $-$ & $4.50\times 10^{-2}$ & $-$ \\
  $\nicefrac{1}{40}$ & $4092$ & $8.00\times 10^{-7}$ & $1.38$ & $6.68\times 10^{-5}$ & $1.54$ & $2.50\times 10^{-2}$ & $0.850$ \\
  $\nicefrac{1}{80}$ & $15594$ & $2.91\times 10^{-7}$ & $1.46$ & $2.47\times 10^{-5}$ & $1.43$ & $1.49\times 10^{-2}$ & $0.741$ \\
  $\nicefrac{1}{160}$ & $60846$ & $1.04\times 10^{-7}$ & $1.49$ & $9.05\times 10^{-6}$ & $1.45$ & $8.67\times 10^{-3}$ & $0.785$ \\
  $\nicefrac{1}{320}$ & $240342$ & $3.10\times 10^{-8}$ & $1.75$ & $2.75\times 10^{-6}$ & $1.72$ & $4.35\times 10^{-3}$ & $0.994$ \\
\hline
\end{tabular}
\caption{Observed rates of convergence for SQGE applied to the Mediterranean Sea
  with forcing function $F = \sin(\frac{\pi}{2} y)$ and true solution obtained
  from fine mesh with $h=\dfrac{1}{640}$.}
  \label{tab:SQGEMedConvergence}
\end{table}

\begin{figure}
  \begin{subfigure}[b]{\textwidth}
  \begin{center}
    \includegraphics{Figures/SQGEMedRe5_27Ro6_051E-4h320sin0_25piy.pdf}
    \caption{``True'' solution of SQGE on the Mediterranean Sea with $Re = 5.27, Ro
      = 6.051\times 10^{-4}, F = \sin( \frac{\pi}{2} y)$, and $955,302$ DoFs
      corresponding to $h = \frac{1}{640}$.}
    \label{fig:SQGEMed}
  \end{center}
  \end{subfigure}
  \begin{subfigure}[b]{\textwidth}
    \begin{center}
      \includegraphics[scale=0.5]{Figures/GyresOfMediterranean.pdf}
      \caption{Large scale currents of the Mediterranean Sea \cite{Ayoub1998}.}
    \label{fig:GyresMed}
    \end{center}
  \end{subfigure}
  \caption{Comparison of large scale structures visible in the numerical results
    of the SQGE applied to the Mediterranean Sea and the observed large scale
    structures of the Mediterranean Sea.}
  \label{fig:MedCompare}
\end{figure}