%%%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a general template file for the LaTeX package SVJour3
% for Springer journals.          Springer Heidelberg 2006/03/15
%
% Copy it to a new file with a new name and use it as the basis
% for your article. Delete % signs as needed.
%
% This template includes a few options for different layouts and
% content for various journals. Please consult a previous issue of
% your journal as needed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
% your LaTeX will extract the file if required
%\begin{filecontents*}{example.eps}
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
%%EndComments
%gsave
%newpath
%  20 20 moveto
%  20 220 lineto
%  220 220 lineto
%  220 20 lineto
%closepath
%2 setlinewidth
%gsave
%  .4 setgray fill
%grestore
%stroke
%grestore
%\end{filecontents*}
%
\documentclass{svjour3}                     % onecolumn (standard format)
%\documentclass[smallextended]{svjour3}     % onecolumn (second format)
%\documentclass[twocolumn]{svjour3}         % twocolumn
%
\smartqed  % flush right qed marks, e.g. at end of proof
\usepackage{graphicx}

% \usepackage{mathptmx}      % use Times fonts if available on your TeX system
%
% insert here the call for the packages your document requires
%\usepackage{latexsym}
% etc.
%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
\usepackage[comma,authoryear]{natbib}
\usepackage{ulem}
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}
%
% please place your own definitions here and don't use \def but
% \newcommand{}{}
%Shortcut definitions
\newcommand{\HO}{$(g1,p6)$}
\newcommand{\LO}{$(g4,p1)$}
\newcommand{\jmp}[1]{[\![#1]\!]}
\newcommand{\mean}[1]{\{\!\{#1\}\!\}}
\newcommand{\nrml}{\hat{\mathbf n}}
\newcommand{\ds}{\displaystyle}
\newcommand{\sref}[1]{\S \ref{sec:#1}}
\newcommand{\eref}[1]{(\ref{eq:#1})}
\newcommand{\tref}[1]{Table [\ref{tab:#1}]}
\newcommand{\fref}[1]{Figure [\ref{fig:#1}]}
\newcommand{\der}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\norm}[1]{\left|\!\left|#1\right|\!\right|_2}

% Insert the name of "your journal" with
 \journalname{Ocean Dynamics}
\begin{document}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{High Order Schemes for 2D Unsteady Biogeochemical Ocean Models  %\thanks{Grants or other notes
%about the article that should go on the front page should be
%placed here. General acknowledgments should be placed at the end of the article.}
}
%\subtitle{Do you have a subtitle?\\ If so, write it here}

%\titlerunning{Short form of title}        % if too long for running head

\author{\vskip-8pt M. P. Ueckermann \and  P. F. J. Lermusiaux} %etc.

%\authorrunning{Short form of author list} % if too long for running head

\institute{Mattheus P. Ueckermann and  Pierre F. J. Lermusiaux \at
              Massachusetts Institute of Technology, 77 Mass. Avenue, Cambridge, MA 02139 \\
              Tel.: +1-617-253-7799\\
              \email{mpuecker@mit.edu}, \email{pierrel@mit.edu}     %  \\
%             \emph{Present address:} of F. Author  %  if needed
%           \and
%           S. Author \at
%              MIT
}
\date{\vskip-70pt Received: date / Accepted: date}
% The correct dates will be entered by the editor


\maketitle
\begin{abstract}\vskip-18pt

Accurate numerical modeling of biogeochemical ocean dynamics is essential for numerous applications, including coastal ecosystem science, environmental management and energy, and climate dynamics. Evaluating computational requirements for such often highly nonlinear and multiscale dynamics is critical. To do so, we complete comprehensive numerical analyses, comparing low to high order discretization schemes, both in time and space, employing standard and hybrid discontinuous Galerkin Finite Element Methods, on both straight and new curved elements. Our analyses and syntheses focus on nutrient-phytoplankton-zooplankton dynamics under advection and diffusion within an ocean strait or sill, in an idealized two-dimensional geometry. For the dynamics, we investigate three biological regimes, one with single stable points at all depths and two with stable limit cycles. We also examine interactions that are dominated by the biology, by the advection, or that are balanced. For these regimes and interactions, we study the sensitivity to multiple numerical parameters including: quadrature-free and quadrature-based discretizations of the source terms; order of the spatial discretizations of advection and diffusion operators; order of the temporal discretization in explicit schemes; and resolution of the spatial mesh, with and without curved elements. A first finding is that both quadrature-based and quadrature-free discretizations give accurate results in well-resolved regions, but the quadrature-based scheme has smaller errors in under-resolved regions. We show that low-order temporal discretizations allow rapidly growing numerical errors in biological fields. We find that if a spatial discretization (mesh resolution and polynomial degree) does not resolve the solution, oscillations due to discontinuities in tracer fields can be locally significant for both low and high order discretizations. When the solution is sufficiently resolved, higher-order schemes on coarser grids perform better (higher accuracy, less dissipative) for the same cost than lower-order scheme on finer grids. This result applies to both passive and reactive tracers, and is confirmed by quantitative analyses of truncation errors and smoothness of solution fields. To reduce oscillations in un-resolved regions, we develop a numerical filter that is active only when and where the solution is not smooth locally. Finally, we consider idealized simulations of biological patchiness. Results reveal that higher-order numerical schemes can maintain patches for long-term integrations while lower-order schemes are much too dissipative and cannot, even at very high resolutions. Implications for the use of simulations to better understand biological blooms, patchiness and other nonlinear reactive dynamics in coastal regions with complex bathymetric features are considerable.


\keywords{Ocean straits \and shelfbreaks \and biogeochemical dynamics\and high-order \and curved boundary \and quadrature \and ocean modeling \and Discontinuous Galerkin}
% \PACS{PACS code1 \and PACS code2 \and more}
% \subclass{MSC code1 \and MSC code2 \and more}
\end{abstract}
%\include{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%           INTRODUCTION          %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}\label{sec:intro}
 Accurate modeling of biogeochemical-physical ocean dynamics is required for multiple scientific and societal applications, covering a wide range of time and space scales. With the increased understanding of biogeochemical interactions 
\citep{lalli97,seav12,fennel04}, ecosystems models have substantially improved in the past decades \citep{fasham90,hofmann98,robinson99,hofmann02,lynch09}. Coupled biogeochemical-physical models have been used from coastal regions e.g. \citep{anderson05,spitz05,ji08,stow09} to basins and global ocean domains e.g. \citep{oschlies98,rothstein06,doney09}. However, in light of the strong nonlinearities observed in biological processes, an important subject that has been largely overlooked is the numerical requirements for such simulation studies. One of the major objectives of our work is to address such computational questions for reactive ocean tracers, directly including the latest advances in computational fluid dynamics e.g. \citep{chung02,ferziger02,lomax03,cebeci05,karniadakis05} and multiscale ocean modeling \citep{lermusiaux08}.
          
          Previous numerical ocean studies related to ours have primarily focused on passive or dynamic (density-related) tracer advections.  The most significant progress include the results  of \citet{hect95,hanert04,budgell07}, but none of these advances have dealt with higher order advection of reactive tracers on unstructured meshes with curved geometries. \citet{iskandarani05} applied and studied high-order schemes for  passive tracer and density dynamics in two-dimensions, including \citet{hect95}'s test and the gravitational adjustment of density in a channel of constant depth \citep{haidvogel99}, but they did not consider curved elements. \citet{levy01} assessed five different low-order finite volume advection schemes for biological modeling and found a 30\% difference in new production estimates, highlighting the need for careful numerical studies.  In \citet{bernard09}, high-order discontinuous Galerkin (DG) methods are used to solve tidal flows around shallow water islands with non-trivial geometries and using curved triangular meshes. Here, we are interested in biogeochemical tracers with possibly highly non-linear reactive or source terms, and we compare a set of low to high order schemes, both in time and in space. We employ the DG Finite Element Method \citep{cockburn98c}, using both straight and curved elements, and we study a varied set of numerical properties. As in  previous computational studies, we restrict our numerical analyses to two-dimensional (2D) flows, focusing on coupled dynamics in idealized straits.           
          
Our ultimate dynamics motivation is to allow quantitative 
simulation studies of fundamental nonlinear biological-physical dynamics in coastal regions with complex bathymetric features such as straits, sills, ridges and shelfbreaks. Such features strongly affect flows and, if they are shallow enough, one can expect biological responses in the euphotic zone. Multiple physical scales are possible, from rapid tidal effects to slow water-mass driven overflows, and biological resonances at some of these scales are likely. Our focus is on the numerical requirements prerequisite to such studies. Our work is partly inspired by our experience in coastal regions with complex geometries \citep{haley10}, especially with steep shelfbreaks such as the Massachusetts Bay and Stellwagen Bank \citep{besiktepe02}, Middle Atlantic Bight shelfbreak \citep{lermusiaux99}, Monterey Bay shelfbreak \citep{haley09}, Taiwan region shelfbreak \citep{lermusiaux10} and Philippine Archipelago Straits \citep{haley10}. The latter effort particularly motivated the present work, within the context of the Philippines Experiment (PhilEx) which is a five-year joint research project focused on interdisciplinary modeling, data assimilation, and dynamical studies in the straits regions of the Philippine Archipelago to better understand, model and predict sub-mesoscale and mesoscale physical and biogeochemical dynamics in complex regions. For realistic PhilEx simulations, we employ our MIT Multidisciplinary Simulation, Estimation, and Assimilation Systems \citep{mseas}. It includes a free surface hydrostactic ocean model over complex geometries with novel implicit schemes for telescoping nesting \citep{haley10}. This physical model is coupled to biological models \citep{besiktepe02}, forced with multiscale barotropic tides \citep{logutov08} and initialized with new objective mapping schemes specific for multi-connected domains \citep{arpitthesis,agarwal10}. The multi-resolution nested domains cover very shallow regions with strong tides, steep bathymetries and the deep ocean. The MSEAS system was employed in real-time, assimilating data sets from ships, gliders and satellite remote sensing, and issuing daily physical-biological forecasts with dynamical descriptions and adaptive sampling guidance \citep{philex}. The complex, nonlinear and multiscale biology in the region confirmed the need for the present computational studies.   
          
          Our work is part of an incubation for the next-generation of ocean modeling systems, focusing on key numerical questions for biogeochemical dynamics. The biological model we employ is based on \citet{flierl02,lisathesis} and \citet{ueckermann09}. We restrict ourselves to a relatively simple model to focus on the numerics. However, the model is complex enough to reveal important characteristics and to complete a large number of parameter sensitivity studies which we can synthesize. We study three biological regimes, one with single stable points at all depths and two with stable limit cycles. We examine interactions that are dominated by the biology, by the advection, or that are balanced. We also consider idealized simulations of biological patchiness which is commonly observed in the coastal ocean. For these regimes and interactions, we study a wide range of temporal and spatial discretizations. In what follows, we give our dynamical problem statement, definitions, and notation in \sref{defs}. Our new numerical schemes and discretization are formulated and studied in \sref{nummeth}. The results of our varied numerical and scientific investigations are described in \sref{apps}. Finally, our conclusions are stated in \sref{concl}. 

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%      Numerical Methodology      %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dynamical Problem statement, definitions, and notation}\label{sec:defs}
 	 	
\subsection{Dynamical Problem Statement}

The biological dynamics are governed by the following advection-diffusion-reaction equations:
\begin{eqnarray}
	\frac{\partial \Phi}{\partial t} + \nabla \cdot (\mathbf u \Phi) -  \kappa \nabla^2 \Phi  &=& \mathbf S(\Phi,\mathbf x, t), ~~ 			\mathrm{in} ~ \Omega
	\label{eq:gov}
\end{eqnarray}
with boundary conditions
\begin{eqnarray}
	\begin{array}{rcll} \displaystyle
		\Phi &=& \mathbf g_{D},              &~~ \mathrm{on} ~ \Gamma_{D} \\ \displaystyle
		(\mathbf u \Phi -  \kappa \nabla \Phi) \cdot \nrml  &=&  \mathbf g_{N}, &~~ \mathrm{on} ~ \Gamma_{N}
	\end{array}
	\label{eq:govbcs}
\end{eqnarray}
where $\Phi(\mathbf x,t)=[\phi^1(\mathbf x,t),\ldots,\phi^{N_c}(\mathbf x,t)]$ is the vector of $N_c$ biological components, $\mathbf u$ is the prescribed velocity field, $\kappa$ is a positive diffusivity coefficient, $\mathbf S(\Phi,\mathbf x, t)$ is the biological reaction terms, and $\mathbf g_{D},~  \mathbf g_{N}$ are the boundary conditions for the Dirichlet and Neumann boundaries respectively. Equations \ref{eq:gov}-\ref{eq:govbcs} are solved on the domain $\Omega \in \mathbb R^d$, where $d$ is the dimension of the problem, with boundary $\partial \Omega = \Gamma_{D} \cup \Gamma_{N}$ such that $\Gamma_{D} \cap \Gamma_{N} = \emptyset$.

Since we are interested in strait dynamics, for the flowfield $\mathbf u$, we assume that earth rotational effects are negligible, which is true if the ratio of the strait width to the Rossby radius is small \citep{pedlosky87,signell89,benoit94,bourgault04}. Additionally, for uniform geometry across the strait with a rigid lid approximation, a small Froude number, and a homogeneous density, the velocity field can be approximated as a potential flow field. A similar setup was used by \citet{signell89} for tidal flows. The potential velocity $\mathbf u$ is obtained by solving for the stream function
\begin{eqnarray}
	\nabla^2 \psi  &=& 0, ~~ \mathrm{in} ~ \Omega \label{eq:gov2}\\
	\mathbf u      &=& \nabla \times \psi	
\end{eqnarray}
with boundary conditions
\begin{eqnarray}
	\psi &=& h_{D},              ~~ \mathrm{on} ~ \partial \Omega.
	\label{eq:govbcs2}
\end{eqnarray}

\subsection{Finite Element Definitions and Notation}\label{sec:notat}
We discretize \eref{gov}-\eref{govbcs2} using the Discontinuous Galerkin (DG) Finite Element method. The first reported use of DG  FEM was by \cite{reed73} where DG was used to solve the steady-state neutron transport equation. However, DG drew little attention until a series of papers \cite{cockburn89a,cockburn89b,cockburn90,cockburn98}, where the Runge-Kutta DG methods were described. The extension of DG to higher order derivatives by \cite{bassi97} made the method applicable to solving advection-diffusion equations, which can be extended to solving the Navier Stokes equations. Since the late 90's, DG has seen a number of realistic applications in aerospace, solid mechanics, and electromagnetism to name a few. For a review on the use of such schemes in next-generation physical ocean models, we refer to \citet{pain05,slingo09} and \citet{ueckermann09}. 

In this section, we first describe the notation used for the domain and the discrete elements. Then we define the notation used for the solution on the element interior, and on the element interfaces.  Next, we described the discontinuous polynomial spaces and the necessary inner products for the DG discretization, followed by the discontinuous finite element space and inner product for the HDG method. A set of terms is also defined.

The basic domain notation is illustrated in \fref{notat}. We let $\mathcal T_h = \cup K_i$ be a finite collection of non-overlapping elements, $K_i$, that discretizes the domain $\Omega$, where $h$ denotes the characteristic size of an element. Also, let $\partial \mathcal T_h = \{ \partial K : K \in \mathcal T_h\}$ be the set of interfaces of all elements. For two elements $K^+$ and $K^-$ belonging to $\mathcal T_h$, we define $e=\partial K^+ \cap \partial K^- \neq \emptyset$ as the unique interior interface between elements $K^+$ and $K^-$. For a single element $K$ belonging to $\mathcal T_h$, $e=\partial K \cap \partial \Omega \neq \emptyset$ is a boundary interface. Let $\varepsilon_h^\circ$ and $\varepsilon_h^\partial$ denote the set of  unique interior and boundary interfaces, respectively, such that $\varepsilon_h =  \varepsilon_h^{\partial} \cup \varepsilon_h^\circ$. We note that in the interior $\partial {\mathcal T_h}$ contains two interfaces, $\partial K^+$ and $\partial K^-$, at the same location (one for each element), whereas the set $\varepsilon_h$ only contains a single interface, $e$, at the same location. 

\begin{figure}
	\centering
	\includegraphics[width=.55\columnwidth]{fig1-notat}%
	\caption{Notation for domain and triangular elements.}%
	\label{fig:notat}
\end{figure}

$K^+$ and $K^-$ have outward pointing normals $\nrml^+$ and $\nrml^-$ respectively. We then let vector and scalar quantities $(\mathbf q^{\pm},u^{\pm})$ be the traces of $(\mathbf q,u)$ on the interface $e$ from the interior of $K^\pm$. The mean value $\mean{ . }$ and jumps $\jmp{.}$ on the interior interface $e \in \varepsilon_h^\circ$ for scalar and vector quantities are then defined as
\begin{eqnarray*}
\begin{array}{rclcrcl}
	\mean{ \mathbf q } &=& (\mathbf q^+ + \mathbf q^-)/2 && \mean{u} &=& (u^+ + u^-)/2 \nonumber \\ \displaystyle
	\jmp {\mathbf q \cdot \nrml } &=& \mathbf q^+ \cdot \nrml^{+} + \mathbf q^- \cdot \nrml^- && \jmp{ u\nrml } &=& u^+\nrml^{+} + u^-\nrml^{-}.
	\end{array}
\end{eqnarray*}
On the set of boundary interfaces $e \in \varepsilon_h^\partial$,  (with outward facing normal $\nrml$ on $\partial \Omega$), we set
\begin{eqnarray*}
\begin{array}{rclcrcl}
	\mean{ \mathbf q } &=& \mathbf q && \mean{u} &=& u \nonumber \\ \displaystyle
	\jmp {\mathbf q \cdot \nrml } &=& \mathbf q \cdot \nrml  && \jmp{ u \nrml } &=& u\nrml.
	\end{array}
\end{eqnarray*}
since here $\mathbf q$ and $u$ are single-valued. Note that the jump in a vector is a scalar (involving only the normal component of the vector), whereas the jump in a scalar is a vector. Additionally, the jump will be zero for a continuous function. 

The main difference between continuous Galerkin (CG) and DG lies in the approximation subspaces used. DG uses bases that are in normed space $L^2(\Omega)$ while CG uses bases that are in the Hilbert space $H^1(\Omega)$, that is, the function has to be continuous across elements. For a function $f(\mathbf x)$ to be in $L^2(\Omega)$, it has to satisfy $\int_\Omega f(\mathbf x)^2  d\Omega < \infty$, whereas a function in $H^1(\Omega)$ has to belong to a smaller space satisfying $\int_\Omega f(\mathbf x)^2 + \nabla f(\mathbf x) \cdot \nabla f(\mathbf x)  d\Omega < \infty$. Let $\mathcal P^p(D)$ denote the set of polynomials of maximum degree $p$ existing on a domain $D$. For example, we will be using $p2$ to denote a second degree polynomial basis, which will result in a 3$^{\mathrm{rd}}$ order accurate scheme. We introduce the discontinuous finite element spaces defined as
\begin{eqnarray*}
	W_h^p&=&\left\{w \in L^{2}(\Omega):w \mid_K \in \mathcal P^p(K), \forall K \in \mathcal T_h \right\} \\
\mathbf	V_h^p&=&\left\{\mathbf v \in (L^{2}(\Omega))^d:\mathbf v \mid_K \in  (\mathcal P^p(K))^d, \forall K \in \mathcal T_h \right\} 
\end{eqnarray*}
where $W_h^p$ is a scalar space, $\mathbf V_h^p$ is a vector space of dimension $d$, and $L^2(D)$ is the space of square integrable functions $f(\mathbf x)$ such that $\int_D f(\mathbf x)^2  dD < \infty$ on domain $D$.

Finally we define the inner products over continuous domains $D \in \mathbb R^d$ and $\partial D \in \mathbb R^{d-1}$ as
\begin{eqnarray}
\begin{array}{rclcrcl} \displaystyle
	(\mathbf q, \mathbf v)_D &=& \int_D \mathbf q \cdot \mathbf v dD 
	  && (u, w)_D &=& \int_D uw dD \\ \\\displaystyle 
	\langle \mathbf q, \mathbf v\rangle_{\partial D} &=& \int_{\partial D} \mathbf q \cdot \mathbf v d{\partial D} &&
	 \langle u, w \rangle_{\partial D} &=& \int_{\partial D} uw d{\partial D} \\
\end{array}
\end{eqnarray}
for vector functions $\mathbf q, \mathbf v$ and scalar functions $u,w$. Over discontinuous domains we also define
\begin{eqnarray}
	(q, v)_{\mathcal T_h} = \sum_{K \in \mathcal T_h} (\mathbf q, \mathbf v)_K,
	   ~~~~\langle u, w \rangle_{\partial \mathcal T_h} = \sum_{K \in \mathcal T_h} \langle u,w \rangle_{\partial K},
\end{eqnarray}
for vector or scalar functions $q, v$ defined on $\mathcal T_h$, and $u,v$ defined on $\partial \mathcal T_h$. 

 To use the Hybridized Discontinuous Galerkin (HDG) framework for solving \eref{gov2}-\eref{govbcs2}, we will require the traced finite element space existing on the interfaces $\varepsilon_h$
\begin{eqnarray*}
	M_h ^p &=& \{\mu \in L^2(\varepsilon_h) : \mu |_e \in \mathcal P^p (e), \forall e \in \varepsilon_h \}.
\end{eqnarray*}
We also set $M_h ^p(g_D) = \{\mu \in M_h ^p : \mu = \mathsf P g_D~\mathrm{on} ~\Gamma_D \}$, where $\mathsf P$ is the $L^2$ projection into the space $\{\mu|_{\partial \Omega} \forall \mu \in M_h^p \}$. Note that $M_h^p$ is continuous on the interface, $e$, shared by $K^+$ and $K^-$, but discontinuous at the borders between different interfaces. 
We will also require the additional inner product on this discontinuous domain
\begin{eqnarray}
	   \langle \mu, \eta \rangle_{\varepsilon_h} = \sum_{e \in \mathcal \varepsilon_h} \langle \mu,\eta \rangle_e 
\end{eqnarray}
for vector or scalar functions $\mu, \eta$ defined on $\varepsilon_h$. 

\subsection{Comparing numerical codes: Defining Efficiency, Accuracy, and Performance}\label{sec:perform}
To be clear, we use the term ``efficiency'' or "cost" to refer exclusively to the computational resources (elapsed-time, memory) required for a simulation, and we do not use ``efficiency'' to imply any degree of correctness of the solution. We reserve the term ``accuracy'' to refer to the correctness of the solution. Finally, here we also use the term ``performance'' as the combined consideration between efficiency and accuracy \citep{chapra06}. 

Comparing different numerical schemes is not a straight-forward task (see \citet{kubatko09}). First, results are not universally applicable and are generally problem dependent. After focusing on a particular class of problems, the usual approach is to fix the computational efficiency of both schemes and then compare the accuracy, or vice versa. The scheme that performs better will then have a superior accuracy, since the efficiency will be the same for both. However, the efficiency of the scheme is dependent on its implementation, as well as the computer architecture on which the simulations are performed. A simple approach, then, is to fix the number of degrees of freedom (DOFs) of the different schemes, that is, having the same number of unconstrained parameters in both schemes. Because the DOFs are related to computational efficiency, this approach is useful for comparing similar numerical schemes with different implementations. However, it is not a good approach when comparing different numerical schemes where the computational cost per DOF is inherently and significantly different between the schemes, which is the case for comparisons between high order and low order schemes. Finally, conclusions drawn about the performance is also dependent on the particular definition of accuracy. The accuracy is normally defined in terms of a quantity useful to a particular researcher. Thus, researchers with different quantities of interest may draw different conclusions about the performance of a scheme. We address the efficiency issue by presenting results for multiple efficiencies, and we address the accuracy issue by using generic global error measures (see \sref{errors}) and by using difference plots. 

\section{Numerical Methodology}\label{sec:nummeth}
In this section we first derive the basic FE formulations  for \eref{gov}-\eref{govbcs2}, using the notation from \sref{notat}. Starting with equations \eref{gov}-\eref{govbcs}, we multiply each biological component by test function $w$ and integrate over the domain. We seek approximations $\Phi_{h}=[\phi_h^1 \ldots \phi_h^{N_c}]$ of $\Phi$ such that for all $K \in \mathcal T_h$, 
\begin{eqnarray}
	\left(\frac{\partial \Phi_{h}}{\partial t},w \right)_{K} + \left(\nabla \cdot (\mathbf u \Phi_{h}),w\right)_{K} -  \left(\kappa 	\nabla^2 \Phi_{h},w\right)_{K}  &=& \left(\mathbf S(\Phi_{h},\mathbf x,t),w\right)_{K},  \forall w \in  \mathcal P(K).
	\label{eq:DGgov}
\end{eqnarray}
where each component $\phi_{h}^i \in W_h^p$. We set $\phi_{h}^i = \sum_{j=1}^{N_P} \phi^i_{j}(t) \theta_j(\mathbf x)$ where $\phi_{j}^i(t)$ are $N_P$ time varying coefficients with $N_P$ corresponding spatial basis functions, $\theta_j(\mathbf x) \in \mathcal P(K)$. For convenience we use Einstein summation notation $\Phi_{h} = {\Phi_{j}} \theta_j$, where the sum over the repeated index $j$ is implied and we let $\Phi_j \equiv \Phi_{j}(t)$,  $\theta_j\equiv\theta_j(\mathbf x)$. 

To obtain our finite-element formulation of equations \eref{gov2}-\eref{govbcs2}, we multiply by test functions $w$ and $\mathbf v$ and integrate over the domain. We seek approximations $\psi_h \in W_h^p$ of $\psi$ and $\mathbf u_h \in \mathbf V_h^p$ of $\mathbf u$ such that for all $K \in \mathcal T_h$,
\begin{eqnarray}
	\left(\nabla^2 \psi_{h},w\right)_{K}  &=& 0, ~~~~~~~~~~~~~~~~~~~~\forall w \in  \mathcal P(K), 	\label{eq:DGgov2}\\
	\left(\mathbf u_{h},\mathbf v \right)_{K} &=& \left(\nabla \times \psi_h,\mathbf v \right)_{K}, ~~~ \forall \mathbf v \in  			
		 (\mathcal P(K))^d
\end{eqnarray}

Next, we describe in detail the discretization of the source-terms in \sref{source}, since one of the novel aspects of our work is the high-order schemes for ocean biogeochemical simulations. Our spatial discretizations are derived in \sref{advectDG} for the advection terms, \sref{HDG} for the diffusive terms, and our temporal discretizations given in \sref{rk4}. The calculation of error norms is outlined in \sref{errors}. Finally, brief details of our implementation are presented in \sref{implement}, and our method for generating high-order curved meshes is described in \sref{meshgen}. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%  SOURCE TERM DISCRETIZATION  %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Source term discretization}\label{sec:source}

Quadrature-based integration approximates a definite integral by a weighted sum of function evaluations at discrete points (quadrature points). For example, Gaussian quadrature using $N_g$ points can exactly evaluate polynomials of degree $2N_g-1$. Quadrature-free methods avoid evaluating the weighted sum by using approximations and/or analytically evaluating definite integrals, and the solution to the definite integral is used directly in numerical implementations. For more details, see for example, \citet{hesthaven08a}. We further restrict the definition of ``quadrature-free'' to mean exact integration on the original polynomial basis. It is possible to use an expanded basis to evaluate the source terms, but additional cost is involved to interpolate the solution unto the higher degree basis.  For the source term discretization, either approach could be used, and we will examine the impact of this choice. 

To discretize the source terms with a quadrature-based approach, consider the source term for a single biological component $i$:
\begin{eqnarray}
	\left(S^i(\Phi_{h},\mathbf x_i,t),w_k\right)_{K}
	  &=& \left( S^i (\Phi_j \theta_j,\mathbf x,t ) ,w_k \right)_{K} \nonumber\\
	  &\approx& S^i\left (\Phi_j \theta_j(\mathbf x_i),\mathbf x_i,t \right ) w_k(\mathbf x_i) \omega_i \nonumber\\
	  &=& \mathcal W_{ki} \left [S^i\left (\Phi_j \Theta_{ji},\mathbf x_i,t \right ) \omega_i J_i \right]
\end{eqnarray}
where $\mathbf x_i$ are the locations of the $N_g$ quadrature points with corresponding weights $\omega_i$, $\mathcal W_{ki} = w_k (\mathbf x_i)$ are the values of the test functions evaluated at the quadrature points, and $\Theta_{ji}=\theta_j(\mathbf x_i)$ are the values of the basis functions evaluated at the quadrature points. To evaluate this integral numerically, we pre-compute the matrix $\mathcal W \in \mathbb R^{N_P \times N_g}$. Because a reference element is used, the corresponding Jacobians at the quadrature points $J_i$ also need to be calculated, and these are multiplied together with the weights ($\omega$). The integration over an element is performed as a matrix-vector multiplication with $\mathcal O (2N_g N_P)$ operations for a single biological component. In addition, we need to consider the interpolation of the values of $\Phi_h$ unto the quadrature points, $\Theta_{ji} \Phi_j$, resulting in an additional $\mathcal O (2N_g N_P)$ for each biological component. Finally, we have to evaluate the function describing the source terms at the quadrature points. The quadrature-based algorithm thus has a total of $\mathcal O (4N_cN_g N_P) + N_g C_s$ operations per element, where $C_s$ is the cost of evaluating the source terms. 
  
  Alternatively, to discretize the source terms with a quadrature-free approach using the original basis, we use
\begin{eqnarray}
	\left(S^i(\Phi_{h},\mathbf x,t),w_k\right)_{K}
	  &\approx& \left( S^i \left (\Phi_{j},\mathbf x_j,t \right )\theta_j  ,w_k\right)_{K} \nonumber\\
	  &=&  \left( \theta_j  ,w_k\right)_{K}  S^i \left (\Phi_{j},\mathbf x_j,t \right ) \nonumber\\
	  &=& \mathcal M_{kj}   S^i \left (\Phi_{j},\mathbf x_j,t \right )
	  \label{eq:sqf}
\end{eqnarray}
where $\mathbf x_j$ are the locations of the $N_P$ nodal points, and we can pre-compute the element local mass matrix $\mathcal M_{kj} = (\theta_j,w_k)_K $ such that $\mathcal M \in \mathbb R^{N_P \times N_p}$. Note, for straight-sided elements, the mass matrix computed on the reference matrix can be used and multiplied by the element local scalar Jacobian. However, on elements with curved boundaries, a mass matrix using the element local spatially variable Jacobian is computed. In the quadrature-free case, we approximate the integral by essentially fitting the function $S^i(\Phi,\mathbf x, t)$ on the polynomial space spanned by $\theta$. For a nodal basis, the coefficients of the basis are the values of the source terms evaluated at discrete points. Since these source term values at nodal points do not depend continuously on space, the second equality in \eref{sqf} follows.  This approach introduces an aliasing error since $S^i(\Phi,\mathbf x, t)$ may contain complicated functions which are not captured in $span\{\theta\}$. The operation count for this approach scales as $\mathcal O (2 N_c N_p N_p) + N_p C_s$ if we need to multiply through by the mass matrix. If the number of quadrature points are equal to the number of basis functions, $N_g = N_p$, and the cost of evaluating the source terms is small, then the quadrature-free algorithm is 2 times more efficient in terms of total operations. 

However, we can normally eliminate the mass matrix multiplication in front of the source term since there is also a mass matrix in front of the $\der {\phi_h^i}{t}$ term, $\mathcal M_{kj} \der {\phi^i_j}{t} = \mathcal M_{kj} S^i \left (\Phi_{j},\mathbf x_j,t \right ) \rightarrow \der {\phi_j^i}{t} = S^i \left (\Phi_{j},\mathbf x_j,t \right )$. Note, the mass matrix cannot be eliminated for all the terms in the partial differential equation (PDE), for example a matrix will remain in front of the discretized advection operator. However, the operation count for evaluating the source terms of the quadrature-free algorithm reduces to $N_p C_s$. Therefore, from an efficiency perspective, it is desirable to use a quadrature-free algorithm. 

The accuracy of the quadrature-based integration is limited by the quadrature rule used, and the number of quadrature points $N_g$. Choosing a greater number of quadrature points, any desirable accuracy can be obtained, at the cost of reduced efficiency. Conversely, the accuracy of the quadrature-free integration is limited by the order of basis function used, and potentially large errors can be introduced due to inexact integration. 
  
The biological source terms are the origin of non-linear and non-homogeneous dynamics for the whole PDE. They can lead to high-accuracy discretization requirements in the other terms of the PDE in both space and time. This is discussed in the next sections \sref{advectDG}-\sref{rk4}.
  
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%  ADVECTION TERM DISCRETIZATION  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
\subsection{Spatial discretization of Advection operators} \label{sec:advectDG}

Integrating the advection terms $\left(\nabla \cdot (\mathbf u \Phi_{h}),w\right)_{K}$ by parts and using the divergence theorem, we obtain the weak form of the advection operator
\begin{eqnarray}
   \left(\nabla \cdot (\mathbf u \Phi_{h}),w\right)_{K} &=& -\left(\mathbf u \Phi_{h},\nabla w\right)_{K} + \langle \widehat 		
   	{\mathbf u \Phi_{h}}\cdot \nrml,\nabla w\rangle_{\partial K},
\end{eqnarray}
where the formulation is complete once we specify the value of the flux $\widehat {\mathbf u \Phi_{h}}$. Here we use the upwind flux
\begin{eqnarray}
	\widehat {\mathbf u \Phi_{h}} \cdot \nrml &=& \mathbf u \cdot \nrml \mean{\Phi_{h}} - \frac{1}{2}|\mathbf u|\cdot\jmp{\Phi_{h} 	
		\nrml}
\end{eqnarray}

The same quadrature versus quadrature-free discussion is relevant to the advection terms. However, here we choose to use a quadrature-based scheme to ensure the accuracy of the advection part of the discretization. Since we are focusing on evaluating the accuracy of the source-terms, we do not want the additional consideration about the accuracy of the advection operators to complicate the discussion. Note, however, that considerable efficiency can be gained for the advection terms as well if a quadrature-free scheme is employed \citep{hesthaven08a}. We then let
\begin{eqnarray}
   -\left(\mathbf u \Phi_{h},\nabla w_k \right)_{K} + \langle \widehat {\mathbf u \Phi_{h}}\cdot \nrml,\nabla w_k 
   	\rangle_{\partial K}
   &\approx& -\mathcal D_{ki} \left [\left(\mathbf u_i \Phi_j \Theta_{ji} \right) \omega_i J_i \right] + \nonumber \\ &&
    \sum _{e\in K} \mathcal W^e_{ki} \left [(\widehat {\mathbf u_i \Phi_{j}}\Theta^e_{ji})\cdot \nrml_i  \omega^e_i J^e_i \right]  
\end{eqnarray}
where $\mathcal D_{ki} = \nabla w_k(\mathbf x_i)$, and the superscript $(.)^e$ indicates that the quantity is defined on an element interface, and $\mathbf u_i = \mathbf u_j \Theta_{ji}$ is the flow-field evaluated at the quadrature points. 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%  DIFFUSION TERM DISCRETIZATION  %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Spatial discretization of diffusive operators}\label{sec:HDG}
To discretize the diffusive operator $\left(\kappa \nabla^2 \Phi_{h},w\right)_{K}$, we follow standard practice with DG and consider the canonical problem
\begin{eqnarray*}
	-\nabla \cdot (\kappa \nabla) \phi  &=& s.
\end{eqnarray*}
Introducing the auxiliary variable $\mathbf q = -\kappa \nabla \phi$, this equation is re-written as two coupled first-order equations
\begin{eqnarray*}
	\nabla \cdot \mathbf q  &=& s	 \\
	\mathbf q + \kappa \nabla \phi  &=& 0.
\end{eqnarray*}
Multiplying these equations by the appropriate test functions, integrating by parts, and applying the divergence theorem, we obtain the finite element formulation for all $K \in \mathcal T_h$
\begin{eqnarray}
	-(\mathbf q_h, \nabla w)_K  + \langle \hat {\mathbf q_h}\cdot \nrml,w \rangle_{\partial K} &=& (s_h,w)_K	~~~\forall w \in 
		\mathcal P(K) \label{eq:hdg1}\\
	(\kappa^{-1} \mathbf q_h,w)_K - (\phi,\nabla \cdot \mathbf v)_K + \langle \hat{\phi_h}, \mathbf v \cdot \nrml \rangle_{\partial K}  &=& 0	 ~~~~~~~~~~~~~\forall \mathbf v \in \mathcal (P(K))^d\label{eq:hdg2}
\end{eqnarray}
where we have multiplied by $\kappa^{-1}$ so that it does not appear as part of the flux (or interface) terms. This formulation is complete once we specify the form of the flux terms at the interfaces, $\hat{\mathbf q}$ and $\hat{\phi_h}$. The diffusive fluxes for DG schemes are normally reported in the form
\begin{eqnarray}
	\mathbf{\hat q_h} &=& \mean {\mathbf q_h} - C_{11} \jmp{\phi_h \nrml} + \mathbf C_{12} \jmp{\mathbf q_h \cdot \nrml} \\
	\hat \phi_h &=& \mean{ \phi_h } - \mathbf C_{12} \cdot \jmp {\phi_h \nrml} - C_{22} \jmp{\mathbf q_h \cdot \nrml}
\label{eq:DGdifffluxes}
\end{eqnarray}

{\it Biogeochemical diffusive fluxes}: Using explicit time integration to solve \eref{gov}, we do not need to invert a matrix in \eref{DGgov}, in which case we utilize the Local Discontinuous Galerkin fluxes \citep{cockburn98b} 
\begin{eqnarray}
C_{11}= \tau, ~~~ \mathbf C_{12} = \frac{\nrml^\pm}{2}, ~~~ C_{22}=0,
\label{eq:LDG}
\end{eqnarray}
although many other choices exist. 

{\it Potential Flows}: For solving \eref{gov2}, a matrix inversion is required in \eref{DGgov2}, in which case we discretize the diffusive operators using the novel Hybridized Discontinuous Galerkin (HDG) method. For the full derivation of these equations, including more specific implementation details, refer to \cite{cockburn09,nguyen09}. The premise of HDG is to recognize that one can solve \eref{hdg1} - \eref{hdg2} locally on an element as long as the flux quantities are known. Within an HDG framework, the local element unknowns are parameterized in terms of a new variable $\lambda_h \in M^p_h(0)$, where the notation $M_h^p (0)$ refers to the space $M_h^p$ that is zero-valued on the boundaries of the domain. The fluxes are expressed as
\begin{eqnarray}
	\hat \phi_h &=& \left \{ \begin{array}[]{rl} 
						\mathsf P g_D, & \mathrm{on}~\varepsilon_h^{\partial}		\\
						\lambda_h , & \mathrm{on}~\varepsilon_h^{\circ}		
						\end{array} \right. \\
	\hat \mathbf q_h &=& \mathbf q_h + \tau (\phi_h - \hat \phi_h) \nrml, ~~ \mathrm{on}~\partial \mathcal T_h
\label{eq:hdguq}
\end{eqnarray}
where $\tau$ is a tune-able stabilization parameter, and $\mathsf P$ is, again, the $L^2$ projection into the space $\{\mu|_{\partial \Omega} \forall \mu \in M_h^p \}$. Now, once $\lambda_h$ is known, \eref{hdg1} - \eref{hdg2} can be solved efficiently on each element independently. What remains is an equation for $\lambda_h$, which can be found by enforcing continuity of the normal diffusive flux
\begin{eqnarray}
	\langle \jmp{\hat {\mathbf q}\cdot \nrml},\mu \rangle_{\varepsilon_h} = \langle g_N, \mu \rangle_{\Gamma_N}.
\label{eq:glbalhdg}
\end{eqnarray}
Note this is an equation with globally coupled unknowns. However, the number of unknowns is greatly reduced compared to the original system, since only unknowns on the interfaces $\varepsilon_h$ are involved. 

 This solution method involves three steps:
\begin{enumerate}
	\item The inversion of local operators on each element to form both the right-hand-side vector and the global matrix.
	\item The global solution to find $\lambda_h$.
	\item The local reconstruction of the solution on the element.
\end{enumerate} 
The local operations are efficient because inversions are done on matrices which are of dimension $\mathbb R^{(1+d)N_p\times (1+d)N_p}$. This procedure dramatically increases the efficiency of solving elliptic problems with DG where implicit time integration is required. Additionally, when the stabilization parameter for this choice of fluxes is chosen optimally ($\tau \sim \mathcal O(1)$), the optimal convergence rate of $\mathcal O(p+1)$ is obtained for both the gradient ($\mathbf q$) and the solution ($\phi$) \citep{nguyen09}. This allows a post-processing procedure that can produce a solution $u_h^*$ which converges at $\mathcal O(p+2)$.  

Reporting the fluxes for HDG in the standard form we have from \cite{nguyen09} the following
\begin{eqnarray}
C_{11}= \frac{\tau^+\tau^-}{\tau^++\tau^-}, ~~~ \mathbf C_{12} = \frac{1}{2} \left( \frac{\jmp {\tau \nrml}}{\tau^++\tau^-} \right), ~~~ C_{22}=\frac{1}{\tau^++\tau^-}.
\label{eq:HDGflux}
\end{eqnarray}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%% TEMPORAL DISCRETIZATION  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Temporal discretization}\label{sec:rk4}
	Motivated by strong nonlinearities in biogeochemical dynamics, an objective of this study is to evaluate effects of temporal discretization errors on the accuracy of numerical simulations. We investigate fourth-order, second-order, and first order schemes in time. Specifically, for the majority of this work, we will use the following four-stage low storage fourth-order-accurate Runge-Kutta scheme for explicit time integration
				\begin{eqnarray*}
				\Phi^a &=& \Phi(t) + \frac{\Delta t}{4} \left . \der{\Phi}{t} \right |_{\Phi(t)}\\
				\Phi^b &=& \Phi(t) + \frac{\Delta t}{3} \left . \der{\Phi}{t} \right |_{\Phi^a} \\
				\Phi^c &=& \Phi(t) + \frac{\Delta t}{2} \left . \der{\Phi}{t} \right |_{\Phi^b} \\
				\Phi(t+\Delta t) &=& \Phi(t) + {\Delta t} \left . \der{\Phi}{t} \right |_{\Phi^c}.
				\end{eqnarray*}
				It is implemented in a four-stage for loop, where the solution at the initial time $\Phi(t)$ is saved in a temporary array, and the array containing the solution is updated using the three intermediate $\Phi^{a-c}$ variables. 
				
				We evaluate the temporal discretization error by considering a 2$^\mathrm{nd}$-accurate explicit Runge-Kutta scheme
				\begin{eqnarray*}
				\Phi^a &=& \Phi(t) + \frac{\Delta t}{2} \left . \der{\Phi}{t} \right |_{\Phi(t)}\\
				\Phi(t+\Delta t) &=& \Phi(t) + {\Delta t} \left . \der{\Phi}{t} \right |_{\Phi^a},
				\end{eqnarray*}
				and the first-order-accurate explicit Euler scheme
				\begin{eqnarray*}
				\Phi(t+\Delta t) &=& \Phi(t) + {\Delta t} \left . \der{\Phi}{t} \right |_{\Phi^t}.
				\end{eqnarray*}
				
				In each case, $\der{\Phi}{t}$ is evaluated using the right-hand-side spatial PDE. 
				
				
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%%%% ERROR NORMS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%			
	
\subsection{Error Norm Calculation}\label{sec:errors}
Unless indicated otherwise, the global domain $L^2$ norm $\norm{e} = \left(\int_\Omega e^2 d\Omega\right)^{\frac{1}{2}}$ is calculated using the quadrature-based approach as described in \sref{source}. That is, the numerical solution is interpolated unto the quadrature points, the error $e = \phi_h-\phi$ is evaluated, and then multiplied by the quadrature weights and summed for an approximate integration. In some cases, we evaluate the global error using an interpolation approach (similar to quadrature-free) and this is mentioned when we do. In these cases, the error is evaluated at the nodal points, then the error is interpolated to the quadrature points, multiplied by quadrature weights, and summed. Where ambiguous, we indicate the quadrature-based error evaluation using $\norm{e}^{qp}$ (quadrature points), and the interpolated error evaluation using $\norm{e}^{nd}$ (nodal points).  

The infinity norm $\left|\!\left| e \right|\!\right|_\infty = \mathrm{max} |e|$ is calculated by evaluating the error at nodal points, and taking the maximum absolute value.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%% IMPLEMENTATION  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				
\subsection{Implementation}\label{sec:implement}
	The discretized equations were implemented for 2D triangular elements. Even though our codes are efficient, our implementation is not fully optimized, in particular, the higher order simulations would most benefit from further optimization (e.g.~see~\cite{lambrechts10}). This is fine for our purposes since if we find that our high-order implementations are more accurate for the same cost/efficiency than lower-order schemes, then further optimization would only accentuate this result.
	
	The correctness of our implementation is verified by performing convergence studies using analytical test cases on curved and straight geometries. The implementation of each discretized operator is verified both separately and collectively. Results of some of the convergence studies are shown in \sref{1Dbio} and \sref{flow}.
	
	 For the polynomial spaces restricted on each element $\mathcal P^p(K)$, a nodal basis with $N_p = \frac{(p+1)(p+2)}{2}$ nodal points is used in two dimensions, and $N_p = (p+1)$ nodal points in one dimension. The total number of DOFs can then be calculated by multiplying the number of nodal points by the number of elements. The node locations are chosen according the method described in \cite{hesthaven08a}. More specifically, we have the basis $\theta_j$ such that
	\begin{eqnarray*}
		\theta_j(\mathbf x_i) = \delta_{ij}
\end{eqnarray*}
where $\mathbf x_i$ are the nodal locations. This basis is constructed for the reference element $[0,~ 0],~ [1,~ 0],~ [0,~ 1]$, and an isoparametric coordinate mapping is used for arbitrary triangles. The isoparametric coordinate mapping is described through the element nodal locations in the problem reference frame. 

To integrate polynomials of degree $p \geq 10$, cubature rules (multi-dimensional quadrature rules) are constructed on triangles by using tensor products of one-dimensional Gauss quadrature rules. For $p<10$ we use tabulated rules from \citet{solin03} and \citet{strang73}. 

In all cases, we utilize the Galerkin approach, that is we choose the test functions to be the same as the basis functions, $w_j=\theta_j$. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%  HIGH ORDER MESH GEENRATION  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
					
\subsection{Higher-Order Mesh Generation}\label{sec:meshgen}

Since higher order DG schemes have more degrees of freedom per element,  a coarse mesh with large elements is required to keep a similar performance across discretizations. To obtain an accurate solution with a coarse, high-order discretization, it is necessary to use curved boundary interfaces, as will be demonstrated in \sref{advect}. Here we describe our new method for creating such a coarse, high-order curved mesh. 

When curving the boundary of an element, care needs to be taken because it is possible to create an element where two of the interfaces cross. The left triangular element shown in \fref{curvedmsh} has the true circular geometry crossing one of the straight interfaces. To avoid this situation, we need to ensure that
\begin{eqnarray}
	h < 2\rho(x) \sin(\theta),
\label{eq:hmin}
\end{eqnarray}
where $h$ is the length of the element side bordering the boundary, $\rho(x) = \frac{|1+f' (x)|^{3/2}}{|f{''}(x)|}$ is the radius of curvature of the boundary described by $f(x)$, and $\theta$ is the minimum angle of the two angles on the edge bordering the boundary. The element shown on the right side of \fref{curvedmsh} illustrates this limiting case for an equilateral triangle, but our condition \eref{hmin} is trivially extended to arbitrary triangles as shown by the dashed lines. 

Using our criteria \eref{hmin}, we define the minimum edge spacing on the boundary as $h_{min} =  2 \rho \sin(30^\circ) = \rho$. Then, we let the minimum edge length grow linearly by a certain percentage (fit to $12\%$ here) away from the boundary upto a specified minimum edge length. Using these criteria, we create coarse base-meshes, then uniformly refine these meshes to obtain finer discretizations. To create the meshes, we primarily used the free mesher {\it Distmesh} \citep{persson04}, but we also used {\it Gmsh} \citep{remacle09}. {\it Distmesh} uses an implicit geometry representation, that is, we define the geometry by a distance function that gives the distance between a queried point and the nearest boundary. Using {\it Distmesh}, we create meshes with straight sides. 

To curve the boundary interfaces, we use the same distance function provided to {\it Distmesh}, and numerically calculate the gradient of the distance function to the boundary. The normalized gradient vector provides the direction of translation, but to determine the magnitude of the translation a weight needs to be applied to the calculated distance. That is,  $p_1^{new} = p_1^{old} +\mathcal W d \frac{\nabla d}{|\nabla d|}$, where $d$ is the distance from the point $p_1^{old}$ to the boundary, and $\mathcal W$ is the weight. Now, points on the straight boundary interface are translated to the true, curved boundary with a weight 1, and points on interior interfaces are not translated, i.e. having weights 0. Points in the volume have weights defined by the same weighting functions used to create the nodal basis, that is
\begin{eqnarray*}
	\mathcal W(e_1) = \left(\frac{2\lambda_3}{2\lambda_3+\lambda_1}\right)\left(\frac{2\lambda_2}{2\lambda_2+\lambda_1}\right)
\end{eqnarray*}
 where the point is defined by the barycentric coordinates $\lambda_i$ corresponding to vertices $i$, and $e_1$ is the curved boundary interface defined by vertices 2 and 3. For details of this blending function, see \cite{hesthaven08a}. 

The base mesh with three mesh refinements is shown in \fref{meshes}, and details of the base mesh for a curved and straight mesh boundary are shown in \fref{meshesI}. Using our criterion $h_{min} = \rho$, the minimum theoretical edge length for our geometry, that is a Gaussian bump defined by $\tilde H(x) = e^{-{x^2}}$, is $h_{min}=0.25$. The mesh shown in \fref{meshes}a) has a minimum edge length of $h_{min}=0.2418$, close to the theoretical value.

\begin{figure}%
\centering
\includegraphics[width=0.45\columnwidth]{fig2-curvedmsh}%
\caption{Minimum triangle angle criterion \eref{hmin} demonstrated on a circle with equilateral triangles. $h_1=2\rho$ does not satisfy the criterion, $h_2=2/3 \rho$ satisfies the criterion, and $h_3 = 2 \rho \sin(\pi/3)$ demonstrates the limiting case. This result can be extended to arbitrary triangles as shown by the dashed lines.}%
\label{fig:curvedmsh}%
\end{figure}

\begin{figure*}%
	  a) \begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig3a-meshes_1} \end{center}
		b) \begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig3b-meshes_2} \end{center}
		c) \begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig3c-meshes_3} \end{center}
		d) \begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig3d-meshes_4}  \end{center}
		e) \begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig3e-meshes_5}  \end{center}
	\caption{a) The base mesh $(g1)$ with 350 elements.  b) First $(g2)$ (1,400 elements) c) second $(g3)$ (5,600 elements), d) third $(g4)$ (22,400 elements), and fourth $(g5)$ (89,600 elements) grid refinements. The more-refined meshes are used for lower-order schemes whereas less-refined meshes are used for higher-order schemes such that the cost of the two schemes are comparable.}%
	\label{fig:meshes}%
\end{figure*}

\begin{figure}%
	a)\begin{center}\vspace{-12pt} \includegraphics[width=0.45\columnwidth]{fig4a-meshesIa}\end{center}
	b)\begin{center}\vspace{-12pt} \includegraphics[width=0.45\columnwidth]{fig4b-meshesIb}\end{center}
\caption{Details of ($g1$) using a a) curved and b) straight mesh for a $p=8$ nodal basis.}%
\label{fig:meshesI}%
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%      Applications and Tests     %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Studies and Scientific Implications}\label{sec:apps}
Biogeochemical models may contain a large number of biological or chemical components \citep{hofmann02}. The simplest models often only use Nutrient, Phytoplankton, and Zooplankton as components, and are commonly called NPZ models. More complicated models \citep{besiktepe02} can be adaptive and contain many components. Each component requires the solution of an advection-diffusion-reaction (ADR) equation of the form \eref{gov}. The source terms describe the commonly nonlinear ``reactions", and may lead to stationary, periodic, or chaotic dynamics. For this numerical work, a non-dimensional version of a NPZ model \citep{flierl02} is used since it contains all characteristics required for our studies:
\begin{eqnarray}
\frac{\partial \phi_N^*}{\partial t^*} +\nabla \cdot (\mathbf u^* \phi_N^*) -\nabla \cdot \frac{1}{P_e} \nabla \phi_N^*
	&=& - \mathcal U^*e^{z^*/h^*} \frac{\phi_P^* \phi_N^*}{\phi_N^* + k_s^*} + d_P^* \phi_P^* + d_Z^* \phi_Z^* \nonumber \\ 
		&& + (1-a) g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*}) \label{eq:NPZ1}\\
\frac{\partial \phi_P^*}{\partial t^*} +\nabla \cdot (\mathbf u^* \phi_P^*) -\nabla \cdot \frac{1}{P_e} \nabla \phi_P^* 
	&=& \mathcal U^*e^{z^*/h^*} \frac{\phi_P^* \phi_N^*}{\phi_N^* + k_s^*} - d_P^* \phi_P^* \nonumber \\
		&&- g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*}) \\
\frac{\partial \phi_Z^*}{\partial t^*} + \nabla \cdot (\mathbf u^* \phi_Z^*) -\nabla \cdot \frac{1}{P_e} \nabla \phi_Z^* 
		&=& - d_Z^* \phi_Z^* + a g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*})
\label{eq:NPZ}
\end{eqnarray}
where $\phi_{(N,P,Z)}^* = \frac{\phi_{(N,P,Z)}}{\mathcal N_t}$, $u_x^*=\frac{u_x}{\bar {\mathbf u}}$,	$u_z^*=\frac{u_z L}{\bar {\mathbf u}H}$, $x^* = \frac{x}{L}$,	$z^* = \frac{z}{H}$,	$t^* = \frac{t}{\bar \tau}$, the parameters are explained in \tref{NPZparams}, the non-dimensional groups with values are given in \tref{NPZparamsvals} with $P_e$ the Peclet number and $D^*$ the aspect ratio, the subscripts $(.)_N, (.)_P, (.)_Z$ refer to Nutrients, Phytoplankton, and Zooplankton respectively, $\nabla=\der{}{x^*}+\der{}{z^*}$, and lowercase $z^*$ refers to the depth coordinate which is positive upward with $z^*=0$ at the surface. Note that not all three equations \eref{NPZ1}-\eref{NPZ} are required since the biological model satisfies the following conservation law for total nutrients, assuming a closed ocean system:
\begin{eqnarray}
\phi_N^* = 1 - \phi_P^* - \phi_Z^*.
\label{eq:NPZconserv}
\end{eqnarray}
The first equation \eref{NPZ1}, for example, could be eliminated in favor of \eref{NPZconserv}; however here we still use \eref{NPZconserv} to check the conservation of the numerical schemes. 
\begin{table}[b]
	\caption{NPZ equation parameter descriptions and units.}
	\label{tab:NPZparams}
	\centering
	{\footnotesize
	\begin{tabular}{lll}
		\hline 
		Parameter & Description  & [units]\\ 
		\hline
		\noalign{\smallskip}
		$\mathcal U$ 	& Phytoplankton uptake rate & [1/day]\\
		$k_s$ & Saturation concentration of phytoplankton  &[$\mu$mol/L]\\
		$d_P$ & Mortality rate of Phytoplankton  &[1/day]\\
		$d_Z$ & Mortality rate of Zooplankton  &[1/day]\\
		$g$ 	& Grazing rate of Zooplankton  &[L/($\mu$mol$\cdot$day]\\
		$a$ 	& Assimilation (efficiency) rate   &[]\\
		$h$ 	& e-folding depth for light (photosynthesis) & [m]\\
		$\nu$ & Parameter for Ivlev form of grazing function  &[$L/\mu \mathrm{mol}$]\\
		$\mathcal N_T$ & Total biomass   & [$\mu$mol/L] \\
		\noalign{\smallskip}
		\hline
		\noalign{\smallskip}
		$\bar {\mathbf u}$ & Average inlet velocity & [km/day] \\
		$H$ & Height of bathymetry & [m] \\
		$D$ & Total maximum depth  & [m] \\
		$L$ & Effective width of bathymetry  & [km] \\ 
		\uuline{$\kappa$} & Diffusion tensor (vertical and horizontal diffusion different) & [m$^2$/s] \\ 
		\noalign{\smallskip}\hline
	\end{tabular}}
\end{table}

\begin{table}
	\caption{Values of the dimensionless numbers entering the NPZ equations \eref{NPZ}, that are used in the examples for this manuscript. Bracketed triplets of  values correspond to the three bio cases [1,2,3]. The other values are the same for the three cases.}
	\label{tab:NPZparamsvals}
	\centering
	{\footnotesize
	\begin{tabular}{lll}
		\hline 
		\noalign{\smallskip}
		Parameter  & Value \\
		\noalign{\smallskip}
		\hline
		\noalign{\smallskip}
		$\mathcal U^* = \mathcal U \bar \tau $ & 7.5  \\
		$k_s^*=\frac{k_s}{\mathcal N_T}$ 			 &  [$\frac{1}{30}$,$\frac{1}{50}$,$\frac{1}{100}$] \\
		$d_P^*=d_P{\bar \tau}$ 								 &  0.2 \\
		$d_Z^*={d_Z}{\bar \tau}$ 							 &  1 \\
		$g_\nu^*=\frac{g \bar \tau}{ \nu}$ 		 &  12.5 \\
		$a^*=a$ 															 &  0.4 \\
		$h^*=\frac{h}{H}$ 									   &  0.34 \\
		$\nu^*=\mathcal N_T \nu$							 &  [$0.3,0.5,1$] \\
		\noalign{\smallskip}
		\hline
		\noalign{\smallskip}
		$P_e=\frac{\bar {\mathbf u} L}{\kappa}$ 	 &  $\infty$ \\
		$D^* = \frac{D}{H}$ 											 & 2 \\
		\noalign{\smallskip}\hline
	\end{tabular}}
\end{table}

The domain setup is depicted in \fref{biodomain} for the geometric parameter values given in \tref{NPZparamsvals}. With this setup, an upwelling of nutrients is created (see \sref{bioadvect}), and the study of idealized biological blooms, which may occur in straits or sills, can be studied. In total, we consider three sets of parameter values, differing by the non-dimensional parameters $\nu^*$ and $k_s^*$. In the absence of advection and diffusion, they lead to equations \ref{eq:NPZ} with at most one physically relevant steady state solution \citep{lisathesis}. The three sets of non-dimensional parameters $\nu^*$ and $k_s^*$ correspond to biological dynamics with single stable points at all depths (bio case 1: $k_s^*=1/30, \nu^*=0.3$), with stable limit cycles for depths around $z^* = 0.4-0.9$ and single stable points elsewhere (bio case 2: $k_s^*=1/50, \nu^*=0.5$), and stable limit cycles everywhere in the euphotic zone (bio case 3: $k_s^*=1/100, \nu^*=0.1$). The middle parameter values, bio case 2, correspond to those values used by \citet{flierl02}: they are idealized and not meant to represent a specific ocean region. We note that biological models with discontinuities in stable solutions are not always representative of nature. However, biology of interest is likely to have intrinsic oscillatory or chaotic time dependence, e.g.~\citet{flierl02}. For our purposes, we address these issues by considering three sets of parameter values and so cover a range of biological dynamics. To handle non-physical negative concentrations due to numerics, we use $\mathrm{max}(0,\phi^*_{(N,P,Z)})$ when evaluating the source terms. In the absence of advection, a timescale of $\bar \tau =1$[days] is used, while in the presence of advection, the advective timescale $\bar \tau_a = L/\bar {\mathbf u}$ is used, where $\bar {\mathbf u} = \int_{Inlet} |\mathbf u| dz$ is the average inlet velocity.


In our numerical study, we need to characterize the three dynamical regimes and their behaviour since these properties affect numerical errors. Specifically, for each dynamical regime, we study the three limiting balances of terms in \eref{NPZ1}-\eref{NPZ}: biological terms dominating; advection terms dominating; advection and biological terms balancing. When biological terms dominate, the advection is slow, and the problem reduces to a one dimensional problem, studied in \sref{1Dbio}. When advection is fast, the biology is unimportant, and we study this case in \sref{advect}, with the generated flowfield studied in \sref{flow}. Finally, the case where the advection and biological terms are approximately balanced is studied in \sref{bioadvect}, with the effect of biological patches demonstrated in \sref{patch}. Since the timescale of biology varies in depth, the advection and biological terms can only be exactly balanced for one depth. While this results in many choices of approximately balanced terms, we will focus on one parameter set where $\bar \tau_a = 12.5$[days]. For more details on biological dynamics in straits, we refer to \citet{lisathesis}. Finally, capabilities of numerical filtering for higher-order schemes are examined in \sref{filter}.

\begin{figure}%
\centering
\includegraphics[width=0.5\columnwidth]{fig5-biodomain}%
\caption{Test case domain with idealized strait bottom geometry described by $\tilde H(x) = He^{-\frac{x^2}{L^2}}$.}%
\label{fig:biodomain}%
\end{figure}

Results in this section will be reported for various grid and polynomial degree combinations, and the notation (grid number, polynomial degree) is used to denote this information. For example, $(g2,p4)$ refers to the second grid (\fref{meshes}b) with a 4$^\mathrm{th}$ degree basis function. We provide a table with the number of DOFs and computational time estimates for the test cases we completed \tref{timings}. In what follows, we show the results of \HO\ and $(g2,p5)$ for our high-order simulations and compare them to \LO. Normalizing by the average computational time of \LO\, these simulations have relative computational times of 0.34, 1, and 1.2 for  \HO, \LO, and $(g2,p5)$, respectively (see \tref{timings}). Following the discussion in \sref{perform} and \sref{implement}, we note that, in terms of efficiency, \LO\ and $(g2,p5)$ are comparable (in fact, $(g2,p5)$ would be cheaper if fully optimized, see \sref{implement}). $(g1,p6)$ is included because it is comparable in accuracy to \LO\, and it highlights the effect of under-resolution (here $g1$) when using higher-order schemes (here $p6$). 


\begin{table}%
\caption{Normalized run-times and DOFs for various grids/polynomial degree basis functions, for the simulations in \sref{bioadvect}. The times are normalized by the $(g4,p1)$ run-time, and numbers in parentheses are the DOFs. }
\label{tab:timings}
\begin{center}
\begin{tabular}{|c|ccccc|}
\hline
 Degree of basis & Grid 1 & Grid 2 & Grid 3 &Grid 4 &Grid 5 \\ \hline
 1        & 0.0014        & 0.014       & 0.12      & \textbf{1.0} & 8.2\\
          & (1,050)       & (4,200)     & (16,800)  & (67,200)     & (268,800) \\\hline
 2        & 0.007         & 0.057       & 0.1       & 4.2          & \\
          & (2,100)       & (8,400)     & (33,600)  & (134,400)    & \\\hline
 3        & 0.026         & 0.21        & 1.8       &              & \\
          & (3,500)       & (14,000)    & (56,000)  &              & \\\hline
 4        & 0.062         & 0.51        & 4.1       &              & \\
          & (5,250)       & (21,000)    & (84,000)  &              & \\\hline
 5        & 0.15          & \textbf{1.2}&           &              & \\ 
          & (7,350)       & (29,400)    &           &              & \\\hline
 6        & \textbf{0.34} & 3.1         &           &              & \\
          & (9,800)       & (39,200)    &           &              & \\ 
          \hline
          \end{tabular}
\end{center}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%  1D BIO SOURCE TERMS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\subsection{One-dimensional Biogeochemical Source terms studies}\label{sec:1Dbio}

In this section we first illustrate the convergence of our numerical implementation. Following this, we examine the numerical behaviour of the biological source terms using three tests: {\it Perturbations of steady-states, Vertical resolution}, and {\it High-order bases}.

{\it Numerical convergence in space and time:} Since an analytical solution to \eref{NPZ} does not exist, we verify the spatial implementation of the quadrature-free and quadrature-based source terms using the analytical test problem $\der{\phi}{t} - 2 \der{^2\phi}{z^2} = S(z,t)$ on $\Omega \in [-100,0]$ integrating until $t=\pi/200$, with solution $\phi = \sin(t)\cos^2(\pi/50z)$ (for an appropriately chosen $S(z,t)$). We use a sufficiently small time-step, such that the errors are dominated by the spatial discretization. The results are shown in \tref{1dconvergs}, with the norm of the error $e = \phi_h - \phi$ calculated as described in \sref{errors}. From \tref{1dconvergs}, we note that our implementation converges at the optimal rates for both the quadrature-based and quadrature-free treatments. While the error numbers are for a special case and not those for \eref{NPZ1}-\eref{NPZ}, they show that the solution using quadratures is more accurate than the solution without quadratures, and this result will be generally expected. 

We verify the implementation of the 4$^{th}$ order low-storage Runge-Kutta time integrator using the ordinary differential equation $\der{\phi}{t} = \phi$ on $\Omega \in [-100,0]$ integrating up to a time $t=1$, with solution $\phi = \phi_0 e^{t}$. Here we choose $\phi_0 = 1$, such that the spatial discretization does not affect the error. The results are given in \tref{1dconvergt}, from which we note that our implementation converges at the optimal rate for this low-storage Runge-Kutta scheme. While this test corresponds to exponential biological growth, as above, the error values are of course not those that would occur for \eref{NPZ1}-\eref{NPZ}.

\begin{table}
	\caption{Spatial convergence of 1D DG solver used to evaluate the source terms using $N_h$ elements. The $L^2$ norm (see \sref{errors}) of the error, $e = \phi_h - \phi$, is smaller for the quadrature-based scheme, but the order of convergence is the same for both. Order of convergence is computed in a standard way, e.g.~\citet{chapra06}.} 
	\label{tab:1dconvergs}
 	\begin{center}
 		\begin{tabular}{|lc|cc|cc|}
 			\hline 
 			&       & $N_h=10$ &  & $N_h=20$ &  \\ 
 			&Degree & $|\!|e|\!|_2$ & Order & $|\!|e|\!|_2$ & Order \\ 
  		\hline 
 			Quadrature-  &p=1 &5.550e-003 &1.9 &1.409e-003 & 2.0 \\ 
 			based        &p=2 &5.901e-004 &2.9 &7.491e-005 & 3.0 \\ 
       				     &p=3 &4.690e-005 &3.9 &3.019e-006 & 4.0 \\ 
             			 &p=4 &2.976e-006 &4.9 &9.673e-008 & 4.9 \\ \hline
 			Quadrature-  &p=1 &1.340e-002 &1.9 &3.435e-003 & 2.0 \\ 
 			free         &p=2 &1.068e-003 &2.9 &1.332e-004 & 3.0 \\ 
       				     &p=3 &7.464e-005 &3.9 &4.724e-006 & 4.0 \\ 
             			 &p=4 &4.282e-006 &5.0 &1.274e-007 & 5.1 \\ 
			\hline 
		\end{tabular} 
	\end{center}
\end{table}  
 
\begin{table}
	\caption{Temporal convergence of 1D DG solver using $N_t$ timesteps (different values of $N_t$ given only to show that the order does not vary with $N_t$, but the absolute error of course changes). Order is computed using \citet{chapra06}.} 
	\label{tab:1dconvergt} 
 	\begin{center}
 		\begin{tabular}{|l|cc|cc|cc|cc|}
 			\hline 
 			Integration & $N_t=16$  &       & $N_t=32$  &       & $N_t=64$  &        \\ 
 			Scheme      & $|\!|e|\!|_2$ & Order & $|\!|e|\!|_2$ & Order & $|\!|e|\!|_2$ & Order  \\ 
  		\hline 
 			RK4         &5.683e-006 &3.9    &3.646e-007 &4.0    &2.308e-008 &4.0     \\ 
			\hline 
		\end{tabular} 
	\end{center}
\end{table} 
 
{\it Perturbations of steady-states:} The purpose of these studies is to examine how the biological dynamics behave as a result of perturbations away from the steady-state. Perturbations will arise due to the forcing and dynamics, and due to numerical reasons in the more complicated tests in \sref{bioadvect}, and it is important to understand how these perturbations will affect the biological dynamics. All three different biological regimes were examined in these perturbation tests. We focus on the behavior of the one-dimensional dynamics for the time interval $t^*=[0,250]$, because this corresponds to the residence time of the biology for the dynamics in \sref{bioadvect}. We initialize all tests using a perturbed or unperturbed exact steady state, which can be found by setting $\der{\Phi^*}{t}+\nabla \cdot (\mathbf u^* \Phi^*) -\nabla \cdot \frac{1}{P_e} \nabla \Phi^*=0$ in \eref{NPZ1}-\eref{NPZ}. The steady state solution is perturbed by setting $\Phi_{(P,Z)perturb}^* = (1 + \gamma) \Phi_{(P,Z)steady}^*$, where $\gamma$ is some constant, and using \eref{NPZconserv}. Where required, we impose $\phi_Z^* + \phi_P^* < 1$, by setting $\Phi_{(P,Z)perturb}^* = \frac{\Phi_{(P,Z)steady}^*}{({\phi_Z^*+\phi_P^*})_{steady}}$. This initialization is done numerically by setting the value of the numerical solution equal to the calculated solution {\it at the nodal points}. 

First we ensure that the exact steady state solution can be maintained, and then we initialize with a perturbation from the exact steady state, and the results are reported in \tref{PertCheck}. For these runs we used 100, second-order accurate linear elements ($p=1$), which roughly corresponds to the resolution at the inlet for $(g5,p1)$. We find that the steady solution can be maintained for all cases up to machine precision for the quadrature-free implementation when evaluating the error at the nodal points. This is because we initialized the numerical simulation using the exactly calculated steady state {\it at the nodal points}. Note that the quadrature-based scheme has a smaller difference than the quadrature-free version when evaluating the error at the quadrature points, except for the case with stable limit cycles in the euphotic zone (bio case 3). Because the quadrature version evaluates the source term at the quadrature points, and the interpolation of the solution onto the quadrature points is not exactly at the analytical steady-state, the source-terms are non-zero, and the solution evolves. If the source-terms were polynomials of lower degree than the basis in the $z$ direction, this would not happen. 

Finally, \tref{PertCheck} gives a rough description of the dynamical properties of the equations. Here the norm of the initial difference, ${\norm{D_i}}$ should be compared to the norm of the final difference at quadrature points for the quadrature-based treatment, and at the nodal points for the quadrature-free implementation. For the case with single stable points (bio case 1), the initial difference of the perturbed solution to the analytical steady state is greater than the final difference, which indicates that the solution is approaching the calculated steady-state value. For the case with stable limit cycles in the euphotic zone (bio case 3), the final difference is greater than the initial perturbed difference, showing the solution is logically not approaching the steady state, but instead spiraling outward toward the stable limit cycles present at each depth. Additionally, plotting the solution (\fref{pertcheck}) profile for the largest perturbation, we can see that the perturbed solution tends towards the steady state for the entire column for bio case 1, only for the top part of the water column for bio case 2, and nowhere for bio case 3. Thus, the parameter set with limit cycles in the euphotic zone (bio case 3) has the most structure in the vertical, and will require the most resolution to model accurately. Also, numerical perturbations will be most important for bio case 3 because the differences will grow away from the calculated steady state, as opposed to decaying. 

\begin{figure}
a) \begin{center} \vspace{-12pt}
\includegraphics[width=0.6\columnwidth]{fig6a-zOrb_NT3_4}
\end{center}
b) \begin{center} \vspace{-12pt}
\includegraphics[width=0.6\columnwidth]{fig6b-zOrb_NT5_4}
\end{center}
c) \begin{center} \vspace{-12pt}
\includegraphics[width=0.6\columnwidth]{fig6c-zOrb_NT10_4}%
\end{center}
\caption{Solution profiles at all depths with $\gamma = 5\%$. Magenta crosses show the analytical steady-state solution, the thick black dashed lines show the initial condition, green circles show the profile at $t^*=250$, and thin blue lines show the profile at $t^*=125$. Plotted for biological dynamics with a) single stable points at all depths, b) stable limit cycles at bottom of euphotic zone, and c) stable limit cycles in entire euphotic zone. The quadrature-based solution is plotted at the quadrature points, whereas the quadrature-free solution is plotted at the nodal points.}%
\label{fig:pertcheck}%
\end{figure}


\begin{table}
	\caption{Difference between analytical steady state solution, and perturbed solution at $t^*=250$. Here $D=\frac{\Phi-\Phi_h}{\int_\Omega 1 d\Omega} \cdot 100\%$ is the percent error per area in the domain. The column $\norm{D_{i}}$ gives the initial difference, $D_q$ indicates using quadratures, $D_{qf}$ indicates quadrature-free, $\norm{.}^{qp}$ indicates the error evaluated at quadrature points, $\norm{.}^{nd}$ indicates the error evaluated at nodal points.} 
	\label{tab:PertCheck} 
 \begin{center}
 \begin{tabular}{|ccll||llll|}
 \hline 
 Stability & $\gamma \cdot 100\%$ & $|\!|D_{i}|\!|_2^{qp}$ & $|\!|D_{i}|\!|_2^{nd}$ & $|\!|D_q|\!|_2^{qp}$ & $|\!|D_{qf}|\!|_2^{qp}$ & $\norm{D_q}^{nd}$ & $\norm{D_{qf}}^{nd}$ \\ 
 \hline 
 Single        & 0    &0.361 &0.000 &0.296 &0.361 &0.117 & 0.000 \\ 
 stable        & 0.05 &0.361 &0.014 &0.296 &0.361 &0.117 & 0.002 \\ 
 points        & 0.50 &0.377 &0.126 &0.295 &0.363 &0.122 & 0.015 \\ 
               & 5.00 &0.906 &0.860 &0.318 &0.402 &0.208 & 0.141 \\  \hline
Stable limit   & 0    &0.364 &0.000 &0.352 &0.364 &0.015 & 0.000 \\ 
cycles at      & 0.05 &0.366 &0.017 &0.353 &0.365 &0.019 & 0.011 \\ 
bottom of      & 0.50 &0.410 &0.169 &0.377 &0.387 &0.116 & 0.112 \\ 
euphotic zone  & 5.00 &1.360 &1.29  &1.01  &0.974 &0.945 & 0.915 \\  \hline
Stable limit   & 0    &0.109 &0.000 &0.112 &0.109 &0.038 & 0.000 \\ 
cycles in whole  & 0.05 &0.111 &0.021 &0.637 &0.517 &0.736 & 0.619 \\ 
euphotic zone  & 0.50 &0.234 &0.205 &1.40  &1.20  &1.51  & 1.36 \\ 
               & 5.00 &1.880 &1.88  &2.39  &2.17  &2.51  & 2.37 \\ 
 \hline 
\end{tabular} 
\end{center}
\end{table} 
 
 {\it Vertical resolution:} By varying the resolution of the problem, we found that a minimum of 25 degrees of freedom were necessary to roughly capture the vertical structure of the biological model dynamics at the final time. For the tests in \sref{bioadvect}, \HO\ has approximately 21 degrees of freedom, indicating that it will be under-resolved.  

\begin{figure}
a) \begin{center} \vspace{-12pt}
\includegraphics[width=0.6\columnwidth]{fig7a-zOrbho_NT5_1}\end{center}
b) \begin{center} \vspace{-12pt}
\includegraphics[width=0.6\columnwidth]{fig7b-zOrbhoU_NT5_1}\end{center}%
\caption{Solution profiles for all depths at $t^*=500$ using a 15$^{th}$ degree polynomial and 3 elements with $\gamma = 5\%$ for dynamics with stable limit cycles at the bottom of the euphotic zone. As in \fref{pertcheck}, the magenta crosses show the analytical steady-state solution, thick black dashed lines show the initial condition, green circles show the profile at $t^*=250$, and thin blue lines give the profile at $t^*=125$. The solution is plotted at the quadrature points for the quadrature version, and at the nodal points for the quadrature-free (i.e. where the source terms are evaluated). a) Uses well-behaved (Gauss-Lobatto) nodal points, b) uses uniform nodal points.}%
\label{fig:1djumps}%
\end{figure}

{\it High-order bases:} We find that the quadrature-based treatment of the source terms results in large jumps of the solution between elements. This is illustrated in \fref{1djumps} after 500 time units of integration using a 15$^{th}$ degree polynomial and 3 elements. The problem is amplified when using a uniform nodal spacing, due to a larger interpolation error. Also, increasing the number of quadrature points used for integration did not solve this problem. The problem originates from the discontinuous jump in the solution, causing oscillations known as Gibbs phenomena. Note that both simulations are initialized in the same manner, but the Gibbs oscillations can only be ``seen'' when evaluating the initial condition at points other than the nodes. The quadrature-based integration, then, ``sees'' these oscillations because the source terms are evaluated at the quadrature points. Using the quadrature-free approach for this one-dimensional problem essentially decouples the vertical nodes, so numerically, the quadrature-free version does not ``see'' the oscillations. The Gibbs oscillations would have occurred in the quadrature-free scheme if the initialization was done at the quadrature-points instead. This example illustrates one of the drawbacks of using increasingly higher order schemes, that is, without special treatment, large oscillations occur for non-smooth functions. Using lower-order but on a finer discretization (more elements) can be a better strategy if special treatment is not used. This issue is further addressed in \sref{bioadvect} and \sref{filter}.

In this section we showed that, with our implementation, both the quadrature-based and quadrature-free treatment of the source-terms give accurate, convergent results (see \tref{1dconvergs}), although the absolute error of the quadrature-based implementation is smaller than the quadrature-free implementation. Then we showed that the analytical steady state solution could be maintained, and illustrated the dynamical behaviour of three different biological parameter sets through perturbations of the steady state solutions. With the vertical resolution tests we found a minimum of 25 degrees of freedom necessary to roughly capture the vertical solution features of our particular setup. Finally, we showed that oscillations can occur solely due to numerics for a high-order discretization. While the quadrature-based algorithm was shown to be more accurate, the oscillations at element interfaces and the added numerical cost need to be considered. The additional accuracy may be warranted when a bifurcation of the solution could occur, or when the solution is under-resolved. As a whole, a key result is that, for any numerical scheme, careful numerical studies should be performed in one-dimension to understand the potential errors arising from the nonlinear source term discretization before proceeding with advective models.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%  FLOW FIELD CONVERGENCE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Flow field convergence}\label{sec:flow}
	A potential flow-field is calculated by solving \eref{gov2} using HDG as described in \sref{HDG}. Once $\psi$ is found, we take $\mathbf u = [\psi_z,-\psi_x]$. The value of $\psi$ is specified on all boundaries. The top and bottom boundaries are specified as constants $\psi=\psi_{top}$ and $\psi=\psi_{bot}$. The inlet and outlet stream functions are specified to vary linearly in $z^*$, i.e. $\psi = \psi_{top}+z^*(\psi_{bot}-\psi_{top})/D^*$. For the advection-dominating cases in \sref{advect}, periodic boundaries are used, in which case the values of $\psi$ are equal at the inlet and outlet, with no need for boundary conditions. A flowfield specific to the grid and polynomial degree is used for all simulations. 

	We perform a convergence study on steady flowfields to verify that we are indeed obtaining near-optimal $\mathcal O(p+1)$ convergence \citep{cockburn09,nguyen09}, for the gradients of $\psi$ used for the flowfield $\mathbf u$. To evaluate the convergence, reference solutions using $(g5,p2)$ and $(g2,p8)$ are calculated, both giving similar results. The error was evaluated by considering the point-wise solution at all interior vertices of grid 1. That is, the solution was averaged across all element-local solutions touching the vertex and compared to the reference solutions. The point-wise error calculated using $(g2,p8)$ as the true solution is plotted in \fref{flowconverg} for multiple grids and polynomial bases. Using this point-wise error calculation, we obtained near-optimal rate of convergence. When the domain boundary nodes are included, non-optimal (smaller and larger) rates of convergence were found for $p>1$. This may be due to variations in the discretizations of the domain boundary for the different grids/polynomial bases. Similarly, when the $\norm{.}$ or $|\!|.|\!|_\infty$ norms are considered, we find near-optimal convergence for $p=1$, but not for $p>1$. When the domain boundary is not curved, the flowfields with $p>1$ have both the same rate of convergence and error magnitude as $p=1$ using any of the error metrics mentioned. This shows that curved boundaries are essential for an accurate high-order solution. Also, by using HDG, we have a velocity field defined on the same grid with the same order of accuracy as our complete solution scheme for the ADR equations. 
	
	Therefore, by using HDG with boundary-curved elements we obtain an accurate, high-order convergent potential flowfield for our geometry. 

\begin{figure}%
\centering
\includegraphics[width=0.6\columnwidth]{fig8-flow_converg}%
\caption{Illustrating the convergence of the flow-field error. The point-wise error calculated using $(g2,p8)$ as the true solution for multiple grids and polynomial bases. Note that HDG gives near-optimal convergence for the derivative quantities $\mathbf u = [\psi_z,-\psi_x]$. }%
\label{fig:flowconverg}%
\end{figure}
		
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%%%  TRACER ADVECTION OVER BUMP  %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
				
\subsection{Tracer advection over bump test case results}\label{sec:advect}

We study tracer advection without source terms to evaluate the behaviour of our numerical scheme when advection dominates. We do not illustrate here the convergence of the resolution (our results as shown in \sref{1Dbio}-\sref{flow}). Instead we focus on numerical advection artifacts that may affect the behaviour of the biology. We examine three cases: a uniform tracer; a tracer with a discontinuous jump in the vertical; and a tracer with a discontinuous jump in the horizontal and a linear horizontal gradient. The first case examines the divergence of the flowfield. The second case examines the behaviour of the scheme in the presence of a horizontal jump, which occurs at the edge of the euphotic zone for our choices of parameters. Note that because we initialize by setting the numerical initial condition at the nodes equal to an analytical initial condition with a jump, and the jump does not necessarily occur at element boundaries, the assigned numerical initial condition contains oscillations. The number of oscillations increases with the polynomial degree, but is present at all orders. Also, the straightness of the interface is affected by the grid resolution. The third case examines the behaviour of the scheme in the presence of a horizontal jump and horizontal gradient, which is a numerical test for frontal dynamics. Unlike the second case, the assigned initial condition for the third case does not contain oscillations because the jump occurs at element boundaries. All cases evolve in a periodic domain. In all cases the flow is left to right, and the duration of the simulation is determined by the time the mean inlet velocity would take to travel through the domain. One of the objectives of these tests is to compare schemes at the same overall cost, for example, a lower order scheme is used on a higher resolution grid. 

The results for the reference solution $(g1,p5)$, low order solution, \LO, and high order solutions, $(g2,p5)$ and \HO\ on curved meshes, are shown in \fref{advection}. In the top plot (case 1), we find that oscillations exceeding $10^{-4}$ are less localized for the low-order schemes than for the high-order schemes. From the middle and bottom plots (cases 2-3) we see that the higher order schemes have larger magnitude oscillations around the jumps (as expected from the initialization). However, the sharpness of the interfaces are comparable between $(g1,p5)$ and $(g2,p5)$, and also between $(g4,p1)$ and $(g1,p6)$. Running a lower order case ($g3,p1)$ with 16,800 DOFs (not shown) resulted in a much more diffuse solution compared to $(g1,p6)$, even though $(g1,p6)$ is as under-resolved. The higher-order discretization with fewer degrees of freedom is therefore less numerically dissipative, illustrating one of the advantages of a high-order scheme. 

The locations where the solution is outside of the intervals  [0.9999,1.0001], [0, 1], and [-2, 3] for case 1, 2, and 3 respectively are plotted in \fref{advecterr}. From this figure, we see	that the smallest magnitude and most localized errors for case 1 happens with high-order schemes using curved boundaries. The largest magnitude errors occur for high-order schemes using straight boundaries, and the least localized errors occur for low order schemes. From case 2, we see that for both high and low-order schemes, the initial oscillations do not remain close to the jump, but spread through the domain. Although, the amplitude of the radiated oscillations are at least an order of magnitude smaller than the initialized oscillations. In case 3 we see that the high and low-order schemes have similar performance, both developing numerical oscillations with the same order of magnitude around the jump. The absolute magnitude of the oscillations are approximately half the size with low-order compared to high-order, and can be explained by the greater number of degrees of freedom with the low-order case. From these plots we note that the magnitude of the oscillations are within 20\% of the solution, and this could have a significant effect on the biology. This illustrates that discontinuous or highly non-smooth functions caused by physics or biology are problematic. In the context of this work, the jump will be smoother for the test cases in \sref{bioadvect} where we examine the case for approximately balanced advection and biological source terms. Alternatively, a shock-capturing (see for example \cite{hoteit04,persson06,krivodonova07}), filtering (for example \cite{hesthaven08b}), or postprocessing technique could be used to handle the discontinuity (see for example \cite{cockburn03,qiu05}). In \sref{filter} we illustrate how a filtering approach can be used to damp oscillations for high-order schemes. Note that special treatment is required for both high and low-order schemes, with a slope-limiting procedure more often used for the latter. 

We note that the large oscillations in case 1 for the high-order scheme happen only down-stream from the peak. This is because we do not ensure that the definition for the discrete divergence is the same in the equations calculating the potential flow field and the tracer/biological ADR equations (see for example \citet{dawson04}). Therefore, we expect to see discretization errors appearing in the ADR discrete divergence operator, which are advected downstream. Examining the normalized discrete divergence as defined in the ADR equations,  $\left(\frac{\nabla\cdot \mathbf u}{|\mathbf u|}\right)$, of curved and straight meshes with \HO\ and \LO, we found that largest divergence errors occurs for the higher-order scheme on the straight mesh. The normalized divergence, in this case, was of $\mathcal O(1)$ near the peak, and $\norm{\frac{\nabla\cdot \mathbf u}{|\mathbf u|}} \approx 0.56$ for the domain. The normalized divergence for the low-order and high-order curved schemes were of $\mathcal O(10^{-2})$ near the peak, with the maximum error for the low-order scheme approximately half the size of the high-order scheme. However, the error for the high-order scheme was more localized, resulting in a lower volume-averaged error $\norm{\frac{\nabla\cdot \mathbf u}{|\mathbf u|}} \approx 0.044$, than that of the low-order scheme which was $\norm{\frac{\nabla\cdot \mathbf u}{|\mathbf u|}} \approx 0.11$.  This reveals that it is important to use a curved mesh for higher-order schemes: the resultant flow field is then less numerically divergent than a flow field solved with a low-order finely resolved scheme. Without using a curved mesh, the divergence can be of the same order as the velocity near the geometry, which may excite non-physical biological dynamics downstream of the peak. In our case (in \sref{bioadvect}), the peak was beneath the euphotic zone and thus using a straight boundary representation would not cause  problems, however, in the general case this result cannot be overlooked. Of course, our advection scheme is conservative, but to also be constancy preserving, the numerical flow field needs to be discretely divergence-free. 

 Since the \HO\ simulations took less computational time than \LO, and because the $(g2,p5)$ simulation had a sharper interface, these results suggest that higher-order schemes performs better than the lower-order schemes for the advection dominated case, as long as curved boundaries are used. 

\begin{figure*}%
%\centering
\begin{center} \includegraphics[width=0.95\textwidth]{fig9-ADVECT_grouped2} \end{center}
\caption{Solution of three passive tracer cases, on a) $(g5,p1)$ with 268,800 DOFs, b) \LO\ with 67,200 DOFs, c) $(g2,p5)$ with 29,400 DOFs, and d) \HO\ with 9,800 DOFs with curved boundaries. Top plots shows advection of a constant tracer, middle is advection of tracer with vertical jump, and bottom is advection with a horizontal jump and horizontal gradient. Note that b) \LO and c) $(g2,p5)$ have similar costs while d) \HO\ is much cheaper and under-resolved (see Table \tref{timings} for costs).}%
\label{fig:advection}%
\end{figure*}

\begin{figure*}%
%\centering
\begin{center} \includegraphics[width=0.95\textwidth]{fig10-ADVECTg5_errgrp}\end{center}
\vspace{-12pt}\caption{Field values outside initially specified fields with a) $(g5,p1)$, b) \LO , c) $(g2,p5)$ on curved mesh, d)  \HO\ on curved mesh, e) $(g2,p5)$ on straight mesh, and f) \HO\ on straight mesh. For each case a-e), the top plot shows advection of a constant tracer with interval [0.9999,1.0001], middle is advection of tracer with vertical jump and interval [0,1], and bottom is advection with a horizontal gradient and horizontal jump with interval [-2,3].}%
\label{fig:advecterr}%
\end{figure*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%%%  FULL NPZ EQUATIONS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	

\subsection{Full NPZ equations} \label{sec:bioadvect}
In this section we explore the case where the advection and biological source terms are approximately balanced. We examine effects of low-order and high-order temporal discretizations in \sref{timecomp}. In \sref{qqf}, we illustrate the difference between using a quadrature-based and quadrature-free scheme to discretize the nonlinear biological source terms. Finally, in \sref{holo}, we study effects of spatial resolution, through both grid resolution and polynomial degree.

We still study the three biological parameter sets: single stable points, stable limit cycles at the bottom of the euphotic zone, and stable limit cycles for the entire euphotic zone, as given in \tref{NPZparamsvals}. Since the timescale of biology varies in depth, the advection and biological source terms can only be balanced for one depth. While this results in many choices of approximately balanced parameter sets, we focus on one where $\bar \tau_a = 12.5$ [days]. For these tests, the inlet is specified as the steady state solution with a smoothed discontinuity. The discontinuity is smoothed by fitting it with a cubic polynomial which can be resolved on $(g1,p6)$. The fit is biased such that 3/4 of the polynomial is below the euphotic zone. For the outlet boundary we use $\der{\Phi}{\nrml}=0$.  Also these results are compared to a $(g5,p1)$ simulation, which is taken as the true solution. 

The final solution fields for the three different regimes of biological dynamics (from \tref{NPZparamsvals}) and using quadrature-based source terms for $(g5,p1)$ (the reference solution) is plotted in \fref{npzfnl}. The results show that idealized strait bathymetry effectively perturbs the biology away from the inlet conditions. The case with single stable points (bio case 1) adjusts back to the stable equilibrium, whereas the two cases with limit cycles show complex structures in the vertical. In all cases, a Phytoplankton bloom over the bump is observed.

To qualitatively evaluate the effect of refining the grid or polynomial degree, we show the solution field for Phytoplankton for $(g3,p1)$, $(g3,p2)$, $(g4,p1)$, and $(g4,p2)$ in \fref{bioconverg}, and these discretizations have 16,800, 33,600, 67,200, and 134,400 DOFs respectively. This figure shows that the solution is converging with increased resolution. More quantitative comparisons are completed next. 

\begin{figure*}%
\begin{center} \includegraphics[width=0.95\textwidth]{fig11-BIOg5_grp2}\end{center}
\vspace{-12pt}\caption{ Biological dynamics at $t^*=20$ (with $\bar \tau_a=12.5$[days]) using $(g5,p1)$ . Biolocial dynamics with a) single stable points, b) stable limit cycles for depths $z^* = 0.4-0.9$, and c) stable limit cycles in whole euphotic zone. This is the reference solution against which all other solutions are compared.}%
\label{fig:npzfnl}%
\end{figure*}

\begin{figure*}%
%\centering
a)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig12a-BIO_convg3p1}\end{center}
b)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig12b-BIO_convg3p2}\end{center}
c)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig12c-BIO_convg4p1}\end{center}
d)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig12d-BIO_convg4p2}\end{center}
\caption{Phytoplankton fields at time $t^*=20$ (with $\bar \tau_a=12.5$[days]), as computed using four different spatial resolutions and order of the FE scheme: a) $(g3,p1)$, (16,800 DOFs) b) $(g3,p2)$ (33,600 DOFs), c) $(g4,p1)$ (67,200 DOFs), and d) $(g4,p2)$ (134,400 DOFs). All fields are for biological dynamics with stable limit cycles in the euphotic zone (bio case 3 in \tref{NPZparamsvals}).}%
\label{fig:bioconverg}%
\end{figure*}

\subsubsection{Comparing low-order and high-order temporal discretizations}\label{sec:timecomp}

We compare the solutions using 4$^{\mathrm {th}}$ order Runge-Kutta, 2$^{\mathrm {nd}}$ order Runge-Kutta, and 1$^{\mathrm {st}}$ order Explicit Euler on $(g2,p4)$ for the biology with stable limit cycles in the euphotic zone (bio case 3). The differences of the lower order schemes compared to 4$^{\mathrm {th}}$ order Runge-Kutta at  $t^*=40$ is plotted in \fref{timediff} for the Phytoplankton field. Note, the timestep size for the 1$^{\mathrm {st}}$ order scheme is half of the 2$^{\mathrm {nd}}$ order scheme, such that the cost of the two are the same.  For this test case we used periodic boundary conditions. From the figure, we note that the major differences occur within the euphotic zone. The stable explicit timestep for the 2$^{\mathrm {nd}}$ order scheme is set by the Courant condition for the advection discretization, and since the largest velocity occurs in the smallest element for this discretization, the timestep size is approximately four orders of magnitude smaller than the biological timescale.  Therefore, it is expected that temporal errors in the source-term should be small even for the low order scheme. Nonetheless, we still observe differences between the 1$^{\mathrm {st}}$, 2$^{\mathrm {nd}}$ and 4$^{\mathrm {th}}$ order schemes. We found that the difference at $t^*=40$ is approximately two orders of magnitude larger than at $t^*=20$, which indicates the errors are growing quickly. For the 1$^{\mathrm {st}}$ order scheme, the maximum error is of $\mathcal O(1)$ at $t^*=40$. This suggests that a low-order time discretization may result in significant errors when long integration times or fast biological timescales are involved. For example, the latter occurs in coastal applications. As another example, for stiff biogeochemical source terms, \citet{burchard05} found that even 4$^{\mathrm {th}}$ Runge-Kutta integration is insufficient to maintain the non-negativity of the biological components. They suggest that positivity preserving Patankar-Runge-Kutta schemes should be used to obtain a non-negative, conservative solution. 


\begin{figure*}%
a)\begin{center} \vspace{-20pt} \includegraphics[width=0.95\textwidth]{fig13a-BIO_RK4EulDiff}\end{center}
b)\begin{center} \vspace{-20pt} \includegraphics[width=0.95\textwidth]{fig13b-BIO_RK4RK2Diff}\end{center}%
\caption{Temporal discretization differences for Phytoplankton field with stable limit cycles in euphotic zone at $t^* = 40$ using periodic boundary conditions and on spatial grid $(g4,p2)$. a) ``1$^{\mathrm {st}}$ order Euler'' minus ``4$^{\mathrm {th}}$ order Runge-Kutta'', and b) ``2$^{\mathrm {nd}}$ order Runge-Kutta'' minus ``4$^{\mathrm {th}}$ order Runge-Kutta''.}%
\label{fig:timediff}%
\end{figure*}

\subsubsection{Comparing quadrature-based and quadrature-free source terms}\label{sec:qqf}
In \sref{1Dbio} we found that the greatest difference between the quadrature and quadrature-free treatment of the source-terms occurred for the biological parameter set with stable limit cycles in the euphotic zone (i.e. bio case 3). Here we examine this case for full ADR dynamics using the ($g1,p6$) discretization. Note that we obtained the same results and conclusions with the ($g4,p1$) and $(g5,p2)$ discretizations (not shown). Plotting the difference (quadrature-free minus quadrature-based) of the solution in \fref{npzdiff2} for ($g1,p6$) we see that the largest differences occur near the outlet of the domain where the mesh solution is under-resolved. The quadrature-based solution is more accurate in the under-resolved region because the source-term integral is more accurately evaluated, and this was verified by comparing the errors of the two schemes. However, where the solution is sufficiently resolved, the quadrature-free and quadrature-based treatments of the source terms have similar accuracy, that is, they differ by approximately 0.1\%. From the one-dimensional studies, we found that the quadrature-free algorithm was less oscillatory at element interfaces than the quadrature-based algorithm, and we observed the same effect in these 2D simulations for $p>7$ on $g1$, although the difference between quadrature-free and quadrature-based was less drastic. The largest differences between the quadrature-based and quadrature-free schemes did occur at element boundaries, and the quadrature-based algorithm was more accurate when under-resolved. 

Using equation \eref{NPZconserv}, we verify the conservation of the scheme. The results for the quadrature-free and quadrature-based source terms were similar up to floating point precision. Also, we find that the conservation error is dominated by the flow field divergence error.  Therefore, the conservation properties of the source term discretization does not affect the choice between quadrature-free and quadrature-based algorithms. 

Because the quadrature-free and quadrature-based algorithms had similar accuracy in well-resolved regions, we recommend using the quadrature-free treatment in these regions because of the improved efficiency. However, when the solution is poorly resolved, the quadrature-based treatment of the source terms is more accurate. Now, depending on the total solution cost of a particular numerical scheme, a finer resolution quadrature-free scheme may be more efficient for the same accuracy than a quadrature-based scheme. 


\begin{figure*}%
%\centering
a) \begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig14a-BIO_QZ}\end{center}
b) \begin{center} \vspace{-12pt}  \includegraphics[width=0.95\textwidth]{fig14b-BIO_QFZ}\end{center}
c) \begin{center} \vspace{-12pt} 	\includegraphics[width=0.95\textwidth]{fig14c-BIO_QF_diffZ}\end{center}
\caption{Zooplankton fields at $t^*=20$ computed using ($g1,p6$) and a) quadrature-based source terms, b) quadrature-free source terms. c) The difference between the quadrature-free and quadrature-based source-term simulations. The biological dynamics used has stable limit cycles within the euphotic zone (bio case 3).}%
\label{fig:npzdiff2}%
\end{figure*}

\subsubsection{Comparing low-order and high-order spatial discretizations}\label{sec:holo}
\fref{npzdiffstable} shows the differences between the reference solution and the solutions using other grids and polynomial degrees in \fref{npzdiffstable} for Zooplankton. We see that the $(g2,p5)$ simulation has the smallest differences, and is therefore the most accurate. This is a key result since it indicates that when the solution is resolved, for the same cost/efficiency, a higher order scheme on a coarser grid performs better than a lower-order scheme on a finer grid. Results for the biological dynamics with limit cycles at the bottom of the euphotic zone (bio case 2) are similar, but for the biology with limit cycles in the entire euphotic zone (bio case 3), both high and low-order schemes are under-resolved for $x^*>7$, especially for the \HO\ scheme. The differences between a fine grid solution ($g5,p1$) and the low-order and high-order schemes are plotted in \fref{npzdiff} for Zooplankton. From \fref{npzdiff}, we note that the errors in the low-order scheme are more localized in the $x^*>7$ region. However, the differences for \LO\ and $(g2,p5)$ are similar in the $x^*>7$ region. The \HO\ scheme has the least localized and largest magnitude errors in the $x^*>7$ region. However, as plotted in \fref{npzdiffdet} where the solution is smooth and the biology has less structure in the vertical, both the high-order schemes are more accurate than the low order scheme. Particularly, note the solution near the surface for $x^*<7$ in \fref{npzdiffdet}.

\begin{figure*}%
  a)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig15a-BIO_NT3_g4_P1}\end{center}
  b)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig15b-BIO_NT3_g2_P5}\end{center}
  c)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig15c-BIO_NT3_g1_P6}\end{center}
\caption{Difference between Zooplankton fields at $t^*=20$ (with $\bar \tau_a=12.5$[days]) computed using ($g5,p1$) and a) ($g4,p1$) , b) ($g2,p5$), and c) ($g1,p6$). This shows the locations of the largest numerical errors for the high-order and low-order schemes. The biological dynamics used have single stable points at all depths (bio case 1).}%
\label{fig:npzdiffstable}%
\end{figure*}

\begin{figure*}%
%\centering
  a)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig16a-BIO_NT10_g4_P1}\end{center}
  b)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig16b-BIO_NT10_g2_P5}\end{center}
  c)\begin{center} \vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig16c-BIO_NT10_g1_P6}\end{center}
\caption{As \fref{npzdiffstable}, but for  the biological dynamics with stable limit cycles within the euphotic zone (bio case 3).}%
\label{fig:npzdiff}%
\end{figure*}

\begin{figure}%
%\centering
  a) \begin{center} \vspace{-12pt} \includegraphics[width=0.45\textwidth]{fig17a-BIO_NT10_g4_P1_Detail}\end{center}
  b) \begin{center} \vspace{-12pt} \includegraphics[width=0.45\textwidth]{fig17b-BIO_NT10_g2_P5_Detail} \end{center}
  c) \begin{center} \vspace{-12pt} \includegraphics[width=0.45\textwidth]{fig17c-BIO_NT10_g1_P6_Detail} \end{center}
\caption{As \fref{npzdiff}, but zoomed in the region above the bathymetry. The difference between Zooplankton fields using  ($g5,p1$) and a) ($g4,p1$), b) ($g2,p5$), and c) ($g1,p6$).}%
\label{fig:npzdiffdet}%
\end{figure}


We examine the error characteristics of these fields more closely by considering the truncated Taylor expansions of the true solution. By the mean value theorem, the truncation error for \LO\ is $\frac{h^2}{2!} \left(\der{}{x}+\der{}{y} \right)^2 \phi(\mathbf x_\eta)$ for some unknown point $\mathbf x_\eta$, and for $(g2,p5)$ and \HO\ these terms are $\frac{h^7}{6!}\left(\der{}{x}+\der{}{y} \right)^6 \phi(\mathbf x_\zeta)$ and $\frac{h^7}{7!}\left(\der{}{x}+\der{}{y} \right)^7 \phi(\mathbf x_\xi)$ for unknown points $\mathbf x_\zeta$ and $\mathbf x_\xi$, where $h$ is the characteristic discretization length of an element. Now, we can examine the approximate truncation error by running simulations ($g4,p2$), ($g2,p6$), and ($g1,p7$), and evaluating the highest order non-zero derivatives of the approximate solution $\phi_h$. To evaluate the derivatives, we interpolate the solution onto an orthogonal modal polynomial basis, that is we find the coefficients $a_{ij}$ such that $\phi_h = \sum_i \phi_i\theta_i = \sum_{ij} a_{ij} P_{ij}$, where $P_{ij}$ is a modal orthogonal polynomial with maximum degree of $i$ on $x$ and $j$ on $y$, for a maximum total degree of $i+j$. The derivatives then evaluate as $\frac{1}{2!} \left(\der{}{x}+\der{}{y} \right)^2 \phi_h(\mathbf x_\eta) = \sum_{i+j=2} a_{ij}$, $ \frac{1}{6!}\left(\der{}{x}+\der{}{y} \right)^6 \phi_h(\mathbf x_\xi)= \sum_{i+j=6} a_{ij}$, and $ \frac{1}{7!}\left(\der{}{x}+\der{}{y} \right)^7 \phi_h(\mathbf x_\xi)= \sum_{i+j=7} a_{ij}$, that is, we simply need to sum the coefficients of the modal orthogonal polynomial basis which correspond to terms with total degree of 2, 6, and 7, respectively. Since the coefficients of the polynomials are evaluated on the reference element, $h \approx 1$ will be the same for all elements. Also, while this approach gives an estimate of the leading order truncated term, it does not give an exact value. In our case, we are not interested in a rigorous error estimator, but instead we only require an estimate of the error to aid the discussion. 

Our approach is similar to that followed by \citet{mavriplis89}, where Legendre polynomials were used instead. The author proposed a smoothness estimator, where the coefficients $a_{ij}$ are fit to the exponentially decaying function $a(i+j) = Ce^{\sigma (i+j)}$. There the author claims that $\sigma<-1$ indicates good resolution or smooth functions, and $\sigma>-1$ indicates poor resolution or non-smooth functions. The adaptive strategy used was to increase the polynomial degree for elements with $\sigma<-1$, and to refine the mesh for elements with $\sigma>-1$, if the error level in that element was insufficient. We evaluate this smoothness indicator $\sigma$ on $(g1,p7)$ by doing a least squares fit of the coefficients to $Ce^{\sigma (i+j)}$. In regions where the magnitude of the solution is close to 0, that is below the euphotic zone for the Zooplankton field, $a_{ij}\; \forall i,j$ will be small, and the smoothness indicator $\sigma$ will not be accurate. The approximate size of the truncated derivative terms along with the smoothness indicator are plotted in \fref{truncerrs}. Only the smoothness indicators calculated on $(g2,p5)$ and $(g1,p7)$ are plotted since the accuracy of the smoothness indicator improves with the number of terms in the polynomial expansion, and is not accurately represented on $(g4,p2)$ \citep{mavriplis89}.

\begin{figure*}%
a)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig18a-BIO_NT10_g4_P2_truncder}\end{center}
b)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig18b-BIO_NT10_g2_P6_truncder}\end{center}
c)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig18c-BIO_NT10_g1_P7_truncder}\end{center}
d)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig18d-BIO_NT10_g2_P6_truncsmooth}\end{center}%
e)\begin{center}\vspace{-12pt} \includegraphics[width=0.95\textwidth]{fig18e-BIO_NT10_g1_P7_truncsmooth}\end{center}%
\caption{Approximate truncation errors for Zooplankton fields at $t^*=20$ (with $\bar \tau_a=12.5$[days]). Calculated on a) ($g4,p2$) using $log_{10}(\sum_{i+j=2} a_{ij})$, on b) $g2,p6$ using $log_{10}(\sum_{i+j=6} a_{ij})$, and on c) $g1,p7$ using $log_{10}(\sum_{i+j=7} a_{ij})$. d-e) Smoothness indicator $\sigma$ calculated on d) $(g2,p6)$ and e) $(g1,p7)$ . The biological dynamics used has stable limit cycles within the euphotic zone (bio case 3).}%
\label{fig:truncerrs}%
\end{figure*}

From \fref{truncerrs}, we note that the largest differences in \fref{npzdiff} correspond to the regions with the largest approximate truncations errors in \fref{truncerrs}. Also, in the region $x^*>7$ where the low-order solution is more accurate than \HO, we have $\sigma > -1$, which suggests that refining the elements instead of the order of accuracy is more appropriate. After one level of refinement on $(g2,p5)$, we see that the smoothness indicator shows a smaller region of non-smooth elements. This illustrates that the smoothness is defined in terms of the numerical discretization, and is not solely a function of the solution field. Also note in the region where the high-order solution is more accurate (see \fref{npzdiffdet}), the approximate derivative of the truncation term is small in both fields, and $\sigma < -1$, suggesting a higher degree polynomial basis is more appropriate in this region. This shows that our mesh is not optimized in terms of the solution field, and highlights the importance of using both mesh refinement and polynomial basis adaptation to generate an optimal discretization for complex biological ocean dynamics. Also, this shows that whether a coarsely discretized higher order scheme is better than a finely discretized lower order scheme depends on the smoothness of the solution and can vary spatially across the solution. The benefit from a higher-order solution is as follows: when the solution is smooth, increasing the polynomial degree causes the error to decrease exponentially, whereas the error would only decrease algebraically if decreasing the element size. The cost of increasing the polynomial degree also scales algebraically, and because of this, a higher-order scheme performs better for smooth or well-resolved fields. Using our implementation, the $(g1,p6)$ simulation took approximately 0.34 of the time taken by the $(g4,p1)$ simulation. We also ran $(g3,p1)$, which was approximately three times more efficient than $(g1,p6)$, but this solution (not shown) was less than 1\% accurate for the majority of the domain. When the solution is not resolved (i.e.~not smooth for the grid resolution or polynomial degree chosen), higher-order schemes will lead to Gibbs oscillations and filtering is required (see \ref{sec:filter}), while lower-order schemes may "look good" but will be very dissipative. When the solution is resolved (i.e.~smooth enough for the grid resolution or polynomial degree chosen), higher-order discretizations perform better than lower-order ones: they are more accurate and less dissipative for the same cost.  

Finally, we note that the approximate truncation error and smoothness metrics were different for the different biological components. Therefore, the optimal discretization for one component is not the same as the optimal discretization for another component. \citet{ueckermann09} proposed a scheme that uses a different order basis function for different biological components, but also cautions that an incurred interpolation cost needs to be considered for adaptation strategies. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%% EVOLUTION OF BIOLOGICAL PATCH %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
\subsection{Evolution of Biological Patch}\label{sec:patch}
In this section we demonstrate how biological activity can enhance the differences between low order and high order discretizations beyond the effect of numerical dissipation alone. For this example, we modify bio case 1 (single stable points at all depths) from \sref{bioadvect} by introducing a vertical column, or ``patch'', of biology that uses the parameters from bio case 2 (stable limit cycles at depths $z^*	=0.4 -0.9$ and single stable points elsewhere). This is easily done in the dimensional form of the equations by increasing the value of $\mathcal N_T$ locally in the patch. Such situations occur frequently in nature, e.g.~ an eddy or front upwelling additional nutrients locally towards the surface. The initial condition and boundary condition is the same as in bio case 1 (the steady state solution with smoothed discontinuity), except inside the patch where the initial conditions for bio case 2 are used instead, that is: 
\begin{eqnarray}
	\phi^*_{(\mathrm{patch})} &=& \phi^*_{(\mathrm{bio~case~1})} + \left [\phi^*_{(\mathrm{bio~case~2})}-\phi^*_{(\mathrm{bio~case~1})}\right ] \cdot e^{-\frac{(x^*+6.4)^4}{2 \cdot (8^4)}},
\label{eq:patchinit}
\end{eqnarray}
where $\phi^*_{(\mathrm{patch})}$ is the initial condition used for this example, $\phi^*_{(\mathrm{bio~case~1})}$ is the steady state with smoothed discontinuity for bio case 1, and $\phi^*_{(\mathrm{bio~case~2})}$ is the steady state with smoothed discontinuity for bio case 2. Note, for this example, we non-dimensionalize $\phi^* = \frac{\phi}{\mathcal N_{T,\mathrm{bio~case~1}}}$ by the total biomass for bio case 1. In addition we superimpose a periodic velocity onto the mean velocity, 
\begin{eqnarray}
	\mathbf u^* = \mathbf u^*_{\mathrm{mean}}  [1+5\cdot \mathrm{sign} \{\cos(0.16 \pi t^*)\}],
\label{eq:patchu}
\end{eqnarray}
where $\mathbf u^*$ is now the velocity used for this example, and $\mathbf u^*_{\mathrm{mean}}$ is the potential flowfield solved from \sref{flow}. The superimposed velocity increases the distance traveled, as well as the number of time integration steps (due to the CFL condition), and therefore has the effect of increasing the numerical dissipation. 

\begin{figure*}%
\includegraphics[width=0.95\columnwidth]{fig19-Blob_Detail_grp}
\caption{Detail around the biological patch with stable limit cycles at the bottom of the euphotic zone at time $t^* = 14.4$ for a) the phytoplankton fields and b) the total biomass.  The solution for $(g2,p5)$ is plotted on the left, $(g4,p1)$ in the middle, and the difference between the solutions, $[(g2,p5)-(g4,p1)]$, is plotted on the right. This shows that $(g2,p5)$ correctly maintains the full peak of the biological patch, while $(g4,p1)$ does not, leading to large differences in the phytoplankton fields.}%
\label{fig:patch}%
\end{figure*}

\begin{figure*}%
\includegraphics[width=0.95\columnwidth]{fig20-Blob_Diff_Time}
\caption{The relative normed difference between the total biomass of the two solutions, ($Q_1$, \eref{patch1}), the sum of relative normed differences between the biological components ($Q_2$ \eref{patch2}), the relative normed difference in production ($Q_3$ \eref{patch3}), and the relative normed difference in grazing ($Q_4$ \eref{patch4}) over time from $t^*=0$ to $t^*=14.4$. This shows that the difference in biological components is amplified beyond the effect of numerical dissipation due to differences in the source terms such as the production and grazing.}%
\label{fig:patch2}%
\end{figure*}


\fref{patch} plots the phytoplankton fields and total biomass for $(g2,p5)$ and $(g4,p1)$ around the patch, as well as the difference of the solution, $(g2,p5)$ minus $(g4,p1)$, at $t^* =14.4$. We do not use the $(g5,p1)$ solution (as was done in \sref{bioadvect}) for the difference plots here because we found that even ($g5,p1$) is more dissipated than $(g2,p5)$, and therefore $(g2,p5)$ is more accurate inside the biological patch where our calculations are performed.  These results show that the total biomass peak is not maintained by the low order scheme, $(g4,p1)$, and the details in the phytoplankton fields are also dissipated. Since these simulations do not contain physical diffusion, any diffusion is due to the numerical scheme, and therefore the $(g2,p5)$ solution is more accurate than $(g4,p1)$ because it does maintain the total biomass peak. Apart from the effects of the periodic velocity, the solution inside the patch should resemble that of \fref{npzfnl}b), and $(g2,p5)$ resembles this solution more closely than $(g4,p1)$.

While some of the differences between the $(g2,p5)$ and $(g4,p1)$ simulations can be accredited solely to the numerical dissipation, the error due to numerical dissipation is amplified by the change in biological activity. To illustrate this point, we show in \fref{patch2}: the relative normed difference between the total biomass of the two solutions, ($Q_1$, \eref{patch1}), the sum of relative normed differences between the biological components ($Q_2$ \eref{patch2}), the relative normed difference in production ($Q_3$ \eref{patch3}), and the relative normed difference in grazing ($Q_4$ \eref{patch4})
\begin{eqnarray}
    Q_1 &=& \frac{\norm{\{\phi_N^*+\phi_P^*+\phi_Z^*\}_{(g2,p5)} - \{\phi_N^*+\phi_P^*+\phi_Z^*\}_{(g1,p4)}}^{\mathrm{patch}}}{\norm{\{\phi_N^*+\phi_P^*+\phi_Z^*-1\}}^{\mathrm{patch}}} \label{eq:patch1}\\
Q_2 &=& \frac{\sum_{I=(N,P,Z)}\norm{\phi^*_{I,(g2,p5)}- \phi^*_{I,(g1,p4)}}^{\mathrm{patch}}}{\norm{\{\phi_N^*+\phi_P^*+\phi_Z^*-1\}_{(g4,p1)}}^{\mathrm{patch}}} \label{eq:patch2}\\
    Q_3 &=&\frac{\norm{\left \{\mathcal U^*e^{z^*/h^*} \frac{\phi_P^* \phi_N^*}{\phi_N^* + k_s^*}\right\}_{g2p5} - \left \{ \mathcal U^*e^{z^*/h^*} \frac{\phi_P^* \phi_N^*}{\phi_N^* + k_s^*}\right \}_{g4p1} }^{\mathrm{patch}}}{\norm{\left \{\mathcal U^*e^{z^*/h^*} \frac{\phi_P^* \phi_N^*}{\phi_N^* + k_s^*}\right\}_{g4p1}}^{\mathrm{patch}}} \label{eq:patch3}\\
    Q_4 &=&\frac{\norm{\left\{a g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*})\right\}_{g2p5} - \left\{a g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*})\right\}_{g4p1} }^{\mathrm{patch}}}{\norm{\left\{a g_\nu^*\phi_Z^*(1-e^{-\nu^* \phi_P^*})\right\}_{g4p1}}^{\mathrm{patch}}}, \label{eq:patch4}
\end{eqnarray} 
where $\norm{e}^{\mathrm{patch}}= \left(\int_\mathrm{patch} e^2 dx^*dz^*\right)^{\frac{1}{2}}$ with the patch area is determined from $(g4,p1)$, and the quantity $\norm{\{\phi_N^*+\phi_P^*+\phi_Z^*-1\}}^{\mathrm{patch}}$ gives the size of the difference between the base solution and the solution inside the patch since the base number of nutrients (non-dimensionalized to 1) is subtracted out. 

Since our numerical scheme conserves the total biomass, the first quantity, $Q_1$, gives a quantitative estimate of the numerical dissipation error only. The second quantity, $Q_2$, should be the same as $Q_1$ if the only difference is due to numerical dissipation. However, from \fref{patch2}, we note that $Q_2 > Q_1$, which means the differences in dissipation is amplified by nonlinearities and the biology. 
This is explained by the differences in biological terms in the
two simulations, for example in the production and
grazing terms, $Q_3$ and $Q_4$. Also, note that these differences are growing over time, and for longer integration periods, the differences will be even greater. As a final note, the initial differences between the two solutions are due to interpolation errors, since the polynomial representation and number of degrees of freedom are not the same for the two simulations.

The example shows that the numerical dissipation due to a lower-order numerical scheme can be amplified by the biological reaction terms. 
This is significant since for accurate biological ocean
science through numerical simulations,
it is important to 
maintain the amplitudes of biological patches. This is particularly true for biology with multiple attractors, where relatively small perturbations can lead to vastly different solutions. The conclusion is that for the same cost,
higher-order schemes on coarser grids are more accurate
than lower-order schemes on finer grids.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%%%%%%%%%%%%%  FILTERING  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	

\subsection{Filtering Based on Smoothness Index} \label{sec:filter}
Based on the results in \sref{bioadvect}, and on the paper by \citet{hesthaven08b}, we created a selective exponential filter. We found that the default filter described in \citet{hesthaven08b}, which is applied at every time step, to result in a diffuse interface, even for high-order filters. This prompted us to develop the following filter. Consider a modal representation of the numerical solution, $\phi_h = \sum_{ij} a_{ij} P_{ij}$
(see \citet{ueckermann09}). After each time integration step, the solution on every element is modified as follows:
\begin{eqnarray}
\phi_h^F &=& \sum_{ij} \sigma(\eta_{ij}) a_{ij} P_{ij},
\label{eq:expfilt}
\end{eqnarray}
where
\begin{eqnarray}
 \sigma &=& \left \{
\begin{array}{cl}
	\mathrm{exp}(-\alpha \eta_{ij}^s), & \mathrm{if}~SI > 0 \\
	1, & \mathrm{if}~SI < 0 
\end{array} \right .\\
 SI &=& \frac{a_{ij}}{\max_{i+j \ge 1}(a_{ij})} - \exp(-(i+j)+1)\\
 \eta &=& \frac{i+j}{p+1} \\
 \alpha &=& - \log(0.01),
\end{eqnarray}
with $s$ being the order of the filter and $SI$ the smoothness index. 

This results in a filter that is applied only when the smoothness index indicates that the solution is not smooth. The smoothness index relies on information about the decay rate of the modal coefficients, so it cannot be used for lower-order schemes. Using this new filter, we obtain a solution where spurious numerical oscillations are no longer generated by the numerical scheme, and the interfaces are not significantly diffused. The new errors plots are shown in \fref{advectfilt}. We also verified that this filter does not affect smooth regions of the domain with other idealized test cases (not shown). 

While this filter looks promising, additional verification is needed, and a number of improvements are also possible. In particular, determining the strength or order, $s$, of the filter in an optimal manner is an open question. Also, selectively applying the filter only in directions where the solution is non-smooth is a topic of future research.

\begin{figure*}%
a) \begin{center} \includegraphics[width=0.95\textwidth]{fig21-ADVECT_filt_grp} \end{center}
\vspace{-12pt}\caption{Solution of three passive tracer cases as Fig.~9, for a) $(g2,p5)$, and b) $(g1,p6)$ but using our selective exponential filter \eref{expfilt} with $s=5$. Field values outside initially specified fields as Fig.~10, with c) $(g2,p5)$, and d) $(g1,p6)$ using the same filter. For each case a), b), c) and d), the filter damps the initialized oscillations of cases 1-2
and no oscillations are created in case 3.}%
\label{fig:advectfilt}%
\end{figure*}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%%%%%%%%%%%            CONCLUSIONS          %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}\label{sec:concl}
We completed a set of computational studies for the modeling of multi-scale biogeochemical dynamics in coastal ocean regions with complex bathymetric features, utilizing recent advances in computational fluid dynamics. Specifically, we compared low to high order discretization schemes, both in time and space, employing standard and hybrid discontinuous Galerkin Finite Element Methods, on both straight and curved elements. We studied the effects of a varied set of numerical properties including: quadrature-free and quadrature-based discretizations of the source terms; order of the spatial discretizations of advection and diffusion operators; order of the temporal discretization in explicit schemes; and, resolution of the spatial mesh, with and without our new curved elements. We verified the convergence of our numerical schemes for both the biology and flow fields, validated the codes on analytical solutions, and completed a rigorous truncation error analysis.

Our numerical analyses concentrated on the nonlinear nutrient-phytoplankton-zooplankton dynamics under advection and diffusion within an ocean strait or sill, in an idealized 2D geometry. We first non-dimensionalized the PDEs, evaluated stability regions and selected three biological dynamical regimes: single stable points at all depths, stable limit cycles at the bottom of the euphotic zone, and stable limit cycles within the whole euphotic zone (the latter two cases have limit cycles that are depth and light dependent). We evaluated the effects of numerical parameters on the three biological regimes, but illustrated only the most relevant results. In addition, for each of these biological regimes, we examined three types of coupled physics-biology interactions: biological terms dominating; advection terms dominating; and advection and biological terms balancing. For the advection-dominating case, we studied the advection over a strait of a uniform tracer, horizontal front and vertical front. For the balanced situation, relatively common in the real ocean, we considered biological dynamics that were either as fast as (e.g. coastal ocean) or slower than advection time-scales.


	In the regime where biological terms dominate, we found that both the quadrature-based and quadrature-free treatment of the source-terms give accurate, convergent results, although the quadrature-based algorithm had slightly smaller errors. We also showed that oscillations can occur solely due to numerics (Gibbs-like phenomena) for a high-order discretizations. A key result is that, for any numerical scheme, careful one-dimensional studies should be performed to understand the potential errors from the nonlinear source-term discretization.
	
	For the advection-dominating regime, we confirmed the flow field convergence and using passive tracers studied numerical advection artifacts that would also affect the biology. We found that for discretizations that do not resolve the solution, oscillations due to discontinuities in the tracers could be large for both low-order and high-order discretizations, but can be damped using a filtering approach for the high-order case. However, our results suggested that when the solution is resolved enough, higher-order schemes on coarser grids perform better (higher accuracy, less dissipative) for the same cost than lower-order scheme on finer grids, as long as curved boundaries were used.
 
 For the case of approximately balanced advection and biological terms, we compared low and high order temporal and spatial discretizations, and studied quadrature-based and quadrature-free discretizations of the source terms. We found that for lower-order temporal discretizations, the errors grew rapidly and would lead to inaccurate solutions for applications with faster biological timescales or longer integration times. We also showed that the quadrature-based source-term discretization was more accurate in regions where the solution was under-resolved, but in well-resolved regions, there was only a 0.1\% discrepancy, and the quadrature-free algorithm could be used for efficiency purposes. By quantitatively evaluating the truncation error and smoothness of the solution fields, we confirmed that higher-order spatial discretizations were more accurate in regions where the solution was smooth (i.e.~resolved enough), but less accurate where non-smooth (un-resolved) due to Gibbs-like oscillations. To reduce these oscillations, we developed a new numerical filter that is active only when and where the solution is not smooth locally, using a smoothness indicator. Finally, we demonstrated the importance of non-dissipative numerical schemes when biological patches are present which is common in the real ocean. First, we found that effects of numerical dissipation were amplified by biological activity, causing dissipation errors to increase faster with integration time. Higher-order spatial discretizations were more accurate when modeling biological patches because they maintained the patches while lower-order schemes did not. For resolved biology (e.g.~as in \fref{npzdiffstable}), higher-order schemes on coarser grids were for the same cost more accurate than lower-order schemes on finer grids. This conclusion is most important for longer-term simulations. It has major implications for fundamental studies of biological blooms, patchiness and other nonlinear dynamics in coastal regions with complex bathymetric features such as straits, sills, ridges and shelfbreaks. One can expect similar implications for longer-term eddy-resolving ecosystem studies or climate applications.
 

Based on our results, future research directions are to further develop schemes to reduce Gibbs-like oscillations without significant loss of accuracy and efficiency (e.g. \cite{persson06}). Without oscillation limiters or filtering, the optimal performance could be obtained by using different polynomial degree basis functions in the domain, where low-order elements could be used in non-smooth regions while high-order elements could be used in smooth regions. Because the smoothness can be determined from the discretization, an adaptive grid and polynomial degree scheme could be developed. Another possibility in this case could be to increase the grid resolution and decrease the order of schemes (e.g. medium order schemes, i.e. ($g3,p3$) or ($g3,p4$)) up to the point when numerical oscillations reach the size of other errors. Another research direction is to develop and evaluate schemes that would preserve the non-negativity of the biological solution. Our results can now be utilized for idealized studies of biological dynamics in straits or sills. Uncertainty quantifications \citep{lermusiaux06,sapsis09} as well as adaptive model learning \citep{lermusiaux07} for biological predictions would also be useful.  Finally, we are now well positioned to implement these new methods in three-dimensional ocean modeling systems (e.g. \citet{mseas}) for realistic coupled biogeochemical-physical ocean science and applications.



\begin{acknowledgements}
We are grateful to Dr. Patrick Haley Jr. and Ms. Lisa Burton for their collaboration on numerics and biology. We thank Dr. N.C. Nguyen for discussions on DG and HDG, and Dr. J. Peraire discussions on numerics. We thank the Office of Naval Research for sustained research support under grants N00014-07-1-1061, N00014-08-1-1097 and N00014-07-1-0473 to the Massachusetts Institute of Technology. We are grateful to the Natural Sciences and Engineering Research Council (NSERC) of Canada for the Postgraduate Scholarship to the first author, and to Mr. and Mrs. Pappalardo for partially supporting this research through a Pappalardo Fellowhip.  We thank Dr.~M.~Iskandarani and the reviewers for detailed and helpful comments and suggestions.  
\end{acknowledgements}

% BibTeX users please use one of
%\bibliographystyle{spbasic}      % basic style, author-year citations
%\bibliographystyle{spmpsci}      % mathematics and physical sciences
%\bibliographystyle{spphys}       % APS-like style for physics
%\bibliography{}   % name your BibTeX data base
\bibliographystyle{spbasic_pjh}

\bibliography{ODYN-D-10-00075R1}

\end{document}
% end of file NPZ_DG.tex
