# A pseudospectral quadrature method for Navier-Stokes equations on rotating
spheres

M. Ganesh  Department of Mathematical and Computer Sciences, Colorado School
of Mines, Golden, CO 80401<EMAIL_ADDRESS>,  Q. T. Le Gia  School of
Mathematics and Statistics, University of New South Wales, Sydney, NSW 2052,
Australia<EMAIL_ADDRESS>and  I. H. Sloan  School of Mathematics and
Statistics, University of New South Wales, Sydney, NSW 2052, Australia
<EMAIL_ADDRESS>
###### Abstract.

In this work, we describe, analyze, and implement a pseudospectral quadrature
method for a global computer modeling of the incompressible surface Navier-
Stokes equations on the rotating unit sphere. Our spectrally accurate
numerical error analysis is based on the Gevrey regularity of the solutions of
the Navier-Stokes equations on the sphere. The scheme is designed for
convenient application of fast evaluation techniques such as the fast Fourier
transform (FFT), and the implementation is based on a stable adaptive time
discretization.

###### Key words and phrases:

Navier-Stokes equations, unit sphere, vector spherical harmonics

###### 2000 Mathematics Subject Classification:

Primary 65M12; Secondary 76D05

##  1\. Introduction

In this paper we develop a pseudospectral quadrature method for the surface
Navier-Stokes partial differential equations (PDEs) on the rotating unit
sphere. Whereas the finite element method is best suited for handling non-
smooth processes, the spectral global basis computer models are very efficient
and perform extremely well for processes with smooth regularity. For example,
exponential convergence properties of the global Fourier basis spectral
Galerkin methods (without quadrature) for the Ginzburg-Landau and Navier-
Stokes PDEs on two dimensional periodic cells are based on the Gevrey
regularity of solutions of the PDEs [7, 20].

The complex three dimensional flows in the atmosphere and oceans are
considered to be accurately modeled by the Navier-Stokes PDEs of fluid
mechanics together with classical thermodynamics [22]. Difficulties in
computer modeling in these PDEs resulted in several simplified models for
which spectral approximations are well known [2, 15, 22]. A famous open
problem is to prove the global regularity for the three dimensional
incompressible Navier-Stokes PDEs [26]. However, the precise Gevrey regularity
of the unique solution of the (practically relevant) surface Navier-Stokes
PDEs on the rotating sphere was proved in [5]. (Because the Earth’s surface is
an approximate sphere, a standard surface model, to study global atmospheric
circulation on large planets, is the sphere.)

Consequently, a natural next step is to describe, analyze, and implement an
exponentially converging pseudospectral method for the Navier-Stokes PDEs on
the rotating sphere. In addition to the continuous model regularity results in
[5], this paper is also motivated by the recent work [12], where discrete
computer modeling of the Navier-Stokes PDEs on one-dimensional and toroidal
domains [8] was extended to the unit sphere.

This paper is concerned with both implementation of our algorithm and its
numerical analysis. The main numerical analysis contributions compared to
results in [5], for the continuous problem, and in [12], for a discrete
problem, are as follows.

For the spatially discrete pseudospectral quadrature Galerkin solutions of the
Navier-Stokes equations, we prove (i) the stability (that is, uniform
boundedness of approximate solutions, independent of their truncation
parameter $N$), see Theorem 4.1; and (ii) a spectrally accurate rate of
convergence [that is, $\mathcal{O}(N^{-2s})$ accuracy, with $s$ depending on
the smoothness of input data], see Theorem 4.2. We achieve these results by
first generalizing the main regularity result in [5] to complex valued times
(see Theorem 2.2 and 2.3), and then using this to prove the spectral rate of
convergence of the time derivative of a Stokes projection comparison function
(see Theorem 3.2). This time-derivative error result plays a crucial role in
the proof of spectral convergence of the approximate solutions (see the proof
of Theorem 4.2). We note that the main result in [12, page 978] establishes
only convergence of the semi-discrete Galerkin method without quadrature in
$L^{p}$ norms, but does not establish either stability or rate of convergence
of the scheme.

The rate of convergence results, supported by numerical experiments, formed
the core part of research on the Navier-Stokes equations on two dimensional
domains over the last few decades, see [8, 9, 27] and references therein.
There is a vast literature on numerical methods and analysis for the Navier-
Stokes PDE on bounded Euclidean domains (see [8, 9, 27] and references
therein), but their counterparts on closed manifolds are rarer (see [12] and
references therein). The implementation of the scheme in [12, page 978] is
based on a fixed time-step explicit Runge-Kutta method that has a small
stability region for the systems of ordinary differential equations arising
from the spatial discretization.

The outline of this paper is as follows. In the next section, we recall
various known preliminary results associated with the Navier-Stokes PDEs on
the unit sphere, in strong and weak form. In Section 3, we introduce essential
computational and numerical analysis tools required for the discretization and
analysis of the Navier-Stokes equations. In Section 4, we describe and prove
spectral accuracy of a pseudospectral quadrature method and give
implementation details required to apply the FFT and adaptive-in-time
simulation of the Navier-Stokes equations. In Section 5, we demonstrate
computationally the accuracy and applicability of the algorithm for well known
benchmark examples.

##  2\. Navier-Stokes equations on the rotating unit sphere

The surface Navier-Stokes equations (NSE) describing a tangential,
incompressible atmospheric stream on the rotating two-dimensional unit sphere
$S\subset\mathbb{R}^{3}$ can be written as [5, 16, 17, 19, 28]

(2.1) $\frac{\partial}{\partial t}{\bf u}+\boldsymbol{\mathrm{\nabla}}\,_{\bf
u}{\bf u}-\nu\boldsymbol{\mathrm{\Delta}}\,{\bf
u}+\boldsymbol{\mathrm{\omega}}\,\times{\bf
u}+\frac{1}{\rho}\boldsymbol{\mathrm{Grad}}\,p={\bf
f},\qquad\mathrm{Div}\,{\bf u}=0,\qquad{\bf u}|_{t=0}={\bf
u}_{0}\qquad\text{on}~{}~{}S.$

Here ${\bf u}={\bf u}(\widehat{{\bf x}},t)=\left(u_{1}(\widehat{{\bf
x}},t),u_{2}(\widehat{{\bf x}},t),u_{3}(\widehat{{\bf x}},t)\right)^{T}$ is
the unknown tangential divergence-free velocity field at $\widehat{{\bf x}}\in
S$ and $t\in[0,T]$, $p=p(\widehat{{\bf x}},t)$ is the unknown pressure. The
known components in (2.1) are the constant viscosity and density of the fluid,
respectively denoted by $\nu,\rho$, the normal vector field
$\boldsymbol{\mathrm{\omega}}\,=\boldsymbol{\mathrm{\omega}}\,(\widehat{{\bf
x}})=\omega(\widehat{{\bf x}})\widehat{{\bf x}}$ for the Coriolis acceleration
term, and the external flow driving vector field ${\bf f}={\bf
f}(\widehat{{\bf x}},t)$. The Coriolis function $\omega$ is given by
$\omega(\widehat{{\bf x}})=2\Omega\cos\theta$, where $\Omega$ is the angular
velocity of the rotating sphere, and $\theta$ is the angle between
$\widehat{{\bf x}}$ and the north pole. The vorticity of the flow associated
with the NSE (2.1), in the curvilinear coordinate system, is a normal vector
field, defined, for a fixed $t\geq 0$, by

(2.2) $\boldsymbol{\mathrm{Vort}}\,{\bf u}(\widehat{{\bf
x}},t)=\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf u}(\widehat{{\bf
x}},t)=\widehat{{\bf x}}\Delta\Psi(\widehat{{\bf
x}},t),\qquad\qquad\widehat{{\bf x}}\in S,$

for some scalar-valued vorticity stream function $\Psi$.

All spatial derivative operators in (2.1)-(2.2) are surface differential
operators, obtained by restricting the corresponding domain operators (defined
in a neighborhood of $S$) to the unit sphere, using standard differential
geometry concepts on closed manifolds in $\mathbb{R}^{3}$ [16, 17].

Using the fact that the outward unit normal at $\widehat{{\bf x}}\in S$ is
$\widehat{{\bf x}}$, the $\boldsymbol{\mathrm{Curl}}\,$ of a scalar function
$v$, of a normal vector field ${\bf w}=w\widehat{{\bf x}}$, and of a
tangential vector field ${\bf v}$ on $S$ are respectively defined by

(2.3) $\boldsymbol{\mathrm{Curl}}\,v=-\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Grad}}\,v,\qquad\boldsymbol{\mathrm{Curl}}\,{\bf
w}=-\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Grad}}\,w,\qquad\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}{\bf v}=-\widehat{{\bf x}}~{}\mathrm{Div}\,(\widehat{{\bf x}}\times{\bf
v}).$

The surface diffusion operator acting on tangential vector fields on $S$ is
denoted by $\boldsymbol{\mathrm{\Delta}}\,$ (known as the Laplace-Beltrami or
Laplace-de Rham operator) and is defined as

(2.4) $\boldsymbol{\mathrm{\Delta}}\,{\bf
v}=\boldsymbol{\mathrm{Grad}}\,\mathrm{Div}\,{\bf
v}-\boldsymbol{\mathrm{Curl}}\,\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}{\bf v}.$

The following relations connecting the above operators will be used throughout
the paper:

(2.5)
$\mathrm{Div}\,\boldsymbol{\mathrm{Curl}}\,v=0,\qquad\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}\boldsymbol{\mathrm{Curl}}\,v=-\widehat{{\bf x}}\Delta
v,\qquad\boldsymbol{\mathrm{\Delta}}\,\boldsymbol{\mathrm{Curl}}\,v=\boldsymbol{\mathrm{Curl}}\,\Delta
v,$ (2.6) $2\boldsymbol{\mathrm{\nabla}}\,_{\bf w}{\bf
v}=-\boldsymbol{\mathrm{Curl}}\,({\bf w}\times{\bf
v})+\boldsymbol{\mathrm{Grad}}\,({\bf w}\cdot{\bf v})-{\bf
v}\mathrm{Div}\,{\bf w}+{\bf w}\mathrm{Div}\,{\bf v}-{\bf
v}\times\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf w}-{\bf
w}\times\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf v}.$

In particular, for tangential divergence-free vector fields, such as the
solution ${\bf u}$ of the NSE, using (2.6), the nonlinear term in (2.1) can be
written as

(2.7) $\boldsymbol{\mathrm{\nabla}}\,_{\bf u}{\bf
u}=\boldsymbol{\mathrm{Grad}}\,\frac{|{\bf u}|^{2}}{2}-{\bf
u}\times\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf u}.$

###  2.1. A weak formulation

A standard technique for removing the scalar pressure field from the Navier-
Stokes equations is to multiply the first equation in (2.1) by test functions
${\bf v}$ from a space with elements having properties of the unknown velocity
field ${\bf u}$ (in particular, $\mathrm{Div}\,{\bf v}=0$) and then integrate
to obtain a weak formulation. (The unknown $p$, can then be computed by
solving a pressure Poisson equation, obtained by applying the surface
divergence operator in (2.1).)

To this end, we introduce the standard inner products on the space of all
square integrable (i) scalar functions on $S$, denoted by $L^{2}(S)$; and (ii)
tangential vector fields on $S$, denoted by $L^{2}(TS)$:

(2.8) $\displaystyle(v_{1},\;v_{2})$ $\displaystyle=$
$\displaystyle(v_{1},\;v_{2})_{L^{2}(S)}~{}~{}=\int_{S}v_{1}\overline{v_{2}}~{}dS,\quad\qquad
v_{2},v_{2}\in L^{2}(S),$ (2.9) $\displaystyle({\bf v}_{1},\;{\bf v}_{2})$
$\displaystyle=$ $\displaystyle({\bf v}_{1},\;{\bf
v}_{2})_{L^{2}(TS)}=\int_{S}{\bf v}_{1}\cdot\overline{{\bf
v}_{2}}~{}dS,\qquad{\bf u},{\bf v}\in L^{2}(TS),$

where $dS=\sin\theta d\theta d\phi$. Throughout the paper, the induced norm on
$L^{2}(TS)$ is denoted by $\|\cdot\|$ and for other inner product spaces, say
$X$ with inner product $(\cdot,\;\cdot)_{X}$, the associated norm is denoted
by $\|\cdot\|_{X}$. For example, for $s>0$, standard norms in the scalar and
vector valued functions Sobolev spaces $H^{s}(S)$ and $H^{s}(TS)$ are denoted
by $\|\cdot\|_{H^{s}(S)}$ and $\|\cdot\|_{H^{s}(TS)}$, respectively. Since
$H^{0}(TS)=L^{2}(TS)$, $\|\cdot\|_{H^{0}(TS)}=\|\cdot\|$.

We have the following identities for appropriate scalar and vector fields [16,
(2.4)-(2.6)]:

(2.10) $\displaystyle(\boldsymbol{\mathrm{Grad}}\,\psi,\;{\bf
v})=-(\psi,\;\mathrm{Div}\,{\bf v}),$
$\displaystyle\qquad(\boldsymbol{\mathrm{Curl}}\,\psi,\;{\bf
v})=(\psi,\;\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf v}),$ (2.11)
$\displaystyle(\boldsymbol{\mathrm{Curl}}\,\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}{\bf w},\;{\bf z})$ $\displaystyle=$
$\displaystyle(\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
w},\;\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf z}).$

In (2.10), the $L^{2}(TS)$ inner product is used on the left hand side and the
$L^{2}(S)$ inner product is used on the right hand side. Throughout the paper,
we identify a normal vector field ${\bf w}$ with a scalar field $w$ and hence

(2.12) $(\psi,\;{\bf w}):=(\psi,\;w)_{L^{2}(S)},\qquad{\bf w}=\widehat{{\bf
x}}w,\qquad\psi,w\in L^{2}(S).$

Using (2.5), smooth ($C^{\infty}$) tangential fields on $S$ can be decomposed
into two components, one in the space of all divergence-free fields and the
other through the Hodge decomposition theorem [1]:

(2.13) $C^{\infty}(TS)=C^{\infty}(TS;\boldsymbol{\mathrm{Grad}}\,)\oplus
C^{\infty}(TS;\boldsymbol{\mathrm{Curl}}\,),$

where

(2.14)
$C^{\infty}(TS;\boldsymbol{\mathrm{Grad}}\,)=\\{\boldsymbol{\mathrm{Grad}}\,\psi:\psi\in
C^{\infty}(S)\\},\;C^{\infty}(TS;\boldsymbol{\mathrm{Curl}}\,)=\\{\boldsymbol{\mathrm{Curl}}\,\psi:\psi\in
C^{\infty}(S)\\}.$

For $s\geq 0$, let $C^{\infty,s}(TS;\boldsymbol{\mathrm{Curl}}\,)$ denote the
closure of $C^{\infty}(TS;\boldsymbol{\mathrm{Curl}}\,)$ in the $H^{s}(TS)$
norm. In particular, following [5] we introduce a simpler notation

$\displaystyle H$ $\displaystyle=$ $\displaystyle\mbox{ closure of
}C^{\infty}(TS;\boldsymbol{\mathrm{Curl}}\,)\mbox{ in
}L^{2}(TS)\;=\;C^{\infty,0}(TS;\boldsymbol{\mathrm{Curl}}\,),$ $\displaystyle
V$ $\displaystyle=$ $\displaystyle\mbox{ closure of
}C^{\infty}(TS;\boldsymbol{\mathrm{Curl}}\,)\mbox{ in
}H^{1}(TS)\;=\;C^{\infty,1}(TS;\boldsymbol{\mathrm{Curl}}\,).$

Using the Gauss surface divergence theorem, for any scalar valued function $v$
on $S$ with $\boldsymbol{\mathrm{Grad}}\,v\in L^{2}(TS)$, using (2.10), we
have

(2.15) $(\boldsymbol{\mathrm{Grad}}\,v,\;{\bf
w})=\int_{S}\boldsymbol{\mathrm{Grad}}\,v\cdot\overline{{\bf
w}}~{}dS=-\int_{S}v\cdot\mathrm{Div}\,\overline{{\bf w}}~{}dS=0,\quad{\bf
w}\in V,$

and hence the unknown pressure can be eliminated from the first equation in
(2.1) through the weak formulation.

Following [16, Page 567], for the diffusion part of the NSE, we consider the
Stokes operator

(2.16)
$\boldsymbol{\mathrm{A}}=\boldsymbol{\mathrm{Curl}}\,\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}.$

Using (2.4) and (2.5), it is easy to see that the Stokes operator is the
restriction of the vector Laplace-de Rham operator
$-\boldsymbol{\mathrm{\Delta}}\,$ on $V$;
$\boldsymbol{\mathrm{A}}=-{\mathbf{P}}_{\boldsymbol{\mathrm{Curl}}\,}\boldsymbol{\mathrm{\Delta}}\,$,
where
${\mathbf{P}}_{\boldsymbol{\mathrm{Curl}}\,}:~{}L^{2}(TS)~{}\rightarrow~{}H$
is the orthogonal projection onto the divergence-free tangent space.

For each positive integer $L=1,2,\ldots$, the eigenvalue $\lambda_{L}$ and the
corresponding eigenvectors of the Stokes operator $\boldsymbol{\mathrm{A}}$
are given by

(2.17)
$\lambda_{L}=L(L+1),\qquad{\mathbf{Z}}_{L,m}(\theta,\varphi)=\lambda^{-1/2}_{L}\boldsymbol{\mathrm{Curl}}\,Y_{L,m}(\theta,\varphi),\quad
m=-L,\ldots,L,$

where $Y_{L,m}$ are the scalar orthonormal spherical harmonics of degree $L$,
defined by

(2.18)
$Y_{L,m}(\theta,\varphi)=\left[\frac{(2L+1)}{4\pi}\frac{(L-|m|)!}{(L+|m|)!}\right]^{1/2}P^{m}_{L}(\cos\theta)e^{im\varphi},\quad
m=-L,\ldots,L,$

with $P^{m}_{L}$ being the associated Legendre polynomials so that
$\overline{Y_{L,m}}=(-1)^{m}Y_{L,-m}$.

The spectral property
$\boldsymbol{\mathrm{A}}{\mathbf{Z}}_{L,m}=\lambda_{L}{\mathbf{Z}}_{L,m}$
follows from the fact that $Y_{L,m}$ are eigenfunctions of the scalar Laplace-
Beltrami operator $-\Delta$ with eigenvalues $\lambda_{L}$, the definition of
the Stokes operator $\boldsymbol{\mathrm{A}}$ in (2.16), ${\mathbf{Z}}_{L,m}$
in (2.17), (2.5), and (2.3). Since $\\{Y_{L,m}:L=0,1,\ldots;m=-L,\ldots,L\\}$
is an orthonormal basis for $L^{2}(S)$, it is easy to see that
$\\{{\mathbf{Z}}_{L,m}:L=1,\ldots;m=-L,\ldots,L\\}$ is an orthonormal basis
for $H$. Thus an arbitrary ${\bf v}\in H$ can be written as

(2.19) ${\bf v}=\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\widehat{{\bf
v}}_{L,m}{\mathbf{Z}}_{L,m},\qquad\widehat{{\bf v}}_{L,m}=\int_{S}{\bf
v}\cdot\overline{{\mathbf{Z}}_{L,m}}dS=({\bf v},\;{\mathbf{Z}}_{L,m}).$

We consider a subset of $H$,

(2.20) $\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2})=\left\\{{\bf v}\in H\;:{\bf
v}=\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\widehat{{\bf
v}}_{L,m}{\mathbf{Z}}_{L,m},\quad\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\lambda^{s}_{L}|\widehat{{\bf
v}}_{L,m}|^{2}<\infty\right\\},$

which is the divergence-free subset of the Sobolev space $H^{s}(TS)$. For
every ${\bf v}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2})$, we set

(2.21) $\|{\bf
v}\|_{H^{s}(TS)}=\left[\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\lambda^{s}_{L}|\widehat{{\bf
v}}_{L,m}|^{2}\right]^{1/2},$

and for ${\bf v}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2})$, we define

(2.22) $\boldsymbol{\mathrm{A}}^{s/2}{\bf
v}:=\quad\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\lambda^{s/2}_{L}\widehat{{\bf
v}}_{L,m}{\mathbf{Z}}_{L,m}\quad\in H.$

For a tangential vector field ${\bf v}$ on $S$, we define the Coriolis
operator $\boldsymbol{\mathrm{C}}$,

(2.23) $(\boldsymbol{\mathrm{C}}{\bf v})(\widehat{{\bf
x}})=\boldsymbol{\mathrm{\omega}}\,(\widehat{{\bf x}})\times{\bf
v}(\widehat{{\bf x}})=\omega(\widehat{{\bf x}})(\widehat{{\bf x}}\times{\bf
v}),\qquad\omega(\widehat{{\bf x}})=2\Omega\cos\theta.$

To treat the nonlinear term in (2.1), we consider the trilinear form $b$ on
$V\times V\times V$, defined as

(2.24) $b({\bf v},{\bf w},{\bf z})=(\boldsymbol{\mathrm{\nabla}}\,_{\bf v}{\bf
w},{\bf z})=\int_{S}\boldsymbol{\mathrm{\nabla}}\,_{\bf v}{\bf
w}\cdot{\overline{\bf z}}~{}dS,\qquad{\bf v},{\bf w},{\bf z}\in V.$

Using (2.6) and (2.10), for divergence free fields ${\bf v},{\bf w},{\bf z}$,
the trilinear form can be written as

(2.25) $b({\bf v},{\bf w},{\bf z})=\frac{1}{2}\int_{S}\left[-{\bf v}\times{\bf
w}\cdot\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\overline{\bf
z}}+\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf v}\times{\bf
w}\cdot{\overline{\bf z}}-{\bf
v}\times\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
w}\cdot{\overline{\bf z}}\right]~{}dS.$

Moreover [16, Lemma 2.1]

(2.26) $b({\bf v},{\bf w},{\bf w})=0,\qquad b({\bf v},{\bf z},{\bf w})=-b({\bf
v},{\bf w},{\bf z})\qquad{\bf v},{\bf w},{\bf z}\in V.$

Throughout the paper, the space $V$ is equipped with the norm
$\|\cdot\|^{2}_{V}=\left(\boldsymbol{\mathrm{A}}\cdot,\cdot\right)$.

Thus, using (2.4), (2.10), (2.16), and (2.25), a _weak solution_ of the
Navier-Stokes equations (2.1) is a vector field ${\bf u}\in L^{2}([0,T];V)$
with ${\bf u}(0)={\bf u}_{0}$ that satisfies the weak form

(2.27) $({\bf u}_{t},{\bf v})+b({\bf u},{\bf u},{\bf
v})+\nu(\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
u},\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
v})+(\boldsymbol{\mathrm{C}}{\bf u},{\bf v})=({\bf f},{\bf v}),\qquad{\bf
v}\in V.$

This weak formulation can be written in operator equation form on $V^{*}$, the
adjoint of $V$: Let ${\bf f}\in L^{2}([0,T];V^{*})$ and ${\bf u}_{0}\in H$.
Find a vector field ${\bf u}\in L^{2}([0,T];V)$, with ${\bf u}_{t}\in
L^{2}([0,T];V^{*})$ such that

(2.28) ${\bf u}_{t}+\nu\boldsymbol{\mathrm{A}}{\bf
u}+\boldsymbol{\mathrm{B}}({\bf u},{\bf u})+\boldsymbol{\mathrm{C}}{\bf
u}={\bf f},\qquad{\bf u}(0)={\bf u}_{0},$

where the bilinear form $\boldsymbol{\mathrm{B}}({\bf u},{\bf v})\in V^{*}$ is
defined by

(2.29) $(\boldsymbol{\mathrm{B}}({\bf u},{\bf v}),{\bf w})=b({\bf u},{\bf
v},{\bf w})\qquad{\bf w}\in V.$

In the subsequent error analysis, we need the following estimate for the
nonlinear term (see Lemma 6.1 in Appendix):

(2.30) $\|\boldsymbol{\mathrm{A}}^{-\delta}\boldsymbol{\mathrm{B}}({\bf
u},{\bf v})\|\leq\begin{cases}C\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf
u}\|\|{\bf v}\|\leq C\|\boldsymbol{\mathrm{A}}^{1/2}{\bf u}\|\|{\bf v}\|,&\cr
C\|{\bf u}\|\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf v}\|\leq C\|{\bf
u}\|\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
v}\|,&\end{cases}\quad\delta\in(1/2,1)\quad{\bf u},{\bf v}\in V.$

In (2.30), as throughout the paper, $C$ is a generic constant independent of
${\bf u}$ and ${\bf v}$, (and the discretization parameter $N$ introduced in
Section 3). From (2.30) we deduce the weak Lipschitz continuity property

(2.31) $\|\boldsymbol{\mathrm{A}}^{-\delta}(\boldsymbol{\mathrm{B}}({\bf
v},{\bf v})-\boldsymbol{\mathrm{B}}({\bf w},{\bf w}))\|\leq C\|{\bf v}-{\bf
w}\|,\quad\delta\in(1/2,1),\quad\text{if }\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
v}\|,\|\boldsymbol{\mathrm{A}}^{1/2}{\bf w}\|<C.$

The existence and uniqueness of the solution ${\bf u}\in L^{2}([0,T];V)$ of
the weak formulation (2.27) are discussed in [16, 17, 19]. A regular solution
of the Navier-Stokes equations (2.1) on $[0,T]$ is a tangential divergence-
free velocity field ${\bf u}$ that satisfies the equation obtained by
integrating in time the weak form (2.27), from $t_{0}$ to $t$, for almost
every $t_{0},t\in[0,T]$. In order to recall the existence, uniqueness, and
Gevrey regularity of the regular solution, we need a few more additional
details from [5]. These are also needed as tools for analyzing our
pseudospectral quadrature method.

###  2.2. Gevrey regularity of regular solution

The Gevrey class of functions of order $s>0$ and index $\sigma>0$, associated
with the Stokes operator defined in (2.16), is denoted by $G^{s/2}_{\sigma}$
and is defined as

(2.32) $G^{s/2}_{\sigma}:=\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2}e^{\sigma
A^{1/2}})\subset\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2}).$

Using (2.20), the Gevrey space

(2.33) $G^{s/2}_{\sigma}=\left\\{{\bf
v}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s/2})\ :{\bf
v}=\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\widehat{{\bf
v}}_{L,m}{\mathbf{Z}}_{L,m},\;\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\lambda^{s}_{L}e^{2\sigma\lambda^{1/2}_{L}}|\widehat{{\bf
v}}_{L,m}|^{2}<\infty\right\\}$

is a Hilbert space with respect to the inner product

(2.34) ${\left\langle{{\bf v}},{{\bf
w}}\right\rangle}_{G^{s/2}_{\sigma}}=\sum_{L=1}^{\infty}\sum_{m=-L}^{L}\lambda^{s}_{L}e^{2\sigma\lambda^{1/2}_{L}}\widehat{{\bf
v}}_{L,m}\overline{\widehat{{\bf w}}_{L,m}},\qquad{\bf v},{\bf w}\in
G^{s/2}_{\sigma}.$

First we recall the following result from [5, 29].

######  Theorem 2.1.

If ${\bf u}_{0}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2})$ and ${\bf
f}\in
L^{\infty}((0,\infty);\mathcal{D}(\boldsymbol{\mathrm{A}}^{s}e^{\sigma_{1}\boldsymbol{\mathrm{A}}^{1/2}}))$,
for some $s,\sigma_{1}>0$, then for all $t>0$ there exists a $T^{*}>0$,
depending only on $\nu,{\bf f}$, and $\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}_{0}\|_{L^{2}(TS)}$, such that the NSE (2.1) on $S$ have a unique regular
solution ${\bf u}(\cdot,t)$ and ${\bf u}(\cdot,t)\in G^{s+1/2}_{\sigma(t)}$,
where $\sigma(t)=\min\\{t,T^{*},\sigma_{1}\\}$.

In addition, from the assumption and proof of Theorem 2.1 [5, page 355],

(2.35)
$\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\sigma(t)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(t)\|^{2}\leq M_{0},\qquad t>0,$

where $M_{0}$ depends on $\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}(0)\|$,
$\sup_{t\geq 0}\|\boldsymbol{\mathrm{A}}^{s}{\bf f}(\cdot,t)\|$ and $\nu$ but
not on $t$. The bound in (2.35) is useful for establishing the quality of
approximation of the Stokes projection of ${\bf u}$ in the next section (see
Theorem 3.1). It is also convenient to have a similar bound for the time
derivative of the Stokes projections with $t$ in (2.35) replaced with certain
complex times $\zeta\in\mathbb{C}$, to prove the power of approximation of the
time derivative of the Stokes projection (see Theorem 3.2). To this end, we
consider the NSE extended to complex times $\zeta$,

(2.36) $\frac{d{\bf u}}{d\zeta}+\nu\boldsymbol{\mathrm{A}}{\bf
u}+\boldsymbol{\mathrm{B}}({\bf u},{\bf u})+\boldsymbol{\mathrm{C}}{\bf
u}={\bf f},\quad\mathrm{Div}\,{\bf u}=0,\quad{\bf u}(0)={\bf
u}_{0},\qquad\zeta\in\mathbb{C},\qquad\text{on}~{}~{}S,$

with standard complexification (see [9]) of all the spaces and operators
introduced earlier. In the next theorem, we extend arguments used in [9] for
the solution of the NSE on the plane to the case of the sphere. The arguments
differ in an essential way only for the nonlinear term.

######  Theorem 2.2.

Let ${\bf u}_{0}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2})$ and ${\bf
f}\in
C([0,T];\mathcal{D}(\boldsymbol{\mathrm{A}}^{s}e^{\sigma_{1}\boldsymbol{\mathrm{A}}^{1/2}}))$,
with $s~{}\geq~{}1/4$. Let the domain $\mathbb{T}$ be defined by

$\mathbb{T}:=\\{\zeta=re^{i\theta}:0\leq r\leq T;\ |\theta|\leq\pi/4\\}.$

We assume further that ${\bf f}(\cdot,\zeta)$ is analytic for
$\zeta\in\mathbb{T}$ and that

(2.37)
$K:=\sup\\{\|\boldsymbol{\mathrm{A}}^{s}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f}(\cdot,\zeta)\|^{2}:\zeta=re^{i\theta}\in\mathbb{T}\\}<\infty.$

Then there exists $T^{**}>0$ such that

(2.38)
$\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(\zeta)\|^{2}\leq M_{1},\qquad\zeta\in\mathbb{T}\mbox{ and }|\zeta|\leq
T^{**},$

where $M_{1}$ depends on $\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}(0)\|$ and
hence ${\bf u}(\cdot,\zeta)\in G^{s+1/2}_{\psi(r\cos\theta)}$, where
$\psi(x)~{}:=~{}\min\\{x,T^{**},\sigma_{1}\\}$.

###### Proof.

Let $\zeta=re^{i\theta}$ with $r>0$ and $|\theta|\leq\pi/4$, and let

${\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta):=\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(\zeta).$

The definition of $\psi$ gives $\frac{d}{dx}\psi:=\psi^{\prime}(x)\leq 1$, and
$\psi\leq\sigma_{1}$, thus for fixed $\theta$ we find

(2.39) $\displaystyle\frac{d}{dr}{\bf u}_{\boldsymbol{\mathrm{A}}}(\zeta)$
$\displaystyle=$
$\displaystyle\psi^{\prime}(r\cos\theta)\cos\theta\boldsymbol{\mathrm{A}}^{s+1}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(\zeta)$
$\displaystyle+e^{i\theta}\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}\frac{d{\bf
u}}{d\zeta}(\zeta).$

Using (2.39) in

$\displaystyle\frac{1}{2}\frac{d}{dr}\left\|{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta)\right\|^{2}=\Re\left(\frac{d}{dr}{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta),{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta)\right),$

where $\Re(\cdot)$ denotes the real-part function, we get

(2.40) $\displaystyle\frac{1}{2}\frac{d}{dr}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta)\|^{2}$ $\displaystyle=$
$\displaystyle\psi^{\prime}(r\cos\theta)\cos\theta\;\Re\;(\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta),{\bf u}_{\boldsymbol{\mathrm{A}}}(\zeta))$
$\displaystyle+\Re\;e^{i\theta}(\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}\frac{d{\bf
u}}{d\zeta}(\zeta),{\bf u}_{\boldsymbol{\mathrm{A}}}(\zeta)).$

Using (2.36) for the last term in (2.40) together with the fact that

$(\boldsymbol{\mathrm{A}}^{s}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}\boldsymbol{\mathrm{C}}{\bf
u},\boldsymbol{\mathrm{A}}^{s+1}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u})=0$

(see [5, Lemma 1]), we find

(2.41) $\displaystyle\frac{1}{2}\frac{d}{dr}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta)\|^{2}+\nu\cos\theta\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta)\|^{2}$ $\displaystyle=$
$\displaystyle\psi^{\prime}(r\cos\theta)\cos\theta(\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}},{\bf
u}_{\boldsymbol{\mathrm{A}}})-\Re\;e^{i\theta}(\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}\boldsymbol{\mathrm{B}}({\bf
u},{\bf u}),{\bf u}_{\boldsymbol{\mathrm{A}}})$
$\displaystyle+\Re\;e^{i\theta}(\boldsymbol{\mathrm{A}}^{s}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f},\boldsymbol{\mathrm{A}}^{1/2}{\bf u}_{\boldsymbol{\mathrm{A}}}),$

where in (2.41) and below we write ${\bf u}_{\boldsymbol{\mathrm{A}}}={\bf
u}_{\boldsymbol{\mathrm{A}}}(\zeta),{\bf u}={\bf u}(\zeta),{\bf f}={\bf
f}(\zeta),\psi=\psi(r\cos\theta)$. From [5, Lemma 2], we have (with
$p=\max\\{2-s,7/4\\}=7/4$, since $s\geq 1/4$),

(2.42)
$|(\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}\boldsymbol{\mathrm{B}}({\bf
u},{\bf u}),{\bf u}_{\boldsymbol{\mathrm{A}}})|\leq C\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{3-p}\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{p}.$

Applying $\psi^{\prime}\leq 1$, the Cauchy-Schwarz inequality, (2.42) and
Young’s inequality ($ab\leq a^{q}/q+b^{q^{\prime}}/q^{\prime}$ with
$1/q+1/q^{\prime}=1$) with $q=2$ and $q=2/(2-p)$ in (2.41), we get

$\displaystyle\frac{1}{2}\frac{d}{dr}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}+\nu\cos\theta\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}$ $\displaystyle\leq$
$\displaystyle\cos\theta\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|\|{\bf u}_{\boldsymbol{\mathrm{A}}}\|+C\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{3-p}\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{p}+\|\boldsymbol{\mathrm{A}}^{s}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f}\|\|\boldsymbol{\mathrm{A}}^{1/2}{\bf u}_{\boldsymbol{\mathrm{A}}}\|$
$\displaystyle\leq$
$\displaystyle\frac{\nu\cos\theta}{4}\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}+\frac{\cos\theta}{\nu}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}$
$\displaystyle+C^{\frac{2}{2-p}}\frac{(2-p)}{p}\left(\frac{p}{\nu\cos\theta}\right)^{\frac{p}{2-p}}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{\frac{2(3-p)}{2-p}}+\frac{\nu\cos\theta}{2}\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}$
$\displaystyle+\frac{1}{\nu\cos\theta}\|\boldsymbol{\mathrm{A}}^{s}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f}\|^{2}+\frac{\nu\cos\theta}{4}\|\boldsymbol{\mathrm{A}}^{1/2}{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}.$

Therefore,

(2.43) $\displaystyle\frac{d}{dr}\|{\bf u}_{\boldsymbol{\mathrm{A}}}\|^{2}$
$\displaystyle\leq$ $\displaystyle\frac{2\cos\theta}{\nu}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2}+C\left(\frac{1}{\nu\cos\theta}\right)^{\frac{p}{2-p}}\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{\frac{2(3-p)}{2-p}}+\frac{2}{\nu\cos\theta}\|\boldsymbol{\mathrm{A}}^{s}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f}\|^{2}$ $\displaystyle\leq$ $\displaystyle\frac{2\cos\theta}{\nu}(1+\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2})+C\left(\frac{1}{\nu\cos\theta}\right)^{\frac{p}{2-p}}(1+\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2})^{\frac{3-p}{2-p}}+$
$\displaystyle\frac{2}{\nu\cos\theta}\|\boldsymbol{\mathrm{A}}^{s}e^{\psi\boldsymbol{\mathrm{A}}^{1/2}}{\bf
f}\|^{2}.$

Using $|\theta|\leq\pi/4$, and $\frac{3-p}{2-p}=5$, in (2.43), we get

(2.44) $\frac{d}{dr}\|{\bf u}_{\boldsymbol{\mathrm{A}}}\|^{2}\leq C(1+\|{\bf
u}_{\boldsymbol{\mathrm{A}}}\|^{2})^{5}+\frac{2\sqrt{2}}{\nu}K.$

With $|\theta|\leq\pi/4$ fixed, let

$y(r)=1+\|{\bf u}_{\boldsymbol{\mathrm{A}}}\|^{2}.$

Then on using (2.44) and $y\geq 1$ we obtain

$\frac{d}{dr}y\leq Cy^{5},$

and hence on integrating the inequality we find that

$y(r)\leq 2y(0)$

provided that

$0\leq
r\leq\frac{15}{64C}\frac{1}{(y(0))^{4}}=\frac{15}{64C}\frac{1}{(1+\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}(0)\|^{2})^{4}}.$

By setting

$T^{**}:=\frac{15}{64C}\frac{1}{(1+\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}(0)\|^{2})^{4}},\quad M_{1}=1+2\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}(0)\|^{2},$

we deduce that (2.38) holds for $0\leq r\leq T^{**}$. ∎

We can extend the bound (2.38) to a larger domain which contains the interval
$[0,T]$ using the following property of the NSE solution on the sphere [29]:

(2.45) $\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}(\cdot,t)\|\leq M_{2}\mbox{
for all }t\in[0,T],$

where the constant $M_{2}$ depends only on
$\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}_{0}\|$, $\sup_{0\leq t\leq
T}\|\boldsymbol{\mathrm{A}}^{s}{\bf f}(\cdot,t)\|$ and $\nu$ but not on $T$.

######  Theorem 2.3.

Suppose ${\bf u}_{0}$ and ${\bf f}$ satisfy all the conditions in the domain
$\mathbb{T}$ as in Theorem 2.2. Then

(2.46)
$\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(\zeta)\|^{2}\leq M_{3},\qquad\zeta\in\mathbb{T}\mbox{ and }|\mbox{Im
}\zeta|\leq T^{**}/\sqrt{2},$

where $M_{3}:=1+2M_{2}^{2}$, $T^{**}$ depends on $M_{2}$, and hence ${\bf
u}(\cdot,\zeta)\in G^{s+1/2}_{\psi(r\cos\theta)}$ for all $\zeta\in\mathbb{T}$
with $|\mbox{Im }\zeta|\leq T^{**}/\sqrt{2}$.

###### Proof.

We proceed as in the proof of Theorem 2.2 to obtain the ordinary differential
equation

$\frac{d}{dr}y\leq Cy^{5},$

where
$y(r)=1+\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(re^{i\theta})\|^{2}$. On integrating the inequality we find that

$y(r)\leq 2y(0)$

provided that

$0\leq r\leq\frac{15}{64C}\frac{1}{(1+\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}(0)\|^{2})^{4}}.$

We define

$T_{1}(\rho):=\frac{15}{64C}\frac{1}{(1+\rho^{2})^{4}},\quad\rho\geq 0.$

If $T_{1}(\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}(0)\|)\geq T$ we have
finished the proof. Otherwise, we let

$T^{**}=T_{1}(M_{2}),$

where $M_{2}$ is given in (2.45). For $\zeta=re^{i\theta}\in\mathbb{T}$,
$0\leq r\leq T^{**}$,

(2.47)
$\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(re^{i\theta})\|^{2}\leq 1+2\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf
u}(0)\|^{2}.$

Consequently, (2.46) holds for $0\leq r\leq T^{**}$ with
$M_{3}=1+2\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf u}(0)\|^{2}$.

Next we consider the case $\zeta=T^{**}+re^{i\theta}$, with $r\in[0,T^{**}]$.
We define, for $|\theta|\leq\pi/4$,

${\bf v}(re^{i\theta}):={\bf u}(T^{**}+re^{i\theta}),\qquad r\in[0,T^{**}].$

Using $\|\boldsymbol{\mathrm{A}}^{s+1/2}{\bf v}(0)\|\leq M_{2}$ we can apply
the previous arguments to obtain (2.46) (with ${\bf u}$ replaced with ${\bf
v}$) for $0\leq r\leq T^{**}$. We complete the proof and obtain the bound
(2.46) by repeating the last argument $n$ times, where $n=\lceil
T/T^{**}\rceil$. ∎

##  3\. Finite dimensional spaces and Stokes projections

Throughout the remainder of the paper, with $s$ and $\sigma_{1}$ as in Theorem
2.1 and Theorem 2.2, we assume that

(3.1) ${\bf u}_{0}\in\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2}),\qquad{\bf
f}\in
C([0,T];\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2}e^{\sigma_{1}\boldsymbol{\mathrm{A}}^{1/2}})),\qquad
s\geq 1/4.$

Natural finite dimensional spaces (depending on a parameter $N>0$) in which to
seek approximations to ${\bf u}(t)$ are

(3.2) $V_{N}:=\mbox{ span
}\\{{\mathbf{Z}}_{L,m}:L=1,\ldots,N;~{}~{}m=-L,\ldots,L\\}.$

The dimension of $V_{N}$ is $N^{2}+2N$. Let ${\mathbf{\Pi}}_{N}:H\rightarrow
V_{N}$ be the orthogonal projection with respect to the $L^{2}(TS)$ inner
product defined by

(3.3) ${\mathbf{\Pi}}_{N}({\bf v})=\sum_{L=1}^{N}\sum_{m=-L}^{L}\widehat{{\bf
v}}_{L,m}{\mathbf{Z}}_{L,m}.$

######  Lemma 3.1.

Let $\alpha>0$ be given. If ${\bf
v}\in{\mathcal{D}}(\boldsymbol{\mathrm{A}}^{\alpha/2})$ then

(3.4) $\|{\bf v}-{\mathbf{\Pi}}_{N}{\bf v}\|\leq N^{-\alpha}\|{\bf
v}\|_{H^{\alpha}(TS)}.$

###### Proof.

Using (2.17), (2.19), and (2.21) we get

$\displaystyle\|{\bf v}-{\mathbf{\Pi}}_{N}{\bf
v}\|^{2}=\sum_{L=N+1}^{\infty}\sum_{m=-L}^{L}|\widehat{{\bf v}}_{L,m}|^{2}$
$\displaystyle\leq
N^{-2\alpha}\sum_{L=N+1}^{\infty}\sum_{m=-L}^{L}\lambda^{\alpha}_{L}|\widehat{{\bf
v}}_{L,m}|^{2}$ $\displaystyle\leq N^{-2\alpha}\|{\bf
v}\|^{2}_{H^{\alpha}(TS)}.$

∎

In particular, using (3.1), we get

(3.5) $\|{\bf f}-{\mathbf{\Pi}}_{N}{\bf f}\|\leq N^{-(2s+1)}\|{\bf
f}\|_{H^{2s+1}(TS)},\qquad t\in[0,T].$

For computer implementation, the Fourier coefficients in (3.3) and all
Galerkin type integrals in computational schemes for the NSE need to be
approximated by cubature rules on the sphere, leading to a pseudospectral
method. To this end, for a continuous scalar field $\psi$ on $S$, we consider
a Gauss-rectangle quadrature sum $Q_{M}(\psi)$ with quadrature points
$\\{\widehat{\boldsymbol{\xi}}_{p,q}=\boldsymbol{p}(\theta_{p},\phi_{q})\\}$
and positive weights ${w_{p}}$ of the form

(3.6)
$Q_{M}(\psi):=\frac{2\pi}{M}\sum_{q=1}^{M}\sum_{p=1}^{M/2}w_{p}\psi(\widehat{\boldsymbol{\xi}}_{p,q})=\frac{2\pi}{M}\sum_{q=1}^{M}\sum_{p=1}^{M/2}w_{p}\psi(\theta_{p},\phi_{q}),$

where $M\geq 2$ is an even integer, $w_{p}$ and $\cos\theta_{p}$ for $\quad
p=1,\ldots,M/2$ are the Gauss-Legendre weights and nodes on $[-1,1]$ and
$\phi_{q}=2q\pi/M$, $q=1,\ldots,M$. The rule (3.6) is exact when $\psi$ is a
polynomial of degree $M-1$ on $S$, that is,

$Q_{M}\psi=\int_{S}\psi\;dS,\qquad\psi\in\mathcal{P}_{M-1}.$

Hence corresponding to (2.8) and (2.9), we define discrete inner products for
scalar and vector fields on the unit sphere as

(3.7) $(v_{1},\;v_{2})_{M}=Q_{M}(v_{1}v_{2}),\qquad\qquad({\bf v}_{1},\;{\bf
v}_{2})_{M}=Q_{M}({\bf v}_{1}\cdot{\bf v}_{2}).$

The choice of $M$ is very important; we choose $M$ such that all Galerkin
integrals with polynomial terms in our scheme are evaluated exactly. In
particular, with the unknown tangential divergence-free velocity field sought
in the polynomial space $V_{N}$, and knowing that the NSE nonlinearity is
quadratic, we choose $M$ such that

(3.8) $(\boldsymbol{\mathrm{B}}({\bf v},{\bf w}),{\bf
z})=(\boldsymbol{\mathrm{B}}({\bf v},{\bf w}),{\bf z})_{M},\qquad{\bf v},{\bf
w},{\bf z}\in V_{N}.$

This holds, for example, if $M=3N+2$. We define a computable counterpart of
(3.3), using ${\mathbf{L}}_{N}:H\cap C(TS)\rightarrow V_{N}$, a discrete
orthogonal projection with respect to the $M^{2}/2$ point discrete inner
product, as

(3.9) ${\mathbf{L}}_{N}({\bf v})=\sum_{L=1}^{N}\sum_{m=-L}^{L}\widehat{{\bf
v}}_{L,m,M}{\mathbf{Z}}_{L,m},\qquad\widehat{{\bf v}}_{L,m,M}=Q_{M}({\bf
v}\cdot\overline{{\mathbf{Z}}_{L,m}})=({\bf v},\;{\mathbf{Z}}_{L,m})_{M}.$

With $M$ chosen to satisfy (3.8), it is easy to see that

(3.10) ${\mathbf{L}}_{N}({\bf v})={\mathbf{\Pi}}_{N}({\bf v})={\bf
v},\qquad{\bf v}\in V_{N},$

and for ${\bf v}\in H\cap C(TS)$ and ${\bf w}\in H\cap C^{k}(TS)$,

(3.11) $\|{\mathbf{L}}_{N}({\bf v})\|\leq C\|{\bf
v}\|_{\infty},\qquad\qquad\|{\bf w}-{\mathbf{L}}_{N}({\bf w})\|\leq
CN^{-k}\|{\bf w}\|_{C^{k}(TS)},$

where the last two inequalities follow from simple arguments used in Theorem
13 and Lemma 14 of [24]. In particular, since
$\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2})\subset H\cap C^{2s}(TS)$, for an
integer $2s$, using (3.1)

(3.12) $\|{\bf f}-{\mathbf{L}}_{N}({\bf f})\|\leq CN^{-2s}\|{\bf
f}\|_{C^{2s}(TS)},\qquad t\in[0,T].$

Next we consider the Stokes projection in $V_{N}$ of the exact unique regular
solution ${\bf u}(t):={\bf u}(.,t)$ of (2.1). For each fixed $t$, the Stokes
projection $\widetilde{{\bf u}}_{N}\in V_{N}$ of ${\bf u}$ is defined by

(3.13) $(\boldsymbol{\mathrm{A}}\widetilde{{\bf u}}_{N},{\bf
v})=(\boldsymbol{\mathrm{A}}{\bf u},{\bf v}),\qquad{\bf v}\in V_{N}.$

Since
${\mathbf{\Pi}}_{N}\boldsymbol{\mathrm{A}}=\boldsymbol{\mathrm{A}}{\mathbf{\Pi}}_{N}$,
it follows that $\widetilde{{\bf u}}_{N}={\mathbf{\Pi}}_{N}({\bf u})$.
Following standard techniques in finite element analysis, the Stokes
projection of ${\bf u}$ plays an important role as a comparison function in
the main analysis in the next section.

######  Theorem 3.1.

Let ${\bf u}_{0}$ and ${\bf f}$ satisfy (3.1). Then

(3.14) $\|{\bf u}-\widetilde{{\bf u}}_{N}\|\leq
C\lambda^{-s-1/2}_{N+1}e^{-\sigma(t)\lambda^{1/2}_{N+1}}\leq
CN^{-2s-1}e^{-\sigma(t)N},\qquad t\in[0,T],$

where $\sigma(t)$ is as in Theorem 2.1.

###### Proof.

Using the fact that $\widetilde{{\bf u}}_{N}={\mathbf{\Pi}}_{N}{\bf u}$, we
have

$\displaystyle\|{\bf u}-\widetilde{{\bf u}}_{N}\|^{2}$ $\displaystyle=$
$\displaystyle\sum_{L>N}\sum_{|m|\leq L}|\widehat{{\bf u}}_{L,m}|^{2}$
$\displaystyle\leq$
$\displaystyle\lambda^{-2s-1}_{N+1}e^{-2\sigma(t)\lambda^{1/2}_{N+1}}\sum_{L>N}\sum_{|m|\leq
L}\lambda^{2s+1}_{L}e^{2\sigma(t)\lambda^{1/2}_{L}}|\widehat{{\bf
u}}_{L,m}|^{2}$ $\displaystyle\leq$
$\displaystyle\lambda^{-2s-1}_{N+1}e^{-2\sigma(t)\lambda^{1/2}_{N+1}}\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\sigma(t)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(t)\|^{2}$ $\displaystyle\leq$ $\displaystyle
C\lambda^{-2s-1}_{N+1}e^{-2\sigma(t)\lambda^{1/2}_{N+1}}$

where in the last step we used (2.35). The last inequality in (3.14) follows
from the fact that $N^{2}\leq\lambda_{N+1}=(N+1)(N+2)$. ∎

######  Theorem 3.2.

Let ${\bf u}_{0}$, ${\bf f}$ satisfy (3.1). We assume further that ${\bf f}$
is analytic in $\mathbb{T}$ and (2.37) holds. Then for $t\in(0,T)$,

(3.15) $\left\|\frac{d}{dt}\left({\bf u}-\widetilde{{\bf
u}}_{N}\right)\right\|\leq
C\lambda^{-s-1/2}_{N+1}e^{-\psi_{1}(t)\lambda^{1/2}_{N+1}}\leq
CN^{-2s-1}e^{-\psi_{1}(t)N},$

where $\psi_{1}(t)=\min\\{(1-1/\sqrt{2})t,T^{**},\sigma_{1}\\}$, and $T^{**}$
is as in Theorem 2.3.

###### Proof.

Let $t\in(0,T)$ be fixed. Let ${\bf p}_{N}(\zeta)=({\bf u}-\widetilde{{\bf
u}}_{N})(\zeta)$ be the standard complexification of ${\bf u}-\widetilde{{\bf
u}}_{N}$ at $\zeta=re^{i\theta}$. Using Theorem 2.3 and the Cauchy integral
formula,

$\frac{d{\bf p}_{N}(t)}{dt}=\frac{1}{2\pi i}\int_{\Gamma}\frac{{\bf
p}_{N}(\zeta)}{(t-\zeta)^{2}}d\zeta,$

where for $t>0$, $\Gamma$ is a circle in the $\zeta$ plane with center $(t,0)$
and radius
$\min\\{t/\sqrt{2},T^{**}/\sqrt{2},T-t\\}$, a condition that ensures that
$\zeta=re^{i\theta}\in\Gamma$ lies in the region $\mathbb{T}$ with $|\mbox{Im
}\zeta|\leq T^{**}/\sqrt{2}$. Using the fact that $\widetilde{{\bf
u}}_{N}={\mathbf{\Pi}}_{N}{\bf u}$, for $\zeta=re^{i\theta}\in\Gamma$ we have,
from Theorem 2.3,

(3.16) $\displaystyle\|{\bf p}_{N}(\zeta)\|$ $\displaystyle=$
$\displaystyle\|{\bf u}(\zeta)-\widetilde{{\bf u}}_{N}(\zeta)\|^{2}$
$\displaystyle=$ $\displaystyle\sum_{L>N}\sum_{|m|\leq L}|\widehat{{\bf
u}}_{L,m}|^{2}$ $\displaystyle\leq$
$\displaystyle\lambda^{-2s-1}_{N+1}e^{-2\psi(r\cos\theta)\lambda^{1/2}_{N+1}}\sum_{L>N}\sum_{|m|\leq
L}\lambda^{2s+1}_{L}e^{2\psi(r\cos\theta)\lambda^{1/2}_{L}}|\widehat{{\bf
u}}_{L,m}|^{2}$ $\displaystyle\leq$
$\displaystyle\lambda^{-2s-1}_{N+1}e^{-2\psi(r\cos\theta)\lambda^{1/2}_{N+1}}\|\boldsymbol{\mathrm{A}}^{s+1/2}e^{\psi(r\cos\theta)\boldsymbol{\mathrm{A}}^{1/2}}{\bf
u}(\zeta)\|^{2}$ $\displaystyle\leq$ $\displaystyle
C\lambda^{-2s-1}_{N+1}e^{-2\psi(r\cos\theta)\lambda^{1/2}_{N+1}}$
$\displaystyle\leq$ $\displaystyle CN^{-2(2s+1)}e^{-2\psi(r\cos\theta)N}.$

For $\zeta=re^{i\theta}\in\Gamma$ it is easily seen that
$r\cos\theta\geq(1-1/\sqrt{2})t$, and hence that

(3.17)
$\psi(r\cos\theta)\geq\min\\{(1-1/\sqrt{2})t,T^{**},\sigma_{1}\\}=:\psi_{1}(t).$

On using (3.16) in (3.14) we get

$\left\|\frac{d}{dt}\left({\bf u}-\widetilde{{\bf
u}}_{N}\right)(t)\right\|\leq\frac{1}{2\pi}\int_{\Gamma}\frac{\|{\bf
p}_{N}(\zeta)\|}{|t-\zeta|^{2}}d\zeta\leq CN^{-2s-1}e^{-\psi_{1}(t)N}.$

∎

##  4\. A pseudospectral quadrature method

We are now ready to describe, analyze, and implement a spectrally accurate
scheme to compute approximate solutions of the NSE (2.1) in $V_{N}$, through
its weak formulation (2.27). The task is then to compute ${\bf
u}_{N}(\cdot,t)\in V_{N}$ for $t\in[0,T]$ with ${\bf u}_{N}(0,\widehat{{\bf
x}})={\mathbf{L}}_{N}{\bf u}_{0}(\widehat{{\bf x}})$, $\widehat{{\bf x}}\in
S$, satisfying the (spatially) discrete system of ordinary differential
equations

(4.1) $\frac{d}{dt}\left({\bf u}_{N},{\bf v}\right)_{M}+b({\bf u}_{N},{\bf
u}_{N},{\bf v})_{M}+\nu(\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
u}_{N},\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf
v})_{M}+(\boldsymbol{\mathrm{C}}{\bf u}_{N},{\bf v})_{M}=({\bf f},{\bf
v})_{M},$

for all ${\bf v}\in V_{N}$ and prove that the scheme is spectrally accurate
with respect to the parameter $N$ (that is, converges with rate determined by
the smoothness of the given data), and demonstrate the theory with numerical
experiments.

Using (3.2), the exactness properties of the discrete inner product, (3.8),
(2.10), (2.16), (2.23), and (2.29), the system (4.1) can be written as

(4.2) $\displaystyle\left(\frac{d}{dt}{\bf
u}_{N}+\nu\boldsymbol{\mathrm{A}}{\bf u}_{N}+\boldsymbol{\mathrm{B}}({\bf
u}_{N},{\bf u}_{N})+\boldsymbol{\mathrm{C}}{\bf
u}_{N},{\mathbf{Z}}_{L,m}\right)$ $\displaystyle=$ $\displaystyle({\bf
f},{\mathbf{Z}}_{L,m})_{M},$ $\displaystyle L=1,\cdots,N;~{}~{}m=-L,\cdots,L.$

###  4.1. Stability and convergence analysis

First we establish the stability of the approximate solution ${\bf u}_{N}$ of
(4.24). That is, similar to [6, Proposition 9.1], we prove that
$\max_{t\in[0,T]}\|{\bf u}_{N}\|_{V}$ is uniformly bounded with the bound
depending only on the initial data, forcing function, and the viscosity term
in (4.24).

######  Theorem 4.1.

Let ${\bf u}_{0}$ and ${\bf f}$ satisfy (3.1). Let $N\geq 1$ be an integer.
Let ${\bf u}_{N}$ be the solution of the pseudospectral quadrature Galerkin
system (4.1). Then there exists a constant $C$ depending on $\nu,\|{\bf
u}_{0}\|_{V}$ and $\|{\bf f}\|_{\infty}:=\max_{t\in[0,T]}\|{\bf
f}(t)\|_{C(TS)}$ so that

$\max_{t\in[0,T]}\|{\bf u}_{N}\|_{V}\leq C.$

###### Proof.

The proof follows by repeating the arguments described in the first four pages
of [6, Section 9] (proving [6, Proposition 9.1]), provided that we establish
[6, Inequalities (9.3) and (9.13)] for our system (4.1) on the spherical
surface with additional Coriolis term and quadrature approximations.

Using (2.26), (2.29) and the exactness of the quadrature rule (given by
(3.7)-(3.8)), we have $\left(\boldsymbol{\mathrm{B}}({\bf u}_{N},{\bf
u}_{N}),{\bf u}_{N}\right)_{M}=0$. Using (4.27), the symmetry of the
coefficients of ${\bf u}_{N}$ in (4.24) and the exactness of the quadrature,
we get

(4.3) $(\boldsymbol{\mathrm{C}}{\bf u}_{N},{\bf
u}_{N})_{M}=(\boldsymbol{\mathrm{C}}{\bf u}_{N},{\bf u}_{N})=(-2\Omega
i)\sum_{L=1}^{N}\lambda^{-1}_{L}\sum_{|m|\leq L}m|\alpha_{L,m}|^{2}=0.$

By taking ${\bf v}$ to be ${\bf u}_{N}$ in (4.1) and using $\|{\bf
u}_{N}\|^{2}=({\bf u}_{N},{\bf u}_{N})_{M}$, $\|{\bf
u}_{N}\|_{V}^{2}=(\boldsymbol{\mathrm{A}}{\bf u}_{N},{\bf u}_{N})_{M}$, (4.3),
Young’s inequality and the fact that all the eigenvalues $\lambda_{J}$ of
$\boldsymbol{\mathrm{A}}$ (corresponding to eigenvectors in ${\bf u}_{N}$)
satisfy $\lambda_{J}\geq\lambda_{1}=2,J=1,\cdots,N$, we obtain

$\frac{1}{2}\frac{d}{dt}\|{\bf u}_{N}\|^{2}+\nu\|{\bf u}_{N}\|_{V}^{2}=({\bf
f},{\bf u}_{N})_{M}\leq\|{\bf f}\|_{\infty}\|{\bf u}_{N}\|\leq\frac{\|{\bf
f}\|^{2}_{\infty}}{4\nu}+\nu\|{\bf u}_{N}\|^{2}\leq\frac{\|{\bf
f}\|^{2}_{\infty}}{4\nu}+\frac{\nu}{2}\|{\bf u}_{N}\|_{V}^{2},$

Hence, for our discrete system (4.1), we obtain [6, Inequality (9.3)]:

(4.4) $\frac{d}{dt}\|{\bf u}_{N}\|^{2}+\nu\|{\bf
u}_{N}\|_{V}^{2}\leq\frac{\|{\bf f}\|^{2}_{\infty}}{\nu\lambda_{1}}.$

Again using (4.27), the exactness of the quadrature rule, eigenfunction
properties of $\boldsymbol{\mathrm{A}}$, and the symmetry of the coefficients
of ${\bf u}_{N}$ in (4.24), we get

(4.5) $(\boldsymbol{\mathrm{C}}{\bf u}_{N},\boldsymbol{\mathrm{A}}{\bf
u}_{N})=(-2\Omega i)\sum_{L=1}^{N}\sum_{|m|\leq L}m|\alpha_{L,m}|^{2}=0,$

By taking ${\bf v}$ to be $\boldsymbol{\mathrm{A}}{\bf u}_{N}$ in equation
(4.1), and using $\|\boldsymbol{\mathrm{A}}{\bf
u}_{N}\|^{2}=(\boldsymbol{\mathrm{A}}{\bf u}_{N},\boldsymbol{\mathrm{A}}{\bf
u}_{N})_{M}$, (3.8), and (4.3), we obtain

(4.6) $\frac{1}{2}\frac{d}{dt}\|{\bf
u}_{N}\|^{2}_{V}+\nu\|\boldsymbol{\mathrm{A}}{\bf u}_{N}\|^{2}+b({\bf
u}_{N},{\bf u}_{N},\boldsymbol{\mathrm{A}}{\bf u}_{N})=({\bf
f},\boldsymbol{\mathrm{A}}{\bf u}_{N})_{M}.$

Using [16, Lemma 3.1], $b({\bf u}_{N},{\bf u}_{N},\boldsymbol{\mathrm{A}}{\bf
u}_{N})=0$. The term $({\bf f},\boldsymbol{\mathrm{A}}{\bf u}_{N})_{M}$ can be
estimated by using the exactness of the quadrature and Young’s inequality:

$({\bf f},\boldsymbol{\mathrm{A}}{\bf u}_{N})_{M}\leq\|{\bf
f}\|_{\infty}\|\boldsymbol{\mathrm{A}}{\bf u}_{N}\|\leq\frac{\|{\bf
f}\|^{2}_{\infty}}{2\nu}+\frac{\nu}{2}\|\boldsymbol{\mathrm{A}}{\bf
u}_{N}\|^{2}$

Hence we obtain a stronger version of [6, Inequality (9.13)] for our
quadrature discrete scheme (4.1):

$\frac{d}{dt}\|{\bf u}_{N}\|^{2}_{V}+\nu\|\boldsymbol{\mathrm{A}}{\bf
u}_{N}\|^{2}\leq\frac{\|{\bf f}\|^{2}_{\infty}}{\nu\lambda_{1}}.$

Thus, the result follows from arguments in [6, Page 74-77]. ∎

Next, using Theorem 3.1, 3.2, and 4.1, we prove the spectral convergence of
the solution ${\bf u}_{N}$ of (4.2) to the solution of ${\bf u}$ of (2.27).

######  Theorem 4.2.

Let ${\bf u}_{0}$, ${\bf f}$ satisfy (2.37), (3.1). Then there exists a
$T^{\\#}>0$, depending only on $\nu,{\bf f}$, ${\bf u}_{0}$ and the uniform
bound in (2.45) (and hence there exists
$0<\mu(t)<\min\\{t,T^{\\#},\sigma_{1}\\}$) such that for all $t\in(0,T)$,

(4.7) $\|{\bf u}-{\bf u}_{N}\|\leq
C\left[N^{-2s-1}e^{-\mu(t)N}+\left\|{\mathbf{\Pi}}_{N}{\bf
f}-{\mathbf{L}}_{N}{\bf f}\right\|\right].$

In particular, with $2s$ being an integer

(4.8) $\|{\bf u}-{\bf u}_{N}\|\leq CN^{-2s}.$

###### Proof.

Let ${\bf w}_{N}=\widetilde{{\bf u}}_{N}-{\bf u}_{N}$, where the comparison
function $\widetilde{{\bf u}}_{N}$ is the solution of (3.13). Since ${\bf
u}-{\bf u}_{N}={\bf p}_{N}+{\bf w}_{N}$, where ${\bf p}_{N}={\bf
u}-\widetilde{{\bf u}}_{N}$, in view of Theorem 3.1 and 3.2, existence of
$T^{\\#}$ and $\mu(t)$ follows and it is sufficient to show that $\|{\bf
w}_{N}\|\leq C\left[N^{-2s-1}e^{-\mu(t)N}+\left\|{\mathbf{\Pi}}_{N}{\bf
f}-{\mathbf{L}}_{N}{\bf f}\right\|\right]$. For any ${\bf v}\in V_{N}$, using
(2.28), (3.13), and (4.2),

$\displaystyle(({\bf w}_{N})_{t},{\bf v})+\nu(\boldsymbol{\mathrm{A}}{\bf
w}_{N},{\bf v})+(\boldsymbol{\mathrm{C}}{\bf w}_{N},{\bf v})$
$\displaystyle=(({\widetilde{{\bf u}}_{N}})_{t},{\bf
v})+\nu(\boldsymbol{\mathrm{A}}{\widetilde{{\bf u}}_{N}},{\bf
v})+(\boldsymbol{\mathrm{C}}{\widetilde{{\bf u}}}_{N},{\bf v})-(({\bf
u}_{N})_{t},{\bf v})-\nu(\boldsymbol{\mathrm{A}}{\bf u}_{N},{\bf
v})-(\boldsymbol{\mathrm{C}}{\bf u}_{N},{\bf v})$
$\displaystyle=(({\widetilde{{\bf u}}_{N}})_{t},{\bf
v})+\nu(\boldsymbol{\mathrm{A}}{\widetilde{{\bf u}}_{N}},{\bf
v})+(\boldsymbol{\mathrm{C}}{\widetilde{{\bf u}}}_{N},{\bf v})-({\bf f},{\bf
v})_{M}+(\boldsymbol{\mathrm{B}}({\bf u}_{N},{\bf u}_{N}),{\bf v})$
$\displaystyle=(({\widetilde{{\bf u}}_{N}})_{t},{\bf
v})+\nu(\boldsymbol{\mathrm{A}}{\bf u},{\bf
v})+(\boldsymbol{\mathrm{C}}{\widetilde{{\bf u}}}_{N},{\bf v})-({\bf f},{\bf
v})_{M}+(\boldsymbol{\mathrm{B}}({\bf u}_{N},{\bf u}_{N}),{\bf v})$
$\displaystyle=(({\widetilde{{\bf u}}_{N}}-{\bf u})_{t},{\bf v})+({\bf f},{\bf
v})-({\bf f},{\bf v})_{M}+(\boldsymbol{\mathrm{C}}{\widetilde{{\bf
u}}_{N}}-\boldsymbol{\mathrm{C}}{\bf u},{\bf v})+(\boldsymbol{\mathrm{B}}({\bf
u}_{N},{\bf u}_{N})-\boldsymbol{\mathrm{B}}({\bf u},{\bf u}),{\bf v}).$

Using the orthogonal projection ${\mathbf{\Pi}}_{N}$ in (3.2), we can write
this relation in functional form as

$\frac{d{\bf
w}_{N}}{dt}=(-\nu\boldsymbol{\mathrm{A}}-\boldsymbol{\mathrm{C}}){\bf
w}_{N}-\mathbf{{\mathbf{\Pi}}_{N}}\left[({\bf
p}_{N})_{t}+\boldsymbol{\mathrm{C}}{\bf p}_{N}\right]+{\mathbf{\Pi}}_{N}{\bf
f}-{\mathbf{L}}_{N}{\bf
f}+{\mathbf{\Pi}}_{N}\left[\boldsymbol{\mathrm{B}}({\bf u}_{N},{\bf
u}_{N})-\boldsymbol{\mathrm{B}}({\bf u},{\bf u})\right].$

Integrating with respect to $t$ and using ${\bf w}_{N}(0)=0$, we have

(4.9) $\displaystyle{\bf w}_{N}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}e^{-(t-s)(\nu\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{C}})}\left[-{\mathbf{\Pi}}_{N}(\frac{d}{ds}{\bf
p}_{N}+\boldsymbol{\mathrm{C}}{\bf p}_{N})+{\mathbf{\Pi}}_{N}{\bf
f}-{\mathbf{L}}_{N}{\bf f}\right](s)\;ds$
$\displaystyle+\int_{0}^{t}e^{-(t-s)(\nu\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{C}})}{\mathbf{\Pi}}_{N}\left[\boldsymbol{\mathrm{B}}({\bf
u}_{N},{\bf u}_{N})-\boldsymbol{\mathrm{B}}({\bf u},{\bf u})\right](s)\;ds.$

Let

(4.10)
$R_{N}(\epsilon,t-s)=\|\nu^{\epsilon}{\mathbf{\Pi}}_{N}\boldsymbol{\mathrm{A}}^{\epsilon}e^{-(t-s)(\nu\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{C}})}\|.$

On using (2.31), with $\delta\in(1/2,1)$, and the uniform boundedness of the
orthogonal projection ${\mathbf{\Pi}}_{N}$, we get

(4.11)
$\displaystyle\|e^{-(t-s)(\nu\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{C}})}{\mathbf{\Pi}}_{N}\left[\boldsymbol{\mathrm{B}}({\bf
u}_{N},{\bf u}_{N})-\boldsymbol{\mathrm{B}}({\bf u},{\bf u})\right]\|$
$\displaystyle\leq\frac{R_{N}(\delta,t-s)}{\nu^{\delta}}\|\boldsymbol{\mathrm{A}}^{-\delta}{\mathbf{\Pi}}_{N}(\boldsymbol{\mathrm{B}}({\bf
u}_{N},{\bf u}_{N})-\boldsymbol{\mathrm{B}}({\bf u},{\bf u}))\|\leq
CR_{N}(\delta,t-s)\|{\bf u}_{N}-{\bf u}\|.$

Taking norms and using $\|{\bf u}_{N}-{\bf u}\|\leq\|{\bf w}_{N}\|+\|{\bf
p}_{N}\|$ together with (4.10), and (4.1) in (4.9), we obtain

$\displaystyle\|{\bf w}_{N}(t)\|$ $\displaystyle\leq$
$\displaystyle\left[\|\frac{d}{dt}{\bf p}_{N}\|+\|\boldsymbol{\mathrm{C}}{\bf
p}_{N}\|+\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\mathbf{L}}_{N}{\bf
f}\right\|\right]\int_{0}^{t}R_{N}(0,t-s)\;ds$
$\displaystyle+C\int_{0}^{t}R_{N}(\delta,t-s)(\|{\bf p}_{N}(s)\|+\|{\bf
w}_{N}(s)\|)ds.$

Using Gronwall’s inequality, we obtain for each $t\in[0,T]$,

(4.12) $\displaystyle\|{\bf w}_{N}(t)\|$ $\displaystyle\leq$ $\displaystyle
C\left\\{\left[\|\frac{d}{dt}{\bf p}_{N}\|+\|\boldsymbol{\mathrm{C}}{\bf
p}_{N}\|+\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\mathbf{L}}_{N}{\bf
f}\right\|\right]\int_{0}^{t}R_{N}(0,t-s)\;ds\right.$
$\displaystyle\qquad\left.+\|{\bf p}_{N}\|\
\int_{0}^{t}R_{N}(\delta,t-s)\;ds\right\\}.$

For $\epsilon\in[0,1]$, we have to bound

(4.13)
$\int_{0}^{t}R_{N}(\epsilon,t-s)\;ds=\int_{0}^{t}R_{N}(\epsilon,r)\;dr.$

Using also (4.10) and (4.27),

$R_{N}(\epsilon,r)\leq\max_{1\leq L\leq N;|m|\leq
L}|(\nu\lambda_{L})^{\epsilon}e^{-\nu\lambda_{L}r}e^{-2\Omega
irm\lambda^{-1/2}_{L}}|\leq\max_{z\in[\nu\lambda_{1},\nu\lambda_{N}]}z^{\epsilon}e^{-rz}.$

Thus

(4.14)
$R_{N}(r)\leq\begin{cases}(\nu\lambda_{N})^{\epsilon}e^{-\nu\lambda_{N}r}&\mbox{
if
}r\leq\epsilon/(\nu\lambda_{N}),\cr\epsilon^{\epsilon}e^{-\epsilon}r^{-\epsilon}&\mbox{
if }\epsilon/(\nu\lambda_{N})\leq
r\leq\epsilon/(\nu\lambda_{1}),\cr(\nu\lambda_{1})^{\epsilon}e^{-\nu\lambda_{1}r}&\mbox{
if }r\geq\epsilon/(\nu\lambda_{1}).\cr\end{cases}$

With

(4.15) $I_{1}=[0,\epsilon/(\nu\lambda_{N})]\cap[0,t],\quad
I_{2}=[\epsilon/(\nu\lambda_{N}),\epsilon/(\nu\lambda_{1})]\cap[0,t],\quad
I_{3}=[\epsilon/(\nu\lambda_{1}),t]\cap[0,t],$

the interval of integration in (4.13) can be subdivided into these three
intervals. In particular, using (4.14) and (4.15), we get

(4.16)
$\int_{I_{1}}R_{N}(r)dr\leq\int_{I_{1}}(\nu\lambda_{N})^{\epsilon}e^{-\nu\lambda_{N}r}dr\leq\frac{(1-e^{-\epsilon})}{(\nu\lambda_{N})^{1-\epsilon}}\leq
C,$ (4.17)
$\int_{I_{2}}R_{N}(r)dr\leq\int_{I_{2}}\epsilon^{\epsilon}e^{-\epsilon}r^{-\epsilon}dr\leq\frac{\epsilon
e^{-\epsilon}}{(1-\epsilon)}\left[\left(\frac{1}{\nu\lambda_{1}}\right)^{(1-\epsilon)}-\left(\frac{1}{\nu\lambda_{N}}\right)^{(1-\epsilon)}\right]\leq
C,$

and

(4.18)
$\int_{I_{3}}R_{N}(r)dr\leq\int_{I_{3}}(\nu\lambda_{1})^{\epsilon}e^{\nu\lambda_{1}r}dr\leq\frac{(e^{-\epsilon}-e^{-\nu\lambda_{1}t})}{(\nu\lambda_{1})^{1-\epsilon}}\leq
C.$

Using (4.16)- (4.18) in (4.13), we get

(4.19) $\int_{0}^{t}R_{N}(\epsilon,t-s)\;ds\leq C,\qquad\epsilon\in[0,1].$

Substituting this in (4.12), we get

(4.20) $\|{\bf w}_{N}(t)\|\leq C\left[\|\frac{d}{dt}{\bf
p}_{N}\|+\|\boldsymbol{\mathrm{C}}{\bf p}_{N}\|+\|{\bf p}_{N}\|\
+\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\mathbf{L}}_{N}{\bf f}\right\|\right].$

The term $\|\boldsymbol{\mathrm{C}}{\bf p}_{N}\|$ in (4.20) can be simplified
using (2.23) and the fact that ${\bf p}_{N}$ is tangential (and hence
$\widehat{{\bf x}}\cdot{\bf p}_{N}(\widehat{{\bf x}})=0$),

$\left[\widehat{{\bf x}}\times{\bf p}_{N}(\widehat{{\bf
x}})\right]\cdot\left[\overline{(\widehat{{\bf x}}\times{\bf
p}_{N}(\widehat{{\bf x}})}\right]=\left[\widehat{{\bf x}}\cdot\widehat{{\bf
x}}\right]\left[{\bf p}_{N}(\widehat{{\bf x}})\cdot\overline{{\bf
p}_{N}(\widehat{{\bf x}})}\right]=\|{\bf p}_{N}\|^{2},$

and hence

(4.21) $\|{\bf w}_{N}(t)\|\leq C\left[\|\frac{d}{dt}{\bf p}_{N}\|+\|{\bf
p}_{N}\|\ +\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\mathbf{L}}_{N}{\bf
f}\right\|\right].$

Hence from Theorem 3.1 and 3.2, for all $t\in(0,T)$ we have

(4.22) $\|{\bf u}-{\bf u}_{N}\|\leq
C\left[N^{-2s-1}e^{-\mu(t)N}+\left\|{\mathbf{\Pi}}_{N}{\bf
f}-{\mathbf{L}}_{N}{\bf f}\right\|\right].$

In particular, using (3.5) and (3.12) with $2s$ being an integer, we get

(4.23) $\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\mathbf{L}}_{N}{\bf
f}\right\|\leq\left\|{\mathbf{\Pi}}_{N}{\bf f}-{\bf f}\right\|+\left\|{\bf
f}-{\mathbf{L}}_{N}{\bf f}\right\|\leq CN^{-2s}.$

Now the result (4.8) follows from (4.22) and (4.23). ∎

###  4.2. Adaptive and fast implementation of the pseudospectral method

Having established spectrally accurate convergence of the spatially discrete
scheme (4.2), for implementation of the scheme (4.2) to simulate stable and
accurate solutions of (2.1) and compare with benchmark random flow simulations
in the literature, we need to discretize the time derivative operator
$\frac{d}{dt}$ in (4.2). Further, at each discrete time step we develop a
(FFT-based) fast evaluation technique to set up the resulting fully discrete
nonlinear system with spatial $\mathcal{O}(N^{4})$ complexity. First we
consider discretization of $\frac{d}{dt}$ in (4.2).

In order to the make (4.2) fully discrete in space and time, and hence compute
the $N^{2}+2N$ unknown time-dependent coefficients in the representation of
the tangential divergence-free approximate real velocity vector field

(4.24) ${\bf u}_{N}(\widehat{{\bf x}},t):=\sum_{L=1}^{N}\sum_{|m|\leq
L}\alpha_{L,m}(t){\mathbf{Z}}_{L,m}(\widehat{{\bf
x}}),\quad\alpha_{L,m}=\overline{\alpha}_{L,-m},\quad\alpha_{L,m}(0)=({\bf
u}_{0},{\mathbf{Z}}_{L,m})_{M},$

for $\widehat{{\bf x}}\in S$ and $t\geq 0$, the standard fixed-time-step
backward-Euler (or Crank-Nicolson) Galerkin approach could be used in (4.1),
leading to a first-order (or a second-order, respectively) in time non-
adaptive scheme [13]. However, due to the complicated unknown flow behavior of
the NSE solutions, when the initial states are random, it is more efficient
instead to integrate (4.1) using a combination of multi-order integration
formulas that allow adaptive choice of time step, leading to computation of
solutions with a specified accuracy in time. In this paper we follow the
latter approach.

For implementation purposes, we first need to substitute (4.24) in (4.2). We
write the resulting $N^{2}+2N$-dimensional system of ordinary differential
equations (ODE), for the unknown $N^{2}+2N$ time dependent coefficients of
${\bf u}_{N}(\widehat{{\bf x}},t)$ in (4.24), as

(4.25)
$\frac{d}{dt}\boldsymbol{\alpha}(t)=\mathbf{F}(t,\boldsymbol{\alpha}(t)).$

It is well known that such nonlinear ODE systems are stiff, and hence it is
important to use time discretization techniques with a large stability region
[14]. Further, it is important to use high-order implicit formulas whenever
possible, but the high-order discretization formulas are appropriate only at
those time steps where the unknown exact solution is smooth.

For practical problems such as the Navier-Stokes equations (with initial
random state) where the necessary spatial discretization (at each time step)
is expensive, it is important to optimize computing time by simulating only up
to a required accuracy by choosing adaptive discrete time steps. Unlike the
adaptive spatial mesh for elliptic PDEs based on the a posteriori estimates,
adaptive time steps can be computed by comparing numerical solutions obtained
using two distinct order formulas [14] and hence simulation using multi-order
formulas is appropriate.

In particular, for practical realization of large stiff nonlinear ODE systems,
multi-order implicit backward differentiation formulas (BDF) and their
generalizations such as the numerical differential formulas (NDF) are most
appropriate. The implicit NDF formula of order $p$ (NDFp) with a parameter
$\kappa_{p}$ (so that $\kappa_{p}=0$ corresponds to BDFp) for the system
(4.25), with $\boldsymbol{\alpha}_{n}\approx\boldsymbol{\alpha}(t_{n})$ and
$\nabla^{m}$ denoting the $m$-th Newton backward difference operator, is

(4.26)
$\sum_{m=1}^{p}\frac{1}{m}\nabla^{m}\boldsymbol{\alpha}_{n+1}=h\mathbf{F}(t_{n+1},\boldsymbol{\alpha}_{n+1})+\kappa_{p}\nabla^{p+1}\boldsymbol{\alpha}_{n+1}\sum_{j=1}^{p}\frac{1}{j}.$

It is well known [14] that BDFp (and hence NDFp) are unstable for $p>6$, and
for $p=6$ the stability region is small and hence not practically useful in
our case. Further the celebrated Dahlquist barrier [14] implies that BDFp (and
hence NDFp) cannot be absolutely stable [that is, $A(\alpha)$-stable with
$\alpha=90\,^{\circ}$] for $p>2$.

Following details in [25], for simulation of (4.25) we use multi-order NDFp
with $p=1,2,3,4,5$ (and respective
$\kappa_{p}=-0.1850,-1/9,-0.0823,-0.0415,0$) and these are
$A(\alpha_{p})$-stable, with respective
$\alpha_{p}=90^{\circ},90^{\circ},80^{\circ},66^{\circ},51^{\circ}$. For
$p=1,2,3,4,5$, NDFp is more accurate than BDFp, however NDFp has slightly
smaller stability angle compared BDFp only for $p=3,4$ (with respective
$\alpha_{p}=86^{\circ},73^{\circ}$) and the same stability angle for
$p=1,2,5$.

For each fixed time discretization step, the computational cost is dominated
by evaluation of $\mathbf{F}(t_{n+1},\boldsymbol{\alpha}_{n+1})$ in (4.26) and
it is important to have an efficient method to set up the spatial part of the
nonlinear system (4.2). Using the spectral properties of the Stokes operator
$\boldsymbol{\mathrm{A}}$ given by (2.17), the linear second term in (4.2) is
trivial to evaluate using the diagonal matrix consisting of the eigenvalues of
$\boldsymbol{\mathrm{A}}$.

The Coriolis term can be evaluated similarly using the identity [5, Equation
(24)]

(4.27)
$(\boldsymbol{\mathrm{C}}\;\boldsymbol{\mathrm{Curl}}\,Y_{J,k},{\mathbf{Z}}_{L,m})=(2\Omega\cos\theta\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Curl}}\,Y_{J,k},{\mathbf{Z}}_{L,m})=-2\Omega
i\frac{m}{\lambda^{1/2}_{L}}\delta_{L,J}\delta_{k,m}.$

For the first component in the nonlinear third term in (4.2), we use (2.7) to
write

(4.28)
$\displaystyle\boldsymbol{\mathrm{B}}({\mathbf{Z}}_{R,s},{\mathbf{Z}}_{J,k})$
$\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{\lambda_{R}\lambda_{J}}}{\mathbf{P}}_{\boldsymbol{\mathrm{Curl}}\,}(\Delta
Y_{R,s}\boldsymbol{\mathrm{Grad}}\,Y_{J,k}),$ (4.29)
$\displaystyle\left(\boldsymbol{\mathrm{B}}({\mathbf{Z}}_{R,s},{\mathbf{Z}}_{J,k}),{\mathbf{Z}}_{L,m}\right)$
$\displaystyle=$
$\displaystyle\sqrt{\frac{\lambda_{R}}{\lambda_{J}\lambda_{L}}}\left(Y_{R,s}\boldsymbol{\mathrm{Grad}}\,Y_{J,k},\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Grad}}\,Y_{L,m}\right).$

It is convenient to write $\boldsymbol{\mathrm{Grad}}\,Y_{J,k}$ and
$\widehat{{\bf x}}\times\boldsymbol{\mathrm{Grad}}\,Y_{L,m}$ in terms of
expressions similar to those in (2.18). Such explicit representations are also
useful for the efficient evaluation of the $N^{2}$ Fourier coefficients $({\bf
f},{\mathbf{Z}}_{L,m})_{M}$ of the source term in (4.2), and eventually for
the computation of the vorticity field.

In order to express the tangential (and normal, needed for computing the
approximate vorticity from ${\bf u}_{N}$) vector harmonics as a linear
combination of the scalar harmonics (2.18), we first recall, from the
classical quantum mechanics literature (see, for example, [30]), the covariant
spherical basis vectors

(4.30) ${\bf e}_{+1}=-\frac{1}{\sqrt{2}}([1,0,0]^{T}+i[0,1,0]^{T}),\quad{\bf
e}_{0}=[0,0,1]^{T},\quad{\bf
e}_{-1}=\frac{1}{\sqrt{2}}([1,0,0]^{T}-i[0,1,0]^{T}),$

and the Clebsch-Gordan coefficients

(4.31)
$C^{j,m}_{j_{1},m_{1},j_{2},m_{2}}:=(-1)^{(m+j_{1}-j_{2})}\sqrt{2j+1}\left(\begin{array}[]{lll}j_{1}&j_{2}&j\\\
m_{1}&m_{2}&-m\end{array}\right),$

where $\left(\begin{array}[]{lll}a&b&c\\\
\alpha&\beta&\gamma\end{array}\right)$ are the Wigner 3j-symbols given, for
example, by the Racah formula,

$\displaystyle\left(\begin{array}[]{lll}a&b&c\\\
\alpha&\beta&\gamma\end{array}\right)$
$\displaystyle=(-1)^{(a-b-\gamma)}\sqrt{T(abc)}\sqrt{(a+\alpha)!(a-\alpha)!(b+\beta)!(b-\beta)!(c+\gamma)!(c-\gamma)!}\times$
$\displaystyle\sum_{t}\frac{(-1)^{t}}{t!(c-b+t+\alpha)!(c-a+t-\beta)!(a+b-c-t)!(a-t-\alpha)!(b-t+\beta)!},$

where the sum is over all integers $t$ for which the factorials in the
denominator all have nonnegative arguments. In particular, the number of terms
in the sum is $1+\min\\{a\pm\alpha,b\pm\beta,c\pm\gamma,a+b-c,b+c-a,c+a-b\\}$.
The triangle coefficient $T(abc)$ is defined by

$T(abc)=\left[\frac{(a+b-c)!(a-b+c)!(-a+b+c)!}{(a+b+c+1)!}\right].$

Below, we require $C^{j,m}_{j_{1},m_{1},j_{2},m_{2}}$ only for some
$j_{2},m_{2}\in\\{-1,0,1\\}$, and using various symmetry and other known
properties (such as $C^{j,m}_{j_{1},m_{1},j_{2},m_{2}}=0$ unless the
conditions $|j_{1}-j_{2}|\leq j\leq j_{1}+j_{2}$ and $m_{1}+m_{2}=m$ hold) of
Wigner 3j-symbols, these coefficients can be efficiently pre-computed and
stored.

In our computation, we used the following basis functions for the tangential
vector fields: $(i)~{}\boldsymbol{\mathrm{Grad}}\,Y_{L,m},(ii)~{}\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Grad}}\,Y_{J,k}$. For the vorticity components of
${\bf u}_{N}$, in addition we used
$(iii)~{}\boldsymbol{\mathrm{Vort}}\,{\mathbf{Z}}_{J,m}=\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf
x}}}\times{\mathbf{Z}}_{J,m}=\lambda_{J}^{1/2}\widehat{{\bf
x}}Y_{J,m}=-\lambda_{J}^{-1/2}\widehat{{\bf x}}\Delta Y_{J,m}$. In particular,
using (4.24) our approximation to the vorticity in (2.2), for a fixed $t\geq
0$ and $\widehat{{\bf x}}\in S$, is

(4.33) $\boldsymbol{\mathrm{Vort}}\,{\bf u}_{N}(\widehat{{\bf
x}},t)=\boldsymbol{\mathrm{Curl}}_{\widehat{{\bf x}}}{\bf u}_{N}(\widehat{{\bf
x}})=\widehat{{\bf x}}\Delta\Psi_{N}(\widehat{{\bf x}}),$

where

$\Psi_{N}(\widehat{{\bf x}},t)=-\sum_{L=1}^{N}\sum_{|m|\leq
L}\lambda_{L}^{-1/2}\alpha_{L,m}(t)Y_{L,m}(\widehat{{\bf x}}).$

To facilitate easy application of fast transforms to evaluate these functions
at the $M=\mathcal{O}(N^{2})$ quadrature points
$\\{\widehat{\boldsymbol{\xi}}_{p,q}=\boldsymbol{p}(\theta_{p},\phi_{q})\\}$,
we represent these three types of fields first as a linear combination of the
covariant vectors in (4.30):

(4.34) $\displaystyle\boldsymbol{\mathrm{Grad}}\,Y_{L,m}$ $\displaystyle=$
$\displaystyle B_{+1,L,m}{\bf e}_{+1}+B_{0,L,m}{\bf e}_{0}+B_{-1,L,m}{\bf
e}_{-1},$ (4.35) $\displaystyle(\widehat{{\bf
x}}\times\boldsymbol{\mathrm{Grad}}\,Y_{J,k})$ $\displaystyle=$ $\displaystyle
D_{+1,J,k}{\bf e}_{+1}+D_{0,J,k}{\bf e}_{0}+D_{-1,J,k}{\bf e}_{-1},$

With $c_{L}=(L+1)\sqrt{\frac{L}{2L+1}},d_{L}=L\sqrt{\frac{L+1}{2L+1}}$, these
coefficients are explicitly given by

$\displaystyle B_{+1,L,m}$ $\displaystyle=$
$\displaystyle\left\\{c_{L}C^{L,m}_{L-1,m-1,1,1}P^{m-1}_{L-1}(\cos\theta)+d_{L}C^{L,m}_{L+1,m-1,1,1}P^{m-1}_{L+1}(\cos\theta)\right\\}e^{i(m-1)\varphi}$
$\displaystyle B_{0,L,m}$ $\displaystyle=$
$\displaystyle\left\\{c_{L}C^{L,m}_{L-1,m,1,0}P^{m}_{L-1}(\cos\theta)+d_{L}C^{L,m}_{L+1,m,1,0}P^{m}_{L+1}(\cos\theta)\right\\}e^{im\varphi}$
$\displaystyle B_{-1,L,m}$ $\displaystyle=$
$\displaystyle\left\\{c_{L}C^{L,m}_{L-1,m+1,1,-1}P^{m+1}_{L-1}(\cos\theta)+d_{L}C^{L,m}_{L+1,m+1,1,-1}P^{m+1}_{L+1}(\cos\theta)\right\\}e^{i(m+1)\varphi}.$
$\displaystyle D_{+1,J,k}$ $\displaystyle=$ $\displaystyle
i\sqrt{\lambda_{J}}C^{J,k}_{J,k-1,1,1}P^{k-1}_{J}(\cos\theta)e^{i(k-1)\varphi},$
$\displaystyle D_{0,J,k}$ $\displaystyle=$ $\displaystyle
i\sqrt{\lambda_{J}}C^{J,k}_{J,k,1,0}P^{k}_{J}(\cos\theta)e^{ik\varphi},$
$\displaystyle D_{-1,J,k}$ $\displaystyle=$ $\displaystyle
i\sqrt{\lambda_{J}}C^{J,k}_{J,k+1,1,-1}P^{k+1}_{J}(\cos\theta)e^{i(k+1)\varphi}.$

Noting (i) the complex azimuthal exponential terms
$e^{ik\varphi},e^{im\varphi}$ in (2.18) and (4.34)-(4.35) (via the above
representations for $B$ and $D$) for $|k|\leq J,|m|\leq L;1\leq L,J\leq N$,
and (ii) the need to evaluate several $\mathcal{O}(N^{2})$ sums, of the form
in (4.24) and (4.28)-(4.28), at the equally spaced $\mathcal{O}(N)$ azimuthal
quadrature points (see (3.6)), we reduce the complexity by $\mathcal{O}(N)$ in
each of these sums, at each adaptive-time step (described below), by using the
FFT for setting up the nonlinear system (4.2), similar to the approach in [4].
In our numerical experiments (see Section 5) for adaptive-time simulation of a
flow induced by random initial states, we observed that such an efficient FFT
based implementation reduced the (non-FFT code) computing time substantially
for the case $N=100$, to simulate from $t=0$ to $t=60$.

In addition, by using the fast Legendre/spherical transforms along the
latitudinal direction (obtained, for example, by modifying the NFFT algorithm
in [21] for evaluation of the Legendre functions in the above terms at
$\mathcal{O}(N)$ non-uniform latitudinal quadrature points), we could reduce
the complexity by $\mathcal{O}(N^{2})$. We did not use the fast
Legendre/spherical transforms in our implementation due to the spectral
convergence of the scheme and the fact that $N\leq 100$ in our simulations.
(In these complexity counts, we ignored $\mathcal{O}(\log N)$ and
$\mathcal{O}(\log^{2}N)$ terms.)

##  5\. Numerical Experiments

We demonstrate the fully discrete pseudospectral quadrature algorithm by
simulating (i) a known solution test case with low to high frequency modes and
(ii) a benchmark example [8, page 305] in which the unknown velocity and
vorticity fields are generated by a random initial state.

The first test example is useful to demonstrate that the pseudospectral
quadrature algorithm reproduces any number of high frequency modes in the
solution (provided $V_{N}$ contains all these modes), with computational error
dominated only by the chosen accuracy for the adaptive time evolution for the
ordinary differential system (4.25).

###  5.1. Example 1.

Our test case first example is (2.1) with

(5.1) ${\bf u}|_{t=0}(\widehat{{\bf x}})={\bf u}_{0}(\widehat{{\bf
x}})=g(0)[{\mathbf{W}}_{1}(\widehat{{\bf x}})-{\mathbf{W}}_{2}(\widehat{{\bf
x}})],$ (5.2) $g(t)=\nu
e^{-t}[\sin(5t)+\cos(10t)],\quad{\mathbf{W}}_{1}={\mathbf{Z}}_{1,0}+2\Re({\mathbf{Z}}_{1,1}),\;{\mathbf{W}}_{2}={\mathbf{Z}}_{2,0}+2\Re({\mathbf{Z}}_{2,1}+{\mathbf{Z}}_{2,2}),$

where ${\mathbf{Z}}_{L,m}$ is given by (2.17), $\Re(\cdot)$ denotes the real-
part function, and the external force ${\bf f}(\widehat{{\bf x}},t)$ in (2.1)
is chosen so that

(5.3) ${\bf u}(\widehat{{\bf x}},t)=tg(t)\sum_{L=1}^{N_{0}}\left[{\bf
Z}_{L,0}+2\sum_{m=1}^{L}\Re({\mathbf{Z}}_{L,m})\right](\widehat{{\bf
x}})+g(t){\mathbf{W}}_{1}(\widehat{{\bf
x}})+(t-1)g(t){\mathbf{W}}_{2}(\widehat{{\bf x}}),$

is the exact tangential divergence-free velocity field, solving the NSE (2.1).
The exact test field (5.2)-(5.3) has high oscillations both in space and time,
and exponentially decays in time. Note the dependence on a parameter $N_{0}$,
the maximum order of the spherical harmonics in the exact solution.

In our calculation of the approximate solution ${\bf u}_{N}$, we chose
$N=N_{0}$, so that all frequencies of the exact solution can be recovered. The
solution (5.3) is then used to validate our algorithm and code by numerical
adaptive time-integration of (4.1), for various values of $N=N_{0}$. In
particular, for a fixed integration tolerance error, we demonstrate in Figure
1 that all $N$ modes in (5.3) can indeed be recovered by the approximate
solution ${\bf u}_{N}$, within the chosen error tolerance, for all
$N=N_{0}=70,80,90,100$.

Figure 1. $\|{\bf u}-{\bf u}_{N}\|$ for Example 1 with a fixed time-
integration error, $N=N_{0}=70,80,90,100$.

###  5.2. Example 2.

Having established the validity of our algorithm for a simple known solution,
we use the same code to simulate unknown velocity and vorticity fields
generated by (2.1) with random initial velocity field as in [8, 12] with
angular velocity of the rotation $\Omega=1$ and $\nu^{-1}=10,000$ in (2.1).

The initial flow and external force in this benchmark example satisfy the main
assumption (3.1) for any $s,\sigma_{1}>0$ and hence, as proved in Theorem 4.2,
the approximate solution ${\bf u}_{N}$ is spectrally accurate and converges
super-algebraically with order given by (4.8) for any $s>0$. As mentioned in
Section 1, this is the main advantage of the present paper over the recent
paper [12], where such spectrally accurate convergence results are neither
discussed nor proved. On the other hand, convergence results for two-
dimensional problems on a Euclidean plane, supported by numerical experiments,
have formed a core part of research on the NSE over the last few decades, see
[8, 9, 27] and references therein.

The random initial tangential divergence-free velocity field, having
properties similar to those considered in [8, page 305] and [12, page 988]
(but not exactly same as in [8, 12], due to randomness), is a smooth function
${\bf u}_{0}={\bf v}\in G^{s/2}_{\sigma}$, a Gevrey class of order $s$ and
index $\sigma$, see (2.33), for any $s,\sigma>0$, with Fourier coefficients
$\widehat{{\bf v}}_{L,m},~{}~{}L=1,2,\cdots,|m|\leq L$, defined by

(5.4) $\widehat{{\bf
v}}_{L,m}=\begin{cases}a_{L}\exp(i\phi_{m})&L=1,\ldots,20;~{}m=0,\ldots,L,\cr
a_{L}(-1)^{m}\exp(-i\phi_{m})&L=1,\ldots,20;~{}m=-L,\cdots,-1,\cr
0,&L>20;~{}m=-L,\ldots,L,\end{cases}$

where $\phi_{0}=0$, and $\phi_{m}\in(0,2\pi)$ are random numbers for $m>0$,
and $a_{L}=b_{L}/||\mathbf{b}||$, with $\mathbf{b}\in\mathbb{R}^{20}$ having
components $b_{L}=2/\left[L+(\nu L)^{2.5}\right],~{}L=1,\ldots,20$. The
vorticity stream function $\Psi$ (see (2.2)) of
$\boldsymbol{\mathrm{Vort}}\,{\bf u}(\cdot,0)$ in Figure 3 demonstrates the
randomness of the field at time $t=0$.

The external force field ${\bf f}={\bf f}(\widehat{{\bf x}},t)$ in (2.1) for
our simulation is motivated by that considered in [8, page 305] and is exactly
same as that in [12, page 988]. The source term ${\bf f}$ is a decaying
tangential divergence-free field which, for any $s,\sigma>0$, belongs to
$C([0,T];\mathcal{D}(\boldsymbol{\mathrm{A}}^{s+1/2}e^{\sigma\boldsymbol{\mathrm{A}}^{1/2}}))$
(for any $T>0$) with the only non-zero Fourier coefficient being
$\widehat{{\bf f}(t)}_{3,0}$. The Fourier coefficient $\widehat{{\bf
f}(t)}_{3,0}$ is a continuous function in time and is defined by

$\widehat{{\bf f}(t)}_{3,0}=\begin{cases}1,&0\leq t\leq 10,\cr\cos(\pi
t/5)\exp(-(t-10)/5)&t>10.\cr\end{cases}$

The impact of the external force on the numerical velocity and its complement
(generating the approximate inertial manifold) is demonstrated in Figure 3
with the time evolution of the velocity field matching that of the external
force, leading to little change in evolution process of the velocity as the
external force gets smaller and smaller.

We chose a fixed relative error tolerance to be of accuracy at least
$\mathcal{O}(10^{-3})$, for adaptive time-integration solving the
$N^{2}+2N$-dimensional system ordinary differential equations (4.1) in time,
using backward differential formulas with variable order (one to five) and
variable adaptive time step sizes that meet the fixed error tolerance.

As discussed in the next subsection, the high-frequency components of the
solution turn out to be unreliable as $t$ increases, [8], presumably because
of the time discretization error. We therefore retained only the frequency
components up to some order $N_{1}\leq N$, correspondingly, we define an
additional approximation of ${\bf u}_{N}$, defined in (4.1)-(4.24), by

(5.5) ${\bf u}_{N_{1};N}(\cdot,t):={\mathbf{\Pi}}_{N_{1}}{\bf
u}_{N}(\cdot,t),\qquad N_{1}\leq N.$

As in Theorem 4.2, for the spatially discrete pseudospectral quadrature
method, assuming (3.1) and exact time integration of (4.1), and hence using
(2.35), (3.4), and (4.8), we get spectral convergence:

$\|{\bf u}-{\bf u}_{N_{1};N}\|\leq\left\|{\bf u}-{\mathbf{\Pi}}_{N_{1}}{\bf
u}\right\|+\left\|{\mathbf{\Pi}}_{N_{1}}\left({\bf u}-{\bf
u}_{N}\right)\right\|\leq C\left[N_{1}^{-2(s+1)}+N^{-2s}\right]\leq
CN_{1}^{-2s}.$

Our simulated approximate velocity fields

${\bf U}_{75}(t):={\bf u}_{75;100}(t),\quad\mbox{ for }t=10,20,30,40,50,60,$

are in Figure 6–9 and the associated vorticity stream function $\Psi_{75}(t)$
of $\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}(t)$ (computed using (4.33)) are
in Figure 12–15. These figures demonstrate that the initial random flow with
several smaller structures evolve into regular flow with larger structures,
similar to those observed in [8, page 307]. The choice of $N_{1}$ will be
discussed in the next subsection.

###  5.3. Energy spectrum of the solution

If ${\bf u}$ (and hence the number of modes in ${\bf u}$) is unknown, as it is
in the case in the benchmark test Example 2, and if ${\bf u}_{N}$ does not
contain all of the modes in ${\bf u}$, the higher modes of ${\bf u}_{N}$ are
usually less accurate than the chosen practical error tolerance and hence can
even violate important physical properties of ${\bf u}$ (because of the error
time-integration tolerance being not very small). In such cases, it is
important to choose $N_{1}<N$, depending on certain known physical properties
of ${\bf u}$.

For a fixed time $t$, the $L$-th mode energy spectrum of a tangential
divergence-free flow ${\bf u}$ on the sphere is defined by

(5.6) $E(L)=E({\bf u},L)=\sum_{|m|\leq L}|\beta_{L,m}(t)|^{2},\qquad\qquad{\bf
u}(\widehat{{\bf x}},t):=\sum_{L=1}^{\infty}\sum_{|m|\leq
L}\beta_{L,m}(t){\mathbf{Z}}_{L,m}(\widehat{{\bf x}}).$

Although the analytical form of the flow in Example 2 is not known, several
investigations have been carried out for such fields with initial spectrum of
the $L$-th mode decaying with order $L^{-1}$ or $L^{-2}$. In particular, it is
well known (see [8]), for this benchmark test case (on periodic two
dimensional geometries), that the energy spectrum of the velocity has a power-
law inertial range and an exponential decay (dissipation range) for wave
numbers larger than the Kraichnan’s dissipation wave number. Further, several
random smaller structures built in the initial random vorticity evolve into
regular flow with larger structures.

With ${\bf u}$ being the unique solution of the NSE (2.1), let us decompose
${\bf u}=\widetilde{{\bf u}}_{N_{1}}+{\bf w}_{N_{1}}$, where $\widetilde{{\bf
u}}_{N_{1}}={\mathbf{\Pi}}_{N_{1}}({\bf u})$ contains all modes lower or equal
$N_{1}$ and ${\bf w}_{N_{1}}:={\bf u}-{\mathbf{\Pi}}_{N_{1}}({\bf u})$
contains all higher modes. The existence of a relation between ${\bf
w}_{N_{1}}$ and $\widetilde{{\bf u}}_{N_{1}}$ of a form ${\bf
w}_{N_{1}}=\Phi(\widetilde{{\bf u}}_{N_{1}})$ was established in [28]. The
graph of $\Phi$ is known as the inertial manifold of (2.1). For computational
purposes, the higher modes can be computed efficiently using an approximate
inertial manifold:

(5.7) $\widetilde{\Phi}_{\widetilde{N_{1}}}(\widetilde{{\bf
u}}_{N_{1}})=(\nu\boldsymbol{\mathrm{A}}+\boldsymbol{\mathrm{C}})^{-1}\left({\mathbf{\Pi}}_{\widetilde{N_{1}}}-{\mathbf{\Pi}}_{N_{1}}\right)\left[{\bf
f}-\boldsymbol{\mathrm{B}}\left(\widetilde{{\bf u}}_{N_{1}},\widetilde{{\bf
u}}_{N_{1}}\right)\right],~{}~{}\widetilde{N_{1}}>N_{1}.$

This well known approximation (without the Coriolis term) was introduced in
[10, 11] for nonlinear dissipative systems, including the NSE on domains and
we choose $\widetilde{N_{1}}=2N_{1}$.

For the benchmark test case, in the dissipation range, even
$L^{4}E(\widetilde{\Phi}_{2N_{1}}(\widetilde{{\bf u}}_{N_{1}}),L)$ decays
exponentially to zero, further justifying restriction of the infinite
dimensional range after certain values of $N_{1}$. As discussed in [8, page
280, 307] (and repeated in [12, page 991]), such a faster decay (one order
higher than the $L^{-3}$ decay known is Kraichnan theory of turbulence) is
expected due to the Reynolds number considered in [8, 12] being much small
than $25,000$. (Extensive study in [3, 23] shows that the turbulence theory
decay can occur only when the Reynolds number is of the order $25,000$.)

Briefly (without repeating technical details in [8]), using the viscosity term
$\nu$ in (2.1), the Reynolds number is $\mathcal{O}(\nu^{-1})$ with the order
constant $r$ given by the product of the mean fluid velocity and the
characteristic length-scale. With $\nu^{-1}=10,000$ and the constant rotation
rate $\Omega=1$, for the Coriolis parameter in (2.23), the decay of energy
spectrum of a numerical velocity and exponential decay of its approximate
complement in Figure 16-17 is well supported by extensive simulations in [3,
23, 8], highlighting further the benchmark applicability of our algorithm,
extending periodic domain results in [3, 23, 8] to a practically relevant
rotating sphere case with Coriolis effect.

Further, the simulated results are substantiated by the well known exponential
decay of $E(\widetilde{\Phi}_{2N_{1}}(\widetilde{{\bf u}}_{N_{1}}),L)$,
observed in Figures 16-17. Finally, the exponential decay of the energy
spectrum show that $N=100,N_{1}=3N/4$ is sufficient to understand the flow
behavior for this benchmark example with our algorithm using adaptive variable
order and variable time-step highly stable backward differentiation formulas
with a practically useful relative error tolerance $\mathcal{O}(10^{-3})$.

Figure 2. Initial random vorticity stream function $\Psi$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf u}$ at $t=0$.

|
---|---
Time evolution of ${\bf U}_{75}(t)$ and $f(t)$. | Time evolution of $\widetilde{\Phi}_{150}({\bf U}_{75}(t))$ and ${\bf f}(t)$.

Figure 3. Impact of the external force on a numerical velocity and its
complement.

Figure 4. Numerical velocity ${\bf U}_{75}(t)$, at $t=10$.

Figure 5. Numerical velocity ${\bf U}_{75}(t)$, at $t=20$.

Figure 6. Numerical velocity ${\bf U}_{75}(t)$, at $t=30$.

Figure 7. Numerical velocity ${\bf U}_{75}(t)$, at $t=40$.

Figure 8. Numerical velocity ${\bf U}_{75}(t)$, at $t=50$.

Figure 9. Numerical velocity ${\bf U}_{75}(t)$, at $t=60$.

Figure 10. Vorticity stream function $\Psi_{75}$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=10$.

Figure 11. Vorticity stream function $\Psi_{75}$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=20$.

Figure 12. Vorticity stream function $\Psi_{75}$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=30$.

Figure 13. Vorticity stream function $\Psi_{75}$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=40$.

Figure 14. Vorticity stream function $\Psi_{75}$ of
$\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=50$.

Figure 15. Vorticity stream function $\Psi_{75}$ of $\boldsymbol{\mathrm{Vort}}\,{\bf U}_{75}$ at $t=60$. |
---|---
Energy spectrum $E(L)$ of ${\bf U}_{75}(t)$. | Energy spectrum $E(L)$ of $\widetilde{\Phi}_{150}({\bf U}_{75}(t))$.
Figure 16. Energy spectra of velocity ${\bf U}_{75}(t)$ and $\Phi_{150}({\bf U}_{75}(t))$. |
---|---
Energy spectrum $L^{4}*E(L)$ of ${\bf U}_{75}(t)$ | Energy spectrum $L^{4}*E(L)$ of $\widetilde{\Phi}_{150}({\bf U}_{75}(t))$.
Figure 17. Scaled energy spectra $L^{4}*E(L)$ of ${\bf U}_{75}(t)$ and
$\widetilde{\Phi}_{150}({\bf U}_{75}(t))$.

##  6\. Appendix

In this section, we generalize certain known domain case estimates, that are
fundamental for the NSE analysis, to the spherical surface case and hence
prove (2.30). Following [16, Page 574], for ${\bf u}\in C^{\infty}(TS)$, we
extend ${\bf u}$ to the spherical layer $S\times I$,
$I=(r_{1},r_{2}),0<r_{1}<1<r_{2}<\infty$ by the formula

(6.1) $\widetilde{{\bf u}}(\widehat{{\bf x}})=\varphi(|\widehat{{\bf x}}|){\bf
u}(\widehat{{\bf x}}/|\widehat{{\bf x}}|),$

where $\varphi(t)\in C_{0}^{\infty}(I)$, $\varphi(t)\geq 0,~{}t\in I$, and
$\varphi(1)=1$. We have

$\int_{S\times I}|{\widetilde{{\bf
u}}}|^{p}=\int_{r_{1}}^{r_{2}}\varphi^{p}(r)\int_{S}|{\bf u}|^{p}dS$

In other words

(6.2) $\|{\widetilde{{\bf u}}}\|_{L^{p}(S\times I)}=c\|{\bf u}\|_{L^{p}(TS)},$

where $c=c(\varphi,r_{1},r_{2})$. Suppose ${\bf u},{\bf v},{\bf w}$ are
extended from $S$ to the spherical layer $S\times I$ by (6.1). Then from [16,
Lemma 4.3] we have

(6.3) $b({\bf u},{\bf v},{\bf w})=cb(\widetilde{{\bf u}},\widetilde{{\bf
v}},\widetilde{{\bf w}}),$

On the sphere $S$, we have the following version of Sobolev embedding
inequality [1]

(6.4) $\|{\bf u}\|_{L^{q}(TS)}\leq C\|{\bf u}\|_{H^{s}(TS)},\qquad
s<1,\quad\frac{1}{q}=\frac{1}{2}-\frac{s}{2}.$

The following nonlinearity estimate is an adaptation of [6, Proposition 6.1]
for $S$.

######  Proposition 6.1.

Let $s_{1},s_{2},s_{3}\geq 0$ be real numbers, and we assume that
$s_{1}+s_{2}+s_{3}\geq 1$ and
$(s_{1},s_{2},s_{3})\neq(0,0,1),(0,1,0),(1,0,0)$. Then there exists a constant
depending on $s_{1},s_{2},s_{3}$ such that

$|b({\bf u},{\bf v},{\bf w})|\leq C\|{\bf u}\|_{H^{s_{1}}(TS)}\|{\bf
v}\|_{H^{s_{2}+1}(TS)}\|{\bf w}\|_{H^{s_{3}}(TS)}$

or in the extrapolated form, writing $H^{s}(TS)=H^{s}_{TS}$, for all ${\bf
u},{\bf v},{\bf w}\in C^{\infty}(TS)$

$\displaystyle|b({\bf u},{\bf v},{\bf w})|\leq C$ $\displaystyle\|{\bf
u}\|^{1+[s_{1}]-s_{1}}_{H^{[s_{1}]}_{TS}}\|{\bf
u}\|^{s_{1}-[s_{1}]}_{H^{[s_{1}]+1}_{TS}}\|{\bf
v}\|^{1+[s_{2}]-s_{2}}_{H^{[s_{2}]+1}_{TS}}\|{\bf
v}\|^{s_{2}-[s_{2}]}_{H^{[s_{2}]+2}_{TS}}\|{\bf
w}\|^{1+[s_{3}]-s_{3}}_{H^{[s_{3}]}_{TS}}\|{\bf
w}\|^{s_{3}-[s_{3}]}_{H^{[s_{3}]+1}_{TS}}.$

###### Proof.

Let ${\bf u},{\bf v},{\bf w}\in C^{\infty}(TS)$ and $\widetilde{{\bf
u}},\widetilde{{\bf v}},\widetilde{{\bf w}}$ be their corresponding extension
to the spherical layer $\widetilde{\Omega}:=S\times I$. Let us consider first
the case $s_{i}<1$ for $i=1,2,3$. Define associated constants
$q_{1},q_{2},q_{2},q_{4}$ so that $\sum_{i=1}^{4}\frac{1}{q_{i}}=1$ and
$\frac{1}{q_{i}}=\frac{1}{2}-\frac{s_{i}}{2}$ for $i=1,2,3$. Then, by Hölder’s
inequality

$\displaystyle|b({\widetilde{{\bf u}}},{\widetilde{{\bf v}}},{\widetilde{{\bf
w}}})|$ $\displaystyle=$
$\displaystyle\left|\sum_{i,j=1}^{3}\int_{\widetilde{\Omega}}\widetilde{u}_{j}\frac{\partial\widetilde{v}_{i}}{\partial
x_{j}}\widetilde{w}_{i}1\right|\leq\|{\widetilde{{\bf
u}}}\|_{L^{q_{1}}(\Omega)}\|\nabla{\widetilde{{\bf
v}}}\|_{L^{q_{2}}(\Omega)}\|{\widetilde{{\bf
w}}}\|_{L^{q_{3}}(\Omega)}\|1\|_{L^{q_{4}}(\Omega)}$

Restricting to the sphere by (6.3) and (6.2) then using the Sobolev embedding
theorem on the sphere (6.4) we have

$\displaystyle|b({\bf u},{\bf v},{\bf w})|$ $\displaystyle\leq$ $\displaystyle
C\|{\bf u}\|_{L^{q_{1}}_{TS}}\|\boldsymbol{\mathrm{Grad}}\,{\bf
v}\|_{L^{q_{2}}_{TS}}\|{\bf w}\|_{L^{q_{3}}_{TS}}\leq C\|{\bf
u}\|_{H^{s_{1}}_{TS}}\|\boldsymbol{\mathrm{Grad}}\,{\bf
v}\|_{H^{s_{2}}_{TS}}\|{\bf w}\|_{H^{s_{3}}_{TS}}.$

∎

######  Lemma 6.1.

Let $\delta\in(1/2,1)$ be given and ${\bf u},{\bf v}\in V$. Then there exists
$C$, independent of ${\bf u}$ and ${\bf v}$, such that

$\|\boldsymbol{\mathrm{A}}^{-\delta}\boldsymbol{\mathrm{B}}({\bf u},{\bf
v})\|\leq C\left\\{\begin{array}[]{cc}\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf
u}\|\|{\bf v}\|&\\\ \|{\bf u}\|\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf
v}\|&\\\ \end{array},\right.\qquad{\bf u},{\bf v}\in V.$

###### Proof.

Let ${\bf u},{\bf v},{\bf w}\in V$. Hence from (2.26), and by using
Proposition 6.1 with $s_{1}=0$, $s_{2}=2\delta-1>0$ and $s_{3}=2-2\delta>0$,

$|b({\bf u},{\bf v},\boldsymbol{\mathrm{A}}^{-\delta}{\bf w})|=|b({\bf
u},\boldsymbol{\mathrm{A}}^{-\delta}{\bf w},{\bf v})|\leq C\|{\bf
u}\|\|\boldsymbol{\mathrm{A}}^{-\delta}{\bf w}\|_{H^{2\delta}(TS)}\|{\bf
v}\|_{H^{2-2\delta}(TS)}.$

Since $\|\boldsymbol{\mathrm{A}}^{-\delta}{\bf
w}\|_{H^{2\delta}(TS)}=\|\boldsymbol{\mathrm{A}}^{\delta}\boldsymbol{\mathrm{A}}^{-\delta}{\bf
w}\|=\|{\bf w}\|$, we get

$|b({\bf u},{\bf v},\boldsymbol{\mathrm{A}}^{-\delta}{\bf w})|\leq C\|{\bf
u}\|\|{\bf w}\|\|{\bf v}\|_{H^{2-2\delta}(TS)},\qquad{\bf w}\in V.$

Since the inequality is true for all ${\bf w}\in V$, we obtain the first bound

$\|\boldsymbol{\mathrm{A}}^{-\delta}\boldsymbol{\mathrm{B}}({\bf u},{\bf
v})\|\leq C\|{\bf u}\|\|{\bf v}\|_{H^{2-2\delta}(TS)}=C\|{\bf
u}\|\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf v}\|.$

To obtain the second bound, we again use (2.26) and Proposition 6.1 but with
$s_{1}=2-2\delta$, $s_{2}=2\delta-1$ and $s_{3}=0$,

$\displaystyle|b({\bf u},{\bf v},\boldsymbol{\mathrm{A}}^{-\delta}{\bf
w})|=|b({\bf u},\boldsymbol{\mathrm{A}}^{-\delta}{\bf w},{\bf v})|$
$\displaystyle\leq$ $\displaystyle C\|{\bf
u}\|_{H^{2-2\delta}(TS)}\|\boldsymbol{\mathrm{A}}^{-\delta}{\bf
w}\|_{H^{2\delta}(TS)}\|{\bf v}\|$ $\displaystyle=$ $\displaystyle
C\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf u}\|\|{\bf w}\|\|{\bf
v}\|,\qquad{\bf w}\in V.$

Since the inequality is true for all ${\bf w}\in H$, we obtain

$\|\boldsymbol{\mathrm{A}}^{-\delta}\boldsymbol{\mathrm{B}}({\bf u},{\bf
v})\|\leq C\|\boldsymbol{\mathrm{A}}^{1-\delta}{\bf u}\|\|{\bf v}\|.$

∎

#### Acknowledgments:

The support of the Australian Research Council under its Discovery and Centre
of Excellence programs is gratefully acknowledged. The authors thank
Professors M. Farge and E. S. Titi [29] for valuable discussions.

## References

  * [1] T. Aubin, Nonlinear analysis on manifolds, Monge-Amp$\grave{e}$re Equations. Springer-Verlag, 1982.
  * [2] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Dover, 2001.
  * [3] M.E. Brachet, M. Meneguzzi, H. Politano, P.L. Sulem, The dynamics of freely decaying two-dimensional turbulence. J. Fluid Mech. 194 (1988), 333-349.
  * [4] M.A.J. Chaplain, M. Ganesh, I.G. Graham, Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth. J. Math. Biology 42 (2001), 387-423.
  * [5] C. Cao, M.A. Rammaha, E.S. Titi, The Navier-Stokes equations on the rotating 2-D sphere: Gevrey regularity and asymptotic degrees of freedom. Zeit. Ang. Math. Phys. 50 (1999), 341–360.
  * [6] P. Constantin, C. Foias, Navier-Stokes equations, Chicago Lectures in Mathematics, The University of Chicago Press, 1988.
  * [7] A. Doelman, E.S. Titi, Regularity of solutions and the convergence of the Galerkin method in the Ginzburg-Landau equation. Numer. Funct. Anal. Optim. 14 (1993), 299–321.
  * [8] A. Debussche, T. Dubois, R. Temam, The nonlinear Galerkin method: a multiscale method applied to the simulation of homogeneous turbulent flows. Theor. Comp. Fluid Dyn. 7 (1995), 279–315.
  * [9] C. Foias, O. Manley, R. Rosa, R. Temam, Navier-Stokes Equations and Turbulence. Cambridge University Press, 2001.
  * [10] C. Foias, M.S. Jolly, I.G. Kevrekidis, G.R. Sell, E.S. Titi, On the computation of inertial manifolds. Phys. Lett. A 131 (1988), 433-436.
  * [11] C. Foias, O. Manley, R. Temam, Modelling of the interaction of small and large eddies in two dimensional turbulence flows, RAIRO Modél. Math. Anal. Numér. 22 (1988), 93–118.
  * [12] M.J. Fengler, W. Freeden, A nonlinear Galerkin scheme involving vector and tensor spherical harmonics for solving the incompressible Navier-Stokes equation on the sphere. SIAM J. Sci. Comp. 27 (2005), 967–994.
  * [13] M. Ganesh, K. Mustapha A fully discrete $H^{1}$-Galerkin method with quadrature for nonlinear advection-diffusion-reaction equations. Numer. Algorithms, 43 (2006), 355–383.
  * [14] E. Hairer, G. Wanner Solving Ordinary Differential Equations II: Stiff and differential-algebraic problems. Springer-Verlag, 1982.
  * [15] J. Hesthaven, S. Gottlieb, D. Gottlieb, Spectral Methods for Time-Dependent Problems. Cambridge University Press, 2007.
  * [16] A.A. Il’in, The Navier-Stokes and Euler equations on two dimensional closed manifolds. Math. USSR Sbornik 69 (1991), 559–579.
  * [17] A.A. Il’in, Partially dissipative semigroups generated by the Navier-Stokes system on two-dimensional manifolds, and their attractors. Russian Acad. Sci. Sbornik Mathematics 78 (1994), 47–76.
  * [18] A.A. Il’in, Navier-Stokes equations on the rotating sphere. A simple proof of the attractor dimension estimate. Nonlinearity 7 (1994), 31–39.
  * [19] A.A. Il’in, A.N. Filatov, On unique solvability of the Navier-Stokes equations on the two dimensional sphere. Soviet Math. Dokl. 38 (1989), 9–13.
  * [20] D.A. Jones, E.S. Titi, A remark on quasi-stationary approximate inertial manifold for the Navier-Stokes equations. SIAM J. Math. Anal. 25 (1994),894–914.
  * [21] S. Kunis, D. Potts, Fast spherical Fourier algorithms. J. Comp. Appl. Math. 161 (2003), 75–98.
  * [22] J. Norbury, I. Roulstone (Eds.), Large-Scale Atmosphere-Ocean Dynamics, Vol. 1& 2, Cambridge University Press, 2002.
  * [23] S. A. Orszag, Turbulence and transition: A progress report. Proc. 5th Int. Conf. Numer. Meth. Fluids Dyn., Lecture Notes in Phys., Vol. 59, Springer, 32–51, 1977.
  * [24] M. Pieper, Vector hyperinterpolation on the sphere, J. Approx. Theory. 156 (2009), 173–186.
  * [25] L. F. Shampine, M. W. Reichlet The Matlab ODE suite. SIAM J. Sci. Comput., 18 (1997), 1–22.
  * [26] T. Tao, Why global regularity for Navier-Stokes is hard. http://terrytao.wordpress.com/2007/03/18/
why-global-regularity-for-navier-stokes-is-hard/.

  * [27] R. Temam, Navier-Stokes Equations, Theory and Numerical Analysis, Amer. Math. Soc., 2001.
  * [28] R. Temam, S. Wang, Inertial forms of Navier-Stokes equations on the sphere. J. Funct. Anal. 117 (1993), 215–242.
  * [29] E. S. Titi, Private communication, (2009).
  * [30] D. A. Varshalovich, A. N. Moskalev, V. K. Khersonskii, Quantum Theory of Angular Momentum, World Scientific, 1988.