# An analytically divergence-free collocation method for the incompressible
Navier-Stokes equations on the rotating sphere

Tino Franz  Applied and Numerical Analysis, University of Bayreuth, D-95440
Bayreuth, Germany<EMAIL_ADDRESS>Orcid-ID:0000-0002-4273-7639

(Date: August 27, 2024)

###### Abstract.

In this work, we develop a high-order collocation method using radial basis
function (RBF) for the incompressible Navier-Stokes equation (NSE) on the
rotating sphere. The method is based on solving the projection of the NSE on
the space of divergence-free functions. For that, we use matrix valued kernel
functions which allow an analytically divergence-free approximation of the
velocity field. Using kernel functions which lead to rotation-free
approximations, the pressure can be recovered by a simple kernel exchange in
one of the occurring approximations, without solving an additional Poisson
problem. We establish precise error estimates for the velocity and the
pressure functions for the semi-discretised solution. In the end, we give a
short estimate of the numerical cost and apply the new method to an
experimental test case.

###### Key words and phrases:

collocation, radial basis function, meshfree methods, Navier-Stokes equation,
vector spherical harmonics

###### MSC Classification Mathematics Subject Classification:

65M12, 65M70, 76D05

The author would like to thank Holger Wendland for his encouragement and
several valuable discussions which lead to this paper.

##  1\. Introduction

In this paper we will develop and analyse a discretisation method for the
incompressible Navier-Stokes equations (NSE) on the rotating unit sphere using
a kernel-based divergence-free approximation method. This method is highly
motivated by and can be seen as an extension of the work of Keim [20], who
constructed a collocation scheme for the time-depended Navier-Stokes equations
on the $d$-dimension torus. See also [21] for the method applied to the time-
dependent Stokes equation.

Kernel-based divergence-free approximation methods have been studied and
analysed in [7, 3, 11, 14, 12, 22, 23, 8], and were successfully applied to
different problems, see for example, [27, 28, 34]. In [25] and [13], the
authors extended and analysed the theory of divergence-free and rotation-free
kernel-based approximations for tangential vector fields on manifolds,
especially the unit sphere. These approximation methods will be used here.

The two dimensional sphere in ${\mathbb{R}}^{3}$ is given by
${\mathbb{S}}^{2}=\\{{{\boldsymbol{x}}}\in{\mathbb{R}}^{3}\,|\,\|{{\boldsymbol{x}}}\|_{2}=1\\}$,
where $\|\cdot\|_{2}$ denotes the standard euclidean norm. A vector field
${{\boldsymbol{u}}}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$ is called to be
tangential, if
${{\boldsymbol{u}}}({{\boldsymbol{x}}})\cdot{{\boldsymbol{x}}}=0$ for all
${{\boldsymbol{x}}}\in{\mathbb{S}}^{2}$. The Navier-Stokes equations on the
sphere are given as follows, see [19, 5]. Given a tangential and divergence-
free vector field ${{\boldsymbol{u}}}_{0}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$
and a tangential and divergence-free forcing term
${{\boldsymbol{f}}}:[0,T)\times{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$, we seek
the solution
${{\boldsymbol{u}}}:[0,T)\times{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$, which
represents the velocity, and $p:[0,T)\times{\mathbb{S}}^{2}\to{\mathbb{R}}$,
which represents the pressure, of

(1)
$\displaystyle\partial_{t}{{\boldsymbol{u}}}-\nu{\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}}+\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})+\mathcal{C}({{\boldsymbol{u}}})+\nabla^{*}p$
$\displaystyle={{\boldsymbol{f}}}$ $\displaystyle\text{on
}{\mathbb{S}}^{2}\times(0,T),$ (2)
$\displaystyle\nabla^{*}\cdot{{\boldsymbol{u}}}$ $\displaystyle=0$
$\displaystyle\text{on }{\mathbb{S}}^{2}\times(0,T),$ (3)
$\displaystyle{{\boldsymbol{u}}}(\cdot,0)$
$\displaystyle={{\boldsymbol{u}}}_{0}$ $\displaystyle\text{on
}{\mathbb{S}}^{2},$

where

$\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}}):=({{\boldsymbol{u}}}\cdot\nabla_{*}){{\boldsymbol{v}}},\quad\mathcal{C}({{\boldsymbol{u}}}):=l{{\boldsymbol{n}}}\times{{\boldsymbol{u}}},$

$\nu>0$ is the diffusion coefficient, ${{\boldsymbol{n}}}$ is the normal
vector and $l({{\boldsymbol{x}}})=2\Omega x_{3}=2\Omega\sin(\phi)$ is the
Coriolis coefficient with angular velocity $\Omega$ and latitude
$\phi=\arcsin(x_{3})$. $\mathcal{B}$ represents the non-linear convection
operator, while $\mathcal{C}$ is the Coriolis operator caused by the rotation
of the sphere. A proper definition of the differential operators on the sphere
will be given in Section 2.1.

To simplify this system, a common approach is to project equation (1) on its
divergence-free part using the Leray operator $\mathcal{P}$ to eliminate the
gradient of the pressure, which leads to

(4)
$\partial_{t}{{\boldsymbol{u}}}-\nu{\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}}=\mathcal{P}({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{C}({{\boldsymbol{u}}})),\quad\text{in
}(0,T)\times{\mathbb{S}}^{2},$

where ${{\boldsymbol{u}}}$ still satisfies (2).

Our method is using a kernel-based interpolation method, which automatically
gives an approximation of the divergence-free as well as the rotation-free
part of a function. More precisely, given a set of points
$\mathcal{X}=\\{{{\boldsymbol{x}}}_{1},\ldots\,{{\boldsymbol{x}}}_{N}\\}\subset{\mathbb{S}}^{2}$
and data sites
${{\boldsymbol{f}}}_{j}={{\boldsymbol{f}}}({{\boldsymbol{x}}}_{j})\in{\mathbb{R}}^{3}$,
$1\leq j\leq N$, with ${{\boldsymbol{f}}}$ being a tangential vector field,
the interpolation of ${{\boldsymbol{f}}}$ is given by

(5)
${{\boldsymbol{f}}}_{h}({{\boldsymbol{x}}})=\sum_{j=1}^{N}\Phi({{\boldsymbol{x}}},{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j},\quad{{\boldsymbol{x}}}\in{\mathbb{S}}^{3},$

where ${\boldsymbol{\beta}}_{j}\in{\mathbb{R}}^{3}$ with
${\boldsymbol{\beta}}_{j}\cdot{{\boldsymbol{x}}}_{j}=0$ for every $1\leq j\leq
N$ and the kernel
$\Phi:{\mathbb{S}}^{2}\times{\mathbb{S}}^{2}\to{\mathbb{R}}^{3\times 3}$ is
given such that the interpolation ${{\boldsymbol{f}}}_{h}$ is tangential. In
our context, the kernel will be given such that $\Phi=\Phi_{div}+\Phi_{curl}$,
i.e. the kernel function is the sum of a kernel function which yields
divergence-free and one that yields curl-free interpolations. Hence, by simply
exchanging the kernel in (5) by $\Phi_{div}$, we obtain an analytically
divergence-free approximation to the divergence-free part of
${{\boldsymbol{f}}}$, which leads to an approximation $\mathcal{P}_{h}$ of the
Leray operator. See [25] and [13] for more details.

To derive an approximation of the divergence-free velocity, we use the method
of lines to separate the variables, such that our approximation has the form

(6)
${{\boldsymbol{u}}}_{h}({{\boldsymbol{x}}},t)=\sum_{j=1}^{N}\Phi_{div}({{\boldsymbol{x}}},{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}(t),\quad{{\boldsymbol{x}}}\in{\mathbb{S}}^{3},t\in(0,T).$

This ensures that (2) is automatically satisfied. Equation (5) is then to be
enforced pointwise on $\mathcal{X}$, which leads to the semi-discrete system
of ODEs

(7)
$\partial_{t}{{\boldsymbol{u}}}_{h}({{\boldsymbol{x}}}_{j},t)-\nu{\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}}_{h}({{\boldsymbol{x}}}_{j},t)=\mathcal{P}_{h}({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h}))({{\boldsymbol{x}}}_{j},t)\
$

for ${{\boldsymbol{x}}}_{j}\in\mathcal{X}$ and $t\in(0,T).$ The main part of
the work consists of analysing the convergence of the solution of (7) to the
solution of (4). Here, the right hand side is given by using the approximated
Leray operator $\mathcal{P}_{h}$, which gives a divergence-free approximation
of the form

(8)
$\mathcal{P}_{h}({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h}))({{\boldsymbol{x}}}_{j},t)=\sum_{k=1}^{N}\Phi_{div}({{\boldsymbol{x}}}_{k},{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}(t)$

for ${{\boldsymbol{x}}}_{j}\in\mathcal{X}$ and $t\in(0,T)$, where the
coefficients ${\boldsymbol{\beta}}_{j}\in{\mathbb{R}}^{3}$ are obtained by
interpolating
${{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})$
on the supporting points $\mathcal{X}$.

Using standard time discretisation methods such as Runge-Kutta methods, the
system given in (7) can be turned into a fully discretised system. Hence, the
linear system to determine ${\boldsymbol{\beta}}_{j}$ for (8) as well as the
linear system to determine the ${\boldsymbol{\alpha}}_{j}$ in (6) have to be
solved in every time-step. However, choosing kernel functions $\Phi$ with
compact support, the matrices in both linear systems can be made sparse, such
that both linear systems can be solved sufficiently fast. As we will
additionally point out, the efficiency can be theoretically further increased
by using Lagrange basis functions.

This method has several advantages. It provides analytically divergent-free
approximations to the velocity, which do not need any spatial numerical
integration for their calculation. Moreover, it can be set up to construct
arbitrarily smooth solution just by the choice of the kernel function, which
leads to high-order spatial approximations.

Finally, it leads easily to a high-order approximation of the pressure just by
exchanging the kernel function in (8), i.e.

$p_{h}({{\boldsymbol{x}}},t)=\sum_{j=1}^{N}\Psi({{\boldsymbol{x}}},{{\boldsymbol{x}}}_{j})^{T}{\boldsymbol{\beta}}_{j}(t).$

for a specific kernel function
$\Psi:{\mathbb{S}}^{2}\times{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$. Hence, we do
not need to solve an additional Poisson equation to obtain an approximation of
the pressure. However, we have to assume a strong solution to the Navier-
Stokes equations to ensure the existence of the point evaluations for the
collocation method.

Even if there is a lot of work for numerical methods for the Navier-Stokes
equations on the Euclidean space or bounded Euclidean domains, see for example
[6, 31], the are only few results about numerical methods and their analysis
for the Navier-Stokes equations on bounded manifolds and the sphere. First of
all, there is the paper of Fengler and Freeden [9], who developed a spectral
non-linear Galerkin method which is based on vector and tensor spherical
harmonics. However, they mainly concentrate on the implementation of their
method, and did not provide a proper mathematical analysis. Next, there is the
paper of Ganesh, Le Gia and Sloan [16], who developed a pseudospectral
quadrature method for the Navier-Stokes equations. Based on the Gevrey-
regularity, they were able to show spectral convergence of their method.
However, their method needs numerical integration of the occurring integrals.
Moreover, both method do not provide a fast and easy possibility to calculate
an approximation of the pressure.

This paper is organized as follows. In the second section we will give a short
introduction to spherical harmonics and Sobolev spaces on the sphere. The
approximation tools we will need for constructing the our approximation method
will be given in the third section. In the fourth section we will shortly
introduce the problem we want to solve, together with a few theoretical
results concerning the non-linear part of the Navier-Stokes equations. After
that, in the fifth section, we will construct our approximation scheme and
give the proof of convergence as the main result of this article. Finally, we
demonstrate the new method numerically.

##  2\. Notation and preliminaries

In this section, we want to provide the basic theory to spherical harmonics
and function spaces on the sphere. We will start with basic notations.

###  2.1. Basic notations and differential operators on the sphere

We will differ between scalar-valued functions
$f:{\mathbb{S}}^{2}\to{\mathbb{R}}$ and vector fields
${{\boldsymbol{u}}}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$, which will be
written in bold letters. Hence, we will also denote spaces, which contain
vector-valued functions, with bold letters. For a
${{\boldsymbol{x}}}\in{\mathbb{S}}^{2}$, the set
$T_{{\boldsymbol{x}}}{\mathbb{S}}^{2}$ contains all
${\boldsymbol{\alpha}}\in{\mathbb{R}}^{3}$ which are tangent to
${\mathbb{S}}^{2}$ at ${{\boldsymbol{x}}}$, i.e.
${\boldsymbol{\alpha}}\cdot{{\boldsymbol{x}}}=0$.

We define the surface gradient $\nabla^{*}$, the surface curl-operator
$\boldsymbol{L}^{*}$ and the surface Laplace-Beltrami operator $\Delta^{*}$
for a sufficiently smooth scalar-valued function
$f:{\mathbb{S}}^{2}\to{\mathbb{R}}$ by

$\displaystyle\nabla^{*}f$ $\displaystyle:=\nabla f-\left<\nabla
f,{\boldsymbol{\nu}}\right>{\boldsymbol{\nu}},$
$\displaystyle{{\boldsymbol{L}}}^{*}f({{\boldsymbol{x}}})$
$\displaystyle:={{\boldsymbol{x}}}\times\nabla^{*}f({{\boldsymbol{x}}}),$
$\displaystyle\Delta^{*}f({{\boldsymbol{x}}})$ $\displaystyle:=\Delta
f({{\boldsymbol{x}}})-\left<{{\boldsymbol{x}}},H(f){{\boldsymbol{x}}}\right>-2\left<\nabla
f({{\boldsymbol{x}}}),{{\boldsymbol{x}}}\right>,$

where $\nabla$ and $\Delta$ are the standard gradient and Laplace operators,
respectively, $\times$ denotes the vector product and $H$ denotes the Hessian
Matrix $H(f)=(\partial_{ij}f)$. For $1\leq j\leq 3$, the $j$-th component of
$\nabla^{*}$ will be denoted by $\partial^{*}_{j}$. The Laplace-Beltrami
operator can also be written as

$\Delta^{*}=\nabla^{*}\cdot\nabla^{*}={{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{L}}}^{*}.$

For vector fields ${{\boldsymbol{u}}}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$, we
define the surface-divergence $\nabla^{*}\cdot$, the vectorial curl-operator
${{\boldsymbol{L}}}^{*}\cdot$ and the vectorial Laplace-Beltrami operator
${\boldsymbol{\Delta}}^{*}$ as

$\displaystyle\nabla^{*}\cdot{{\boldsymbol{u}}}$
$\displaystyle:=\sum_{i=1}^{3}(\nabla^{*}u_{i})\cdot{{\boldsymbol{e}}}_{i},$
$\displaystyle{{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{u}}}({{\boldsymbol{x}}})$
$\displaystyle:=\sum_{i=1}^{3}({{\boldsymbol{L}}}^{*}u_{i})\cdot{{\boldsymbol{e}}}_{i},$
$\displaystyle{\boldsymbol{\Delta}}^{*}$
$\displaystyle:=p_{nor}(\Delta^{*}+2)p_{nor}+p_{tan}\Delta^{*}p_{tan},$

where ${{\boldsymbol{e}}}_{i}$ denotes the $i$-th unit vector,
$p_{nor}({{\boldsymbol{u}}})({{\boldsymbol{x}}})=\langle{{\boldsymbol{u}}}({{\boldsymbol{x}}}),{{\boldsymbol{x}}}\rangle{{\boldsymbol{x}}}$
and $p_{tan}({{\boldsymbol{u}}})={{\boldsymbol{u}}}-p_{nor}{{\boldsymbol{u}}}$
are the normal and tangential part of the function ${{\boldsymbol{u}}}$ with
respect to the sphere, respectively. For further details on differential
operators on the sphere, see for example [10]. Note that we will only be
interested in tangential vector fields ${{\boldsymbol{u}}}$, where the normal
part vanishes. In this case, the vectorial Laplace-Beltrami operator
degenerates to the scalar one applied to the single components $u_{i}$, $1\leq
i\leq 3$, of ${{\boldsymbol{u}}}$. Finally, we will call a vector field
${{\boldsymbol{u}}}$ to be divergence-free if
$\nabla^{*}\cdot{{\boldsymbol{u}}}$ vanishes and curl-free if the rotation
${{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{u}}}$ vanishes.

We will use time-depended function spaces, which are defined as follows. For
an arbitrary Banach space $X$, the space $C([0,T];\\!\\!X)$ consists of all
continuous functions from $[0,T]$ to $X$. For $1\leq p\leq\infty$, we denote
the space of all strongly measurable functions from $[0,T]$ to $X$ such that
$t\mapsto\|f(t)\|_{X}\in L^{p}(0,T)$ by $L^{p}(0,T;\\!\\!X)$, with norm

$\|f\|_{L^{p}(0,T;\\!\\!X)}=\begin{cases}\left(\int_{0}^{T}\|f(t)\|_{X}^{2}dt\right)^{1/p},&1\leq
p<\infty,\\\ \operatorname*{ess\,sup}_{0\leq t\leq
T}\|f(t)\|_{X},&p=\infty.\end{cases}$

For an ${{\boldsymbol{u}}}:{\mathbb{S}}^{2}\times[0,T]\to{\mathbb{R}}^{3}$, we
will sometimes write ${{\boldsymbol{u}}}(t)$ instead of
${{\boldsymbol{u}}}(\cdot,t)$ for the sake of readability.

###  2.2. Spherical Harmonics and the space of square integrable functions

As usual, the Hilbert space of square integrable functions on the sphere is
given by

$L^{2}:=L^{2}({\mathbb{S}}^{2})=\left\\{f:{\mathbb{S}}^{2}\to{\mathbb{R}}\,\Big{|}\,\int_{{\mathbb{S}}^{2}}|f({{\boldsymbol{x}}})|^{2}dS({{\boldsymbol{x}}})<\infty\right\\}$

where the inner product is given by

$\langle
f,g\rangle:=\int_{S^{2}}f({{\boldsymbol{x}}})g({{\boldsymbol{x}}})dS({{\boldsymbol{x}}}),\quad
f,g\in L^{2}.$

Another way to describe $L^{2}$ functions on the sphere is by using spherical
harmonics, see [24, 10] for an introduction. The spherical harmonics
$Y_{\ell,k}$ of order $\ell\in{\mathbb{N}}_{0}$ and $1\leq k\leq 2\ell+1$ are
the eigenfunctions of $-\Delta_{*}$ with corresponding eigenvalues

$\lambda_{\ell}=\ell(\ell+1),\quad\ell\in{\mathbb{N}}_{0}.$

They form an orthonormal basis of $L^{2}$, which means that every function
$f\in L^{2}$ has a Fourier representation

$f({{\boldsymbol{x}}})=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2n+1}\widehat{f}_{\ell,k}Y_{\ell,k}({{\boldsymbol{x}}}),\quad\widehat{f}_{\ell,k}=\langle
f,Y_{\ell,k}\rangle.$

Using Parseval’s identity, the $L^{2}$ inner product can then also be written
in terms of the Fourier coefficients by

$\langle
f,g\rangle=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2n+1}\widehat{f}_{\ell,k}\widehat{g}_{\ell,k},\quad
f,g\in L^{2}.$

The idea can be extended to vector-valued functions. We denote the space of
square-integrable vector-valued functions which are tangential to the sphere
by

${{\boldsymbol{L}}}^{2}:={{\boldsymbol{L}}}^{2}({\mathbb{S}}^{2})=\left\\{{{\boldsymbol{u}}}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}\,\Big{|}\,{{\boldsymbol{u}}}\text{
tangential and
}\int_{{\mathbb{S}}^{2}}\|{{\boldsymbol{u}}}({{\boldsymbol{x}}})\|_{2}^{2}dS({{\boldsymbol{x}}})<\infty\right\\}$

where the inner product is given by

$\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle:=\int_{S^{2}}{{\boldsymbol{u}}}({{\boldsymbol{x}}})\cdot{{\boldsymbol{v}}}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}}),\quad{{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}.$

Again, we can build an orthonormal basis for this function space which is
given by the analogue of the spherical harmonics for vector fields: the
vectorial spherical harmonics. There are three different,
${{\boldsymbol{L}}}_{2}$-orthogonal types: the divergence-free and the curl-
free type, which are both tangential to the sphere, and one type that is
normal to the sphere. Here, we are only interested in the tangential type.

For $\ell\geq 1$ and $1\leq k\leq 2\ell+1$ the divergence-free vector
spherical harmonics are given by
${{\boldsymbol{y}}}_{\ell,k}={{\boldsymbol{L}}}_{*}Y_{\ell,k}/\sqrt{\lambda_{\ell}}$
and the curl-free vector spherical harmonics are given by
${{\boldsymbol{z}}}_{\ell,k}=\nabla^{*}Y_{\ell,k}/\sqrt{\lambda_{\ell}}$ for
all . Together, they form an orthonormal basis for ${{\boldsymbol{L}}}^{2}$
with respect to the ${{\boldsymbol{L}}}^{2}$ norm. Moreover,
${{\boldsymbol{y}}}_{\ell,k}$ and ${{\boldsymbol{z}}}_{\ell,k}$ are the
eigenfunctions to the vectorial Laplace Beltrami operator with corresponding
eigenvalue $\lambda_{\ell}$.

Again, we can write each vector field
${{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}$ in its Fourier form

${{\boldsymbol{u}}}=\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\left(\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}{{\boldsymbol{y}}}_{\ell,k}+\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}{{\boldsymbol{z}}}_{\ell,k}\right),$

where the divergence-free and the curl-free Fourier coefficients are given by

$\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}=\langle{{\boldsymbol{u}}},{{\boldsymbol{y}}}_{\ell,k}\rangle,\quad\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}=\langle{{\boldsymbol{u}}},{{\boldsymbol{z}}}_{\ell,k}\rangle.$

The ${{\boldsymbol{L}}}^{2}$ inner product can again be written as
$\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle=\sum_{\ell=1}^{\infty}\sum_{k=1}^{2n+1}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{div}_{\ell,k}+\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{curl}_{\ell,k}$
for ${{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}$.

The Fourier expansion gives us a direct and unique decomposition of a
tangential vector field ${{\boldsymbol{u}}}$ in a divergence-free part
${{\boldsymbol{u}}}_{div}=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}{{\boldsymbol{y}}}_{\ell,k}$
and a curl-free part
${{\boldsymbol{u}}}_{curl}=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}{{\boldsymbol{z}}}_{\ell,k}$.
The corresponding Hilbert spaces for square-integrable divergence-free and the
curl-free functions are therefore given by

$\begin{split}{{\boldsymbol{L}}}^{2}_{div}&:=\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}\,|\,\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}=0,\;1\leq\ell\leq\infty,1\leq
k\leq 2\ell+1\right\\},\\\
{{\boldsymbol{L}}}^{2}_{curl}&:=\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}\,|\,\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}=0,\;1\leq\ell\leq\infty,1\leq
k\leq 2\ell+1\right\\}\end{split}$

with given inner products

$\begin{split}\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{L}}}^{2}_{div}}&:=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{div}_{\ell,k}=\int_{{\mathbb{S}}^{2}}{{\boldsymbol{u}}}_{div}({{\boldsymbol{x}}}){{\boldsymbol{v}}}_{div}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}}),\quad{{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}_{div},\\\
\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{L}}}^{2}_{curl}}&:=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{curl}_{\ell,k}=\int_{{\mathbb{S}}^{2}}{{\boldsymbol{u}}}_{curl}({{\boldsymbol{x}}}){{\boldsymbol{v}}}_{curl}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}}),\quad{{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}_{curl}.\end{split}$

for ${{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}$. The
spaces ${{\boldsymbol{L}}}^{2}_{div}$ and ${{\boldsymbol{L}}}^{2}_{curl}$ are
a direct decomposition of ${{\boldsymbol{L}}}^{2}$, i.e.
${{\boldsymbol{L}}}^{2}={{\boldsymbol{L}}}^{2}_{div}\oplus{{\boldsymbol{L}}}^{2}_{curl}$.
Note that for sufficiently smooth
${{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}_{div}$ and
${{\boldsymbol{v}}}\in{{\boldsymbol{L}}}^{2}_{curl}$ we have
$\nabla^{*}\cdot{{\boldsymbol{u}}}=0$ and
${{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{v}}}=0$. Finally, we define the
Leray operator
$\mathcal{P}:{{\boldsymbol{L}}}^{2}\to{{\boldsymbol{L}}}^{2}_{div}$,
${{\boldsymbol{u}}}\mapsto{{\boldsymbol{u}}}_{div}$, which maps to the
divergence-free part of a vector field
${{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}$. The operator, which maps to the
curl-free part of a vector field is then simply given by
$(id-\mathcal{P}):{{\boldsymbol{L}}}^{2}\to{{\boldsymbol{L}}}^{2}_{curl}$.

###  2.3. Sobolev spaces on the sphere

Sobolev spaces on the sphere can simply be introduced in terms of (vector)
spherical harmonics. For $\sigma>0$, we define the Sobolev spaces for scalar
valued functions on the sphere by $H^{\sigma}:=\\{f\in
L^{2}\,\big{|}\,\|f\|_{H^{\sigma}}<\infty\\}$, where the norm is induced by
the inner product, which is given by

$\langle
f,g\rangle_{H^{\sigma}}:=\sum_{\ell=0}^{\infty}(1+\lambda_{\ell})^{\sigma}\sum_{k=1}^{N(d,\ell)}\widehat{f}_{\ell,k}\widehat{g}_{\ell,k}.$

In the vector-valued case, we define the divergence-free and the curl-free
Sobolev spaces by

$\begin{split}{{\boldsymbol{H}}}^{\sigma}_{div}&=\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}_{div}\,\big{|}\,\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}<\infty\right\\},\quad{{\boldsymbol{H}}}^{\sigma}_{curl}=\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}_{curl}\,\big{|}\,\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}<\infty\right\\},\end{split}$

where the corresponding inner products are given by

$\begin{split}\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}_{div}}&:=\sum_{\ell=1}^{\infty}(1+\lambda_{\ell})^{\sigma}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{div}_{\ell,k},\quad{{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}_{div},\\\
\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}_{curl}}&:=\sum_{\ell=1}^{\infty}(1+\lambda_{\ell})^{\sigma}\sum_{k=1}^{2\ell+1}\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{curl}_{\ell,k}\quad{{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}_{curl}.\end{split}$

The tangential Sobolev space is then simply given by the direct sum
${{\boldsymbol{H}}}^{\sigma}={{\boldsymbol{H}}}^{\sigma}_{div}\oplus{{\boldsymbol{H}}}^{\sigma}_{curl}$
with inner product
$\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}}=\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}_{div}}+\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}_{curl}}$
for ${{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}$.
Finally, note that if $\sigma\in{\mathbb{N}}_{0}$, the norm for
${{\boldsymbol{H}}}^{\sigma}$ can also be written as

$\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}=\|(1-{\boldsymbol{\Delta}}^{*})^{\sigma}{{\boldsymbol{u}}}\|_{{{\boldsymbol{L}}}^{2}}.$

This also applies to the ${{\boldsymbol{H}}}^{\sigma}_{div}$ and
${{\boldsymbol{H}}}^{\sigma}_{curl}$ norms.

###  2.4. Auxiliary results

For the following analysis, we will need some auxiliary results for fuctions
from the Sobolev space. It is easy to see that if ${{\boldsymbol{u}}}$ is a
tangential vector field, then the zero-th Fourier coefficients of its scalar-
valued components vanishes, i.e. $\widehat{(u_{i})}_{0,0}=0$ for $1\leq i\leq
3$, see [10, Theorem 12.4.1]. The first result, which uses this fact, is a
generalization of [9, Lemma 3.10] and gives us a connection between the
regularity of a vector field and its components. Since the idea of proof does
not change, it will be omitted.

######  Lemma 2.1.

Let $\sigma\geq 0$ and suppose that
${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$. Then,

$\|u_{i}\|_{H^{\sigma}}\leq
C\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}},\quad i=1,2,3.$

The next auxiliary result gives us the regularity of the product of two
Sobolev functions. The proof can be found in [1] and [29].

######  Lemma 2.2.

Let $\sigma>1$ and $f,g\in H^{\sigma}$. Then, $fg\in H^{\sigma}$ and there
exists a constant $C=C(\sigma)>0$ such that

$\|fg\|_{H^{\sigma}}\leq C\|f\|_{H^{\sigma}}\|g\|_{H^{\sigma}}.$

We will need the following extension of Lemma 2.2 for vector fields.

######  Lemma 2.3.

Let $\sigma>1$.

  1. a)

Let ${{\boldsymbol{u}}},{{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}$.
Then, ${{\boldsymbol{u}}}\cdot{{\boldsymbol{v}}}\in H^{\sigma}$ and there
exists a constant $C=C(\sigma)>0$ such that

$\|{{\boldsymbol{u}}}\cdot{{\boldsymbol{v}}}\|_{H^{\sigma}}\leq
C\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.$

  2. b)

Let $f\in H^{\sigma}$ and ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$.
Then, $f{{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$ and there exists a
constant $C=C(\sigma)>0$ such that

$\|f{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\leq
C\|f\|_{H^{\sigma}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.$

###### Proof.

Part a) of this lemma follows easily by the triangle inequality, Lemma 2.1 and
Lemma 2.2, which together yield

$\|{{\boldsymbol{u}}}\cdot{{\boldsymbol{v}}}\|_{H^{\sigma}}\leq\sum_{i=1}^{3}\|u_{i}v_{i}\|_{H^{\sigma}}\leq
C\sum_{i=1}^{3}\|u_{i}\|_{H^{\sigma}}\|v_{i}\|_{H^{\sigma}}\leq
3C\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.$

To prove the second part, we first note that we can bound the norm of a vector
field ${{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}$ by

$\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\leq
C\left(\|\nabla^{*}\cdot{{\boldsymbol{v}}}\|_{H^{\sigma-1}}+\|{{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{v}}}\|_{H^{\sigma-1}}\right).$

Setting ${{\boldsymbol{v}}}=f{{\boldsymbol{u}}}$ and using the product rule
yields

$\begin{split}\|f{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\leq&C\big{(}\|f\nabla^{*}\cdot{{\boldsymbol{u}}}\|_{H^{\sigma-1}}+\|\nabla^{*}f\cdot{{\boldsymbol{u}}}\|_{H^{\sigma-1}}\\\
&+\|f{{\boldsymbol{L}}}^{*}\cdot{{\boldsymbol{u}}}\|_{H^{\sigma-1}}+\|{{\boldsymbol{L}}}^{*}f\cdot{{\boldsymbol{u}}}\|_{H^{\sigma-1}}\big{)}.\end{split}$

Using the first part of this lemma together with Lemma 2.3 finishes the proof.
∎

Finally, to derive high-order estimates for the non-linear operator
$\mathcal{B}$, we will need the following result about the norm of partial
derivatives of a vector field.

######  Lemma 2.4.

Let $\sigma>2$ and ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$. Then,
we have

$\sum_{j=1}^{3}\|\partial_{j}^{*}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}^{2}\leq\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}}^{2}.$

This inequality also holds for the ${{\boldsymbol{H}}}^{\sigma}_{div}$ and
${{\boldsymbol{H}}}^{\sigma}_{curl}$ norms.

###### Proof.

We will show the inequality only for the divergence-free part, the proof for
the curl-free part is nearly the same. With
${{\boldsymbol{v}}}=(1-{\boldsymbol{\Delta}}^{*})^{\sigma}{{\boldsymbol{u}}}$,
we can rewrite the norm as

$\|\partial_{j}^{*}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}^{2}=\langle\partial_{j}^{*}{{\boldsymbol{v}}},\partial_{j}^{*}{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{L}}}^{2}_{div}}=-\langle(\partial_{j}^{*})^{2}{{\boldsymbol{v}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{L}}}^{2}_{div}}$

Summing over $j$ and using that
$\widehat{({\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}})}^{div}_{\ell,k}=-\lambda_{\ell}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}$
then yields

$\begin{split}\sum_{j=1}^{3}\|\partial_{j}^{*}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}^{2}=&-\langle{\boldsymbol{\Delta}}^{*}{{\boldsymbol{v}}},{{\boldsymbol{v}}}\rangle_{{{\boldsymbol{L}}}^{2}_{div}}=-\langle{\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}},{{\boldsymbol{u}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\\
=&-\sum_{\ell=1}^{\infty}(1+\lambda_{\ell})^{\sigma}\sum_{k=1}^{2\ell+1}\widehat{({\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}})}^{div}_{\ell,k}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}\\\
=&\sum_{\ell=1}^{\infty}\lambda_{\ell}(1+\lambda_{\ell})^{\sigma}\sum_{k=1}^{2\ell+1}(\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k})^{2}\leq\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}^{2},\end{split}$

which finishes the proof. ∎

##  3\. Kernel-Based Interpolation and Native Spaces

In this section we will give a brief introduction to the tools needed for the
kernel-based interpolation of vector fields on the sphere. For a general
introduction to kernel-based approximation, see, for example, [33].

###  3.1. Native spaces for scalar valued functions on the sphere

We start with a spherical basis function (SBF)
$\phi:{\mathbb{S}}^{2}\times{\mathbb{S}}^{2}\to{\mathbb{R}}$ of the form

$\phi({{\boldsymbol{x}}},{{\boldsymbol{y}}})=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\widehat{\phi}_{\ell}Y_{\ell,k}({{\boldsymbol{x}}})Y_{\ell,k}({{\boldsymbol{y}}}),\quad{{\boldsymbol{x}}},{{\boldsymbol{y}}}\in{\mathbb{S}}^{2},$

with certain Fourier coefficients $\widehat{\phi}_{\ell}\in{\mathbb{R}}$ such
that $\sum_{\ell=0}^{\infty}(2\ell+1)\widehat{\phi}_{\ell}<\infty$. This
kernel can be seen as the Fourier series of a zonal function, see [33, Ch.
17]. The kernel $\phi$ is positive definite if $\widehat{\phi}_{\ell}>0$ for
all $\ell\geq 0$. In this case the native space of the kernel $\phi$ and its
inner product are given by

$\mathcal{N}_{\phi}:=\left\\{f\in
L^{2}\,\Big{|}\,\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\frac{\widehat{f}_{\ell,k}^{\,2}}{\widehat{\phi}_{\ell}}<\infty\right\\},\quad\langle
f,g\rangle_{\phi}:=\sum_{\ell=0}^{\infty}\sum_{k=1}^{2\ell+1}\frac{\widehat{f}_{\ell,k}\widehat{g}_{\ell,k}}{\widehat{\phi}_{\ell}}.$

This space is a so-called reproducing kernel Hilbert space (RKHS), which on
the one hand means that $\phi(\cdot,{{\boldsymbol{x}}})\in\mathcal{N}_{\phi}$
and on the other hand that $f({{\boldsymbol{x}}})=\langle
f,\phi(\cdot,{{\boldsymbol{x}}})\rangle_{\phi}$ for all
$f\in\mathcal{N}_{\phi}$ and ${{\boldsymbol{x}}}\in{\mathbb{S}}^{2}$. It is a
well known fact that $\mathcal{N}_{\phi}=H^{\sigma}$ if there exists a
constant $C>0$ such that $\widehat{\phi}_{\ell}$ yield a lower and upper bound
of the form

(9) $C^{-1}(1+\lambda_{\ell})^{-\sigma}\leq\widehat{\phi}_{\ell}\leq
C(1+\lambda_{\ell})^{-\sigma},\quad\ell\in{\mathbb{N}}_{0}.$

In this is case, the $H^{\sigma}$-norm and the $\mathcal{N}_{\phi}$-norm are
equivalent.

Spherical basis functions can easily be constructed out of radial basis
functions, since they are zonal functions restricted to the sphere, see [13,
25]. Suitable choices for RBFs are for example the Gaussian or the Wendland
radial basis functions, see [32]. Here, we will use the latter one for the
construction of our numerical scheme. The advantage of these functions is that
they can be constructed of arbitrary smoothness and such that they satisfy (9)
for a given $\sigma>1$, which means that their native space is isomorphic to
the Sobolev space, see [26]. Moreover, they are piecewise polynomials with
compact support, which makes them efficient to calculate. Examples of the
Wendland functions of various smoothness which are suitable for the sphere are
given in Table 1. See also [33] for more details.

Function | $\sigma$
---|---
$\phi_{1}(r)\doteq(1-r)_{+}^{4}(4r+1)$ | $5/2$
$\phi_{2}(r)\doteq(1-r)_{+}^{6}(35r^{2}+18r+3)$ | $7/2$
$\phi_{3}(r)\doteq(1-r)_{+}^{8}(32r^{3}+25r^{2}+8r+1)$ | $9/2$
$\phi_{4}(r)\doteq(1-r)_{+}^{10}(429r^{4}+450r^{3}+210r^{2}+50r+5)$ | $11/2$
Table 1. Compactly supported radial basis functions of various smoothness.

###  3.2. Kernel-based interpolation for tangential vector spaces

Positive definite kernels on the sphere which yield interpolants for
divergence-free vector-valued functions where introduced in [25]. In [13],
positive definite kernels which yield interpolants for rotation free vector
fields were introduced and an error analysis was established. Both, the
divergence-free as well as the rotation free kernel can be constructed from
the scalar-valued kernel $\phi$.

Assuming that
$\sum_{\ell=0}^{\infty}(2\ell+1)\lambda_{\ell}\widehat{\phi}_{\ell}<\infty$,
the kernel which yields divergence-free interpolants is given by the matrix-
valued function

(10)
$\Phi_{div}({{\boldsymbol{x}}},{{\boldsymbol{y}}}):={{\boldsymbol{L}}}^{*}_{{\boldsymbol{x}}}({{\boldsymbol{L}}}^{*}_{{\boldsymbol{y}}})^{T}\phi({{\boldsymbol{x}}},{{\boldsymbol{y}}})=\sum_{\ell=1}^{\infty}\lambda_{\ell}\widehat{\phi}(\ell)\sum_{k=1}^{2\ell+1}{{\boldsymbol{y}}}_{\ell,k}({{\boldsymbol{x}}}){{\boldsymbol{y}}}_{\ell,k}^{T}({{\boldsymbol{y}}}),$

where ${{\boldsymbol{L}}}^{*}_{{\boldsymbol{x}}}$ and
$({{\boldsymbol{L}}}^{*}_{{\boldsymbol{y}}})^{T}$ act on the
${{\boldsymbol{x}}}$ and ${{\boldsymbol{y}}}$ variables, respectively. It is a
well known fact that if $\widehat{\phi}(\ell)>0$ for all $\ell\geq 1$ then
$\Phi_{div}$ is positive definite, which means that given a set of pairwise
distinct points
$\mathcal{X}=\\{{{\boldsymbol{x}}}_{1},\ldots,{{\boldsymbol{x}}}_{N}\\}\subset{\mathbb{S}}^{2}$
and an arbitrary set of tangent vectors ${\boldsymbol{\alpha}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2})\setminus\\{\boldsymbol{0}\\}$, we
have

$\sum_{j,k=1}^{N}{\boldsymbol{\alpha}}_{j}^{T}\Phi_{div}({{\boldsymbol{x}}}_{j},{{\boldsymbol{x}}}_{k}){\boldsymbol{\alpha}}_{k}>0.$

Its curl-free counterpart is given by the matrix-valued kernel

(11)
$\Phi_{curl}({{\boldsymbol{x}}},{{\boldsymbol{y}}}):=\nabla_{{{\boldsymbol{x}}}}^{*}(\nabla_{{{\boldsymbol{y}}}}^{*})^{T}\phi({{\boldsymbol{x}}},{{\boldsymbol{y}}})=\sum_{\ell=1}^{\infty}\lambda_{\ell}\widehat{\phi}(\ell)\sum_{k=1}^{2\ell+1}{{\boldsymbol{z}}}_{\ell,k}({{\boldsymbol{x}}}){{\boldsymbol{z}}}_{\ell,k}^{T}({{\boldsymbol{y}}}).$

Again, if $\widehat{\phi}(\ell)>0$ for all $\ell\geq 1$, then $\Phi_{curl}$ is
positive definite. The native space for both kernels $\Phi_{div}$ and
$\Phi_{curl}$ are given by

$\begin{split}\mathcal{N}_{\Phi_{div}}&=\\!\\!\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}\,\Big{|}\,\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\frac{(\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k})^{2}}{\lambda_{\ell}\widehat{\phi}(\ell)}<\infty\right\\},\\\
\\!\mathcal{N}_{\Phi_{curl}}&=\\!\\!\left\\{{{\boldsymbol{u}}}\in{{\boldsymbol{L}}}^{2}\,\Big{|}\,\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\frac{(\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k})^{2}}{\lambda_{\ell}\widehat{\phi}(\ell)}<\infty\right\\},\end{split}$

where the norms are induced by the inner products, which are given by

$\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{\Phi_{div}}\\!=\\!\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\frac{1}{\lambda_{\ell}\widehat{\phi}(\ell)}\widehat{{{\boldsymbol{u}}}}^{div}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{div}_{\ell,k},\quad\\!\langle{{\boldsymbol{u}}},{{\boldsymbol{v}}}\rangle_{\Phi_{curl}}\\!=\\!\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\frac{1}{\lambda_{\ell}\widehat{\phi}(\ell)}\widehat{{{\boldsymbol{u}}}}^{curl}_{\ell,k}\widehat{{{\boldsymbol{v}}}}^{curl}_{\ell,k}.$

Both spaces are again RKHS with reproducing kernels $\Phi_{div}$ and
$\Phi_{curl}$, respectively, which means in the case of $\Phi_{div}$ that for
every ${{\boldsymbol{u}}}\in\mathcal{N}_{\Phi_{div}}$,
${{\boldsymbol{x}}}\in{\mathbb{S}}^{2}$ and ${\boldsymbol{\alpha}}\in
T_{{\boldsymbol{x}}}({\mathbb{S}}^{2})$ we have
$\Phi_{div}(\cdot,{{\boldsymbol{x}}}){\boldsymbol{\alpha}}\in\mathcal{N}_{\Phi_{div}}$
on the one hand and
${\boldsymbol{\alpha}}^{T}{{\boldsymbol{u}}}({{\boldsymbol{x}}})=\langle{{\boldsymbol{u}}},\Phi_{div}(\cdot,{{\boldsymbol{x}}}){\boldsymbol{\alpha}}\rangle_{\Phi_{div}}$
on the other hand [15]. If the Fourier coefficients behave like
$\widehat{\phi}_{\ell}\sim(1+\lambda_{\ell})^{-(\sigma+1)}$ for a $\sigma>0$,
then $\mathcal{N}_{\Phi_{div}}={{\boldsymbol{H}}}^{\sigma}_{div}$ and
$\mathcal{N}_{\Phi_{curl}}={{\boldsymbol{H}}}^{\sigma}_{curl}$ with equivalent
norms, which means there is a constant $C>0$ with

$\begin{split}C^{-1}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}&\leq\|{{\boldsymbol{u}}}\|_{\Phi_{div}}\leq
C\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}},\quad{{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}_{div},\\\
C^{-1}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}&\leq\|{{\boldsymbol{u}}}\|_{\Phi_{curl}}\leq
C\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}},\quad{{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}_{curl}.\end{split}$

Finally, the kernel $\Phi=\Phi_{div}+\Phi_{curl}$ is also positive definite if
$\Phi_{div}$ and $\Phi_{curl}$ are positive definite. The native space
$\mathcal{N}_{\Phi}$ of $\Phi$ is then given by the direct sum of the native
spaces $\mathcal{N}_{\Phi_{div}}$ and $\mathcal{N}_{\Phi_{curl}}$, i.e.
$\mathcal{N}_{\Phi}=\mathcal{N}_{\Phi_{div}}\oplus\mathcal{N}_{\Phi_{curl}}$.
Hence, $\mathcal{N}_{\Phi}$ can be identified with
${{\boldsymbol{H}}}^{\sigma}$ and their norms are equivalent. The exact
representation for $\Phi_{div}$ and $\Phi_{curl}$ defined by a zonal kernel
$\phi$ can be found in [13].

###  3.3. Interpolation with Kernel Functions

From now on, we assume that we have a positive definite matrix-valued kernel
$\Phi=\Phi_{div}+\Phi_{curl}$ with
$\widehat{\phi}_{\ell}\sim(1+\lambda_{\ell})^{-(\sigma+1)}$ for a $\sigma>0$.
In the following, we define a discrete set
$\mathcal{X}=\\{{{\boldsymbol{x}}}_{1},\ldots,{{\boldsymbol{x}}}_{N}\\}\subset{\mathbb{S}}^{2}$
of $N\in{\mathbb{N}}$ pairwise distinct data points with fill distance

$h_{\mathcal{X}}:=\sup_{{{\boldsymbol{x}}}\in{\mathbb{S}}^{2}}\min_{{{\boldsymbol{x}}}_{j}\in\mathcal{X}}d({{\boldsymbol{x}}},{{\boldsymbol{x}}}_{j}),$

where $d({{\boldsymbol{x}}},{{\boldsymbol{y}}})$ is the geodesic distance
between ${{\boldsymbol{x}}},{{\boldsymbol{y}}}\in{\mathbb{S}}^{2}$.
Furthermore, we define our approximation space

${\boldsymbol{\mathcal{F}}}_{\Phi}(\mathcal{X}):=\operatorname{span}\\{\Phi(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}\>|\>{{\boldsymbol{x}}}_{j}\in\mathcal{X},\,{\boldsymbol{\alpha}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2}),\,1\leq j\leq N\\}.$

The interpolation operator, which maps a function
${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$ to its best approximation
in ${\boldsymbol{\mathcal{F}}}_{\Phi}(\mathcal{X})$, is given by

$\mathcal{I}_{\Phi}:{{\boldsymbol{H}}}^{\sigma}\to{\boldsymbol{\mathcal{F}}}_{\Phi}(\mathcal{X}),\quad{{\boldsymbol{u}}}\mapsto{{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}=\sum_{j=1}^{N}\Phi(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j},$

where the coefficients ${\boldsymbol{\alpha}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2})$ are uniquely determined by
${{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}({{\boldsymbol{x}}}_{j})={{\boldsymbol{u}}}({{\boldsymbol{x}}}_{j})$
for all $1\leq j\leq N$.

Since the kernel $\Phi$ is just the sum of its divergence-free and its curl-
free part, we also have a simple approximation for the divergence-free and the
curl-free part of ${{\boldsymbol{u}}}$. Note that the operator
$\mathcal{I}_{\Phi}$ as well as the following operators depend on the data
sites $\mathcal{X}$. However, since we do not vary the data sites, we will
omit this dependency due to readability.

Projecting an interpolant onto the space of divergence-free functions, we can
use that
$\mathcal{P}\Phi({{\boldsymbol{x}}}_{j},\cdot){\boldsymbol{\alpha}}_{j}=\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}$.
Hence, we define the interpolation operator of the divergence-free part, which
can also be seen as an approximated Leray operator, by

(12)
$\mathcal{P}_{\Phi}:{{\boldsymbol{H}}}^{\sigma}\to{\boldsymbol{\mathcal{F}}}_{\Phi_{div}}(\mathcal{X}),\quad{{\boldsymbol{u}}}\mapsto\sum_{j=1}^{N}\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}.$

In the same way, we can project an interpolant on its curl-free part by using
the operator

$(\operatorname{id}-\mathcal{P})\mathcal{I}_{\Phi}:{{\boldsymbol{H}}}^{\sigma}\to{\boldsymbol{\mathcal{F}}}_{\Phi_{curl}}(\mathcal{X}),\quad{{\boldsymbol{u}}}\mapsto\sum_{j=1}^{N}\Phi_{curl}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}.$

From [13], we will use the following error estimates for the three
interpolation operators, which are given in terms of the fill distance
$h_{\mathcal{X}}$.

######  Theorem 3.1.

Let $\sigma>3/2$ and ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$. Then,
there exists a $C>0$ such that for all $\tau<\sigma$

$\|{{\boldsymbol{u}}}-\mathcal{I}_{\Phi}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau}}\leq
Ch_{\mathcal{X}}^{\sigma-\tau}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.$

In particular, there also exist constants $C_{1},C_{2}>0$ such that

$\begin{split}\|\mathcal{P}{{\boldsymbol{u}}}-\mathcal{P}_{\Phi}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau}}&\leq
C_{1}h_{\mathcal{X}}^{\sigma-\tau}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\\\
\|({\operatorname{id}}-\mathcal{P}){{\boldsymbol{u}}}-({\operatorname{id}}-\mathcal{P}_{\Phi}){{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau}}&\leq
C_{2}h_{\mathcal{X}}^{\sigma-\tau}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\end{split}$

These estimates also hold in the ${{\boldsymbol{H}}}^{\sigma}_{div}$ and
${{\boldsymbol{H}}}^{\sigma}_{curl}$ norms.

###  3.4. Collocation methods for Helmholtz Equation

For the convergence result, we will need an approximation method for the
Helmholtz equation. For a given right hand side
${{\boldsymbol{f}}}:{\mathbb{S}}^{2}\to{\mathbb{R}}^{3}$ the vectorial
Helmholtz equation is given by

(13)
$-{\boldsymbol{\Delta}}^{*}{{\boldsymbol{u}}}+{{\boldsymbol{u}}}={{\boldsymbol{f}}},\quad\text{on
}{\mathbb{S}}^{2}.$

Its solution can easily be obtained via Fourier analysis. Every right hand
side ${{\boldsymbol{f}}}$ with zero divergence admits a divergence-free
solution ${{\boldsymbol{u}}}$.

Let $\widehat{\phi}_{\ell}\sim(1+\lambda_{\ell})^{-(\sigma+1)}$ with
$\sigma>5/2$ and
${{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}=\sum_{j=1}^{N}\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}\in{\boldsymbol{\mathcal{F}}}_{\Phi_{div}}$.
We want to find the coefficient-vectors ${\boldsymbol{\alpha}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2})$, $1\leq j\leq N$, such that
${{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}$ is an approximation in
${\boldsymbol{\mathcal{F}}}_{\Phi_{div}}$ to
${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma+1}_{div}\subset{{\boldsymbol{H}}}^{\sigma}_{div}$
with the property that
$(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}){{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}({{\boldsymbol{x}}}_{j})=(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}){{\boldsymbol{u}}}({{\boldsymbol{x}}}_{j})$
for every ${{\boldsymbol{x}}}_{j}\in\mathcal{X}$.

Using (10) and the fact that
${\boldsymbol{\Delta}}^{*}{{\boldsymbol{y}}}_{\ell,m}=-\lambda_{\ell}{{\boldsymbol{y}}}_{\ell,m}$,
we easily see that

$\left(-{\boldsymbol{\Delta}}^{*}_{{\boldsymbol{x}}}+\operatorname{id}_{{\boldsymbol{x}}}\right)\Phi_{div}({{\boldsymbol{x}}},{{\boldsymbol{y}}})=\sum_{\ell=1}^{\infty}\lambda_{\ell}(1+\lambda_{\ell})\widehat{\phi}(\ell)\sum_{k=1}^{2\ell+1}{{\boldsymbol{y}}}_{\ell,m}({{\boldsymbol{x}}}){{\boldsymbol{y}}}_{\ell,m}^{T}({{\boldsymbol{y}}}).$

If $\Phi_{div}$ is positive definite, so is
$\left(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}\right)\Phi_{div}$, and we
can find uniquely determined ${\boldsymbol{\alpha}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2})$, $1\leq j\leq N$, with

(14)
$\left(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}\right){{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}({{\boldsymbol{x}}}_{j})={{\boldsymbol{f}}}({{\boldsymbol{x}}}_{j})=\left(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}\right){{\boldsymbol{u}}}({{\boldsymbol{x}}}_{j}),\quad{{\boldsymbol{x}}}_{j}\in\mathcal{X}.$

The operator

$\mathcal{R}_{\Phi_{div}}:{{\boldsymbol{H}}}^{\sigma+1}_{div}\to{\boldsymbol{\mathcal{F}}}_{\Phi_{div}},\quad{{\boldsymbol{u}}}\to{{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}},$

which gives the approximation
${{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}$ to the solution
${{\boldsymbol{u}}}$ of the Helmholtz-equation, is called Ritz-projection
operator.

######  Theorem 3.2.

Let $\sigma>5/2$. If
${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma+1}_{div}$ then there exists
exactly one approximation
${{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}=\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{F}}}_{\Phi_{div}}$
which satisfies
$(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}){{\boldsymbol{s}}}_{\mathcal{X},{{\boldsymbol{u}}}}({{\boldsymbol{x}}}_{j})=(-{\boldsymbol{\Delta}}^{*}+\operatorname{id}){{\boldsymbol{u}}}({{\boldsymbol{x}}}_{j})$
for every ${{\boldsymbol{x}}}_{j}\in\mathcal{X}$. Moreover, there exists a
constant $C>0$ such that

$\|{{\boldsymbol{u}}}-\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau_{1}}}\leq
Ch_{\mathcal{X}}^{\tau_{2}-\tau_{1}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau}_{2}},\quad{{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\tau_{2}}_{div}$

for all $5/2\leq\tau_{1}\leq\sigma+1\leq\tau_{2}\leq 2\sigma$.

###### Proof.

As we can easily verify, if
${{\boldsymbol{f}}}\in{{\boldsymbol{H}}}^{\sigma}_{div}$ and
${{\boldsymbol{u}}}$ is a solution of (13), then
${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma+2}_{div}$ and we have

$\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma+2}_{div}}=\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}.$

Hence , we can apply Theorem 3.1, which yields

$\begin{split}\|{{\boldsymbol{u}}}-\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau_{1}}}&=\|{{\boldsymbol{f}}}-\mathcal{I}_{(-{\boldsymbol{\Delta}}^{*}+\operatorname{id})\Phi_{div},\mathcal{X}}{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\tau_{1}-2}}\\\
&\leq
Ch^{(\tau_{2}-2)-(\tau_{1}-2)}\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\tau_{2}-2}}\\\
&=Ch^{\tau_{2}-\tau_{1}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\tau_{2}}}.\end{split}$

The uniquiness of the solution comes from the fact that the kernel
$\left(-{\boldsymbol{\Delta}}^{*}_{{\boldsymbol{x}}}+\operatorname{id}_{{\boldsymbol{x}}}\right)\Phi_{div}$
is positive definite. ∎

##  4\. The Navier-Stokes equations on the Sphere and auxiliary Results

As mentioned in the introduction, we want to solve the incompressible Navier-
Stokes equations on the rotating sphere (1) - (3). We do this by applying the
Leary operator $\mathcal{P}$ to equation (1). By (2), the velocity itself is
divergence-free, so the time derivative of the velocity and the diffusion term
are not affected, but the gradient of the pressure vanishes. We arrive at

(15)
$\displaystyle\partial_{t}{{\boldsymbol{u}}}-\nu\Delta_{*}{{\boldsymbol{u}}}$
$\displaystyle=\mathcal{P}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{C}({{\boldsymbol{u}}})\right)$
$\displaystyle\text{in }{\mathbb{S}}^{2}\times(0,T),$ (16)
$\displaystyle{{\boldsymbol{u}}}(\cdot,0)$
$\displaystyle={{\boldsymbol{u}}}_{0}$ $\displaystyle\text{on
}{\mathbb{S}}^{2}.$

Hence, the velocity can be calculated by finding a solution of (15) - (16). On
the other hand, applying the operator $\mathcal{I}-\mathcal{P}$ to (1) yields
the equation

(17)
$\nabla_{*}p=(id-\mathcal{P})\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{C}({{\boldsymbol{u}}})\right)\quad\text{in
}{\mathbb{S}}^{2}\times(0,T),$

from which the pressure can be calculated. Since the pressure can only be
calculated up to a constant, we will assume the pressure to have mean zero
integral, i.e.
$\int_{{\mathbb{S}}^{2}}p({{\boldsymbol{x}}})dS({{\boldsymbol{x}}})=0$. The
pair $({{\boldsymbol{u}}},p)$ which will be given by the solutions of (15) and
(17) is the a solution of the Navier-Stokes equations (1) - (3).

The existence and uniqueness of the solution of the divergence-free part of
the Navier-Stokes equations has been discussed in [19, 17, 18]. However, the
results are often based on the weak formulation and do not yield smooth
solutions. Since we need high-order regularity of our solution, we will use
the result of [5], which is based on the Gevrey regularity. The following
lemma is a simple deduction of their existence and uniqueness result.

######  Lemma 4.1.

Let $\sigma>0$. Assume that
${{\boldsymbol{u}}}_{0}\in{{\boldsymbol{H}}}^{2\sigma+1}_{div}$ and
${{\boldsymbol{f}}}\in L^{\infty}(0,\infty;{{\boldsymbol{L}}}^{2}_{div})$ such
that

$\sum_{\ell=1}^{\infty}\sum_{k=1}^{2\ell+1}\lambda_{\ell}^{2\sigma}e^{\tau\lambda_{\ell}^{1/2}}(\widehat{{{\boldsymbol{f}}}(\cdot,t)}_{\ell,k}^{div})^{2}<\infty$

for some $\tau>0$ and all $t\in(0,\infty)$. Then, (15) - (16) has a unique
global solution $u\in
L^{\infty}(0,\infty;{{\boldsymbol{H}}}^{2\sigma+1}_{div})$.

Hence, this result ensures that there exists a solution for the divergence-
free part of the Navier-Stokes equations (15) - (16). The existence of the
solution of (17) is simply given by the fact that the right hand side of (17)
is curl-free.

Next we want to gain further information about the two operators $\mathcal{B}$
and $\mathcal{C}$. As one can easily see,
$\mathcal{B}:{{\boldsymbol{H}}}^{\sigma}\times{{\boldsymbol{H}}}^{\sigma+1}\to{{\boldsymbol{H}}}^{\sigma}$
is a bilinear operator. In the following lemma we will show that the operator
is also bounded.

######  Lemma 4.2.

Let $\sigma>1$. Then, there exist a constant $C_{1}>0$, which only depends on
$d$ and $\sigma$, such that

$\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}})\|_{{{\boldsymbol{H}}}^{\sigma}}\leq
C_{1}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\left(\sum_{j=1}^{3}\|\partial_{j}^{*}{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}^{2}\right)^{1/2}\leq
C_{1}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}},\quad,$

for all ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$ and
${{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma+1}$. Moreover, there exist a
constant $C_{2}>0$, which depends on $d$ and $\sigma$, such that

$\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}}),{{\boldsymbol{w}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}}\leq
C_{2}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma-1}}\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\|{{\boldsymbol{w}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}}.$

for all ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma-1}$,
${{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}$ and
${{\boldsymbol{w}}}\in{{\boldsymbol{H}}}^{\sigma+1}$.

###### Proof.

For the first term, the triangle inequality, Lemma 2.2 and the Cauchy-Schwarz
inequality yield

$\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}})\|_{{{\boldsymbol{H}}}^{\sigma}}\leq\sum_{i=1}^{3}\|u_{i}\partial^{*}_{i}{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\leq
C\left(\sum_{i=1}^{3}\|u_{i}\|^{2}_{H^{\sigma}}\right)^{1/2}\left(\sum_{i=1}^{3}\|\partial^{*}_{i}{{\boldsymbol{v}}}\|^{2}_{{{\boldsymbol{H}}}^{\sigma}}\right)^{1/2}.$

Using Lemma 2.1 gives the first inequality of the first statement. The second
one follows from Lemma 2.4. For the second estimate, we use that

$\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}}),{{\boldsymbol{w}}}\rangle_{{{\boldsymbol{H}}}^{\sigma}}\leq\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{v}}})\|_{{{\boldsymbol{H}}}^{\sigma-1}}\|{{\boldsymbol{w}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}}.$

Applying the first estimate finishes the proof. ∎

It is also easy to see that
$\mathcal{C}:{{\boldsymbol{H}}}^{\sigma}\to{{\boldsymbol{H}}}^{\sigma}$ is a
bounded linear operator. This will be given in the following lemma.

######  Lemma 4.3.

Let $\sigma>1$ and ${{\boldsymbol{u}}}\in{{\boldsymbol{H}}}^{\sigma}$. Then,
there exists a constant $C>0$, which only depends on $d$ and $\sigma$ such
that

$\begin{split}\|\mathcal{C}({{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}}&\leq
C\Omega\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.\end{split}$

###### Proof.

Using Lemma 2.3, we have

$\|\mathcal{C}({{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}}=2\Omega\|n_{3}{{\boldsymbol{n}}}\times{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}\leq
2\Omega\|n_{3}\|_{H^{\sigma}}\|{{\boldsymbol{n}}}\times{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}.$

Moreover, using that
${{\boldsymbol{n}}}\times{{\boldsymbol{y}}}_{\ell,k}={{\boldsymbol{z}}}_{\ell,k}$
for $\ell\in{\mathbb{N}}$ and $1\leq k\leq 2\ell+1$, the divergence-free
Fourier coefficients of ${{\boldsymbol{n}}}\times{{\boldsymbol{u}}}$ satisfy

$\begin{split}(\widehat{{{\boldsymbol{n}}}\times{{\boldsymbol{u}}}})_{\ell,k}^{div}=&\int_{{\mathbb{S}}^{2}}\left({{\boldsymbol{n}}}\times{{\boldsymbol{u}}}\right)({{\boldsymbol{n}}})\cdot{{\boldsymbol{y}}}_{\ell,k}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}})\\\
=&\int_{{\mathbb{S}}^{2}}{{\boldsymbol{u}}}({{\boldsymbol{x}}})\cdot\left({{\boldsymbol{n}}}\times{{\boldsymbol{y}}}_{\ell,k}\right)({{\boldsymbol{x}}})dS({{\boldsymbol{x}}})=\widehat{{{\boldsymbol{u}}}}_{\ell,k}^{curl}.\end{split}$

With the same argument, one can also show that
$(\widehat{{{\boldsymbol{n}}}\times{{\boldsymbol{u}}}})_{\ell,k}^{curl}=\widehat{{{\boldsymbol{u}}}}_{\ell,k}^{div}$
which means that
$\|{{\boldsymbol{n}}}\times{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}=\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}}$.
∎

##  5\. The Semi-discrete Problem

We now want to give the semi-discretised collocation scheme for the
incompressible Navier-Stokes equations on the sphere. We will differ between
the divergence-free part to get an approximation of the velocity and the curl-
free part to get an approximation of the pressure.

###  5.1. Approximation Scheme for the Velocity

As mentioned in the introduction, the idea is to approximate the velocity with
the interpolation scheme given in section 3.3 and using the method of lines by

(18)
${{\boldsymbol{u}}}_{h}(\cdot,t)=\sum_{j=1}^{N}\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\alpha}}_{j}(t),\quad
t\in[0,T_{h}),$

where $T_{h}>0$ is the maximal time of existence of ${{\boldsymbol{u}}}_{h}$.
By its construction, the approximation is divergence-free and automatically
satisfies equation (2). Moreover, to get an approximation of the right hand
side in (15), we exchange the Leray operator by its approximated version from
(12). Then, we arrive at the following problem: We search for the solution
${{\boldsymbol{u}}}_{h}\in C^{1}([0,\infty);\mathcal{F}_{\Phi_{div}}(X))$
which satisfies the discretised equations

(19)
$\displaystyle\partial_{t}{{\boldsymbol{u}}}_{h}-\nu\Delta^{*}{{\boldsymbol{u}}}_{h}$
$\displaystyle=\mathcal{P}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right),$
$\displaystyle\text{in }\mathcal{X}\times(0,T_{h}),$ (20)
$\displaystyle{{\boldsymbol{u}}}_{h}(\cdot,0)$
$\displaystyle={{\boldsymbol{u}}}_{0},$ $\displaystyle\text{in }\mathcal{X}.$

To be more precise, we now have a system of ODEs in time, where we have to
find the coefficient vectors ${\boldsymbol{\alpha}}_{j}\in
C^{1}([0,\infty);T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2}))$,
$j=1,\ldots,N$. Standard theory of ordinary differential equation ensures the
local-in-time existence of a solution up to a maximum time which we will
denote by $T_{h}$. This time depends on the initial data, the right hand side
as well as the data sites $\mathcal{X}$. However, analogue to the continuation
criterion for ODEs, one can show that the solution of (19) - (20) has to exist
globally in time if the solution is bounded on the maximum time-interval
$[0,T_{h})$.

We will investigate the difference between the solution ${{\boldsymbol{u}}}$
of the projected Navier-Stokes equations (15) - (16) and our approximated
solution ${{\boldsymbol{u}}}_{h}$ of (19) - (20)

$\|{{\boldsymbol{u}}}(\cdot,t)-{{\boldsymbol{u}}}_{h}(\cdot,t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq\|{\boldsymbol{\rho}}_{h}(\cdot,t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+\|{\boldsymbol{\theta}}_{h}(\cdot,t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}},$

where we split up the error into the projection error
${\boldsymbol{\rho}}_{h}={{\boldsymbol{u}}}-\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}$
and the stability error
${\boldsymbol{\theta}}_{h}=\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}$.
The projection error can easily be estimated by Theorem 3.2, so it remains to
investigate the stability error. The procedure will be as follows. To find an
estimate for the stability error ${\boldsymbol{\theta}}_{h}$, we will
investigate

(21)
$\begin{split}\frac{1}{2}\partial_{t}\|{\boldsymbol{\theta}}_{h}\|^{2}_{\Phi_{div}}-\langle\nu\boldsymbol{\Delta}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\left\langle-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}+{{\boldsymbol{f}}}-\mathcal{I}_{\Phi_{div}}{{\boldsymbol{f}}},{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}}\\\
&-\left\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}}\\\
&-\left\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}}.\end{split}$

The identity (21) will be shown in Lemma 5.1. Estimates for the last two terms
of the right hand side will be given separately in Lemma 5.2 and Lemma 5.3.
After that, we will give a local-in-time error estimate for the approximated
solution in Lemma 5.4, before we will use a bootstrap argument to show that
the estimate is valid on the whole time interval where the analytical solution
exists. Theorem 5.5 will be the main convergence result of this paper. Note
that, according to Lemma 4.1, the analytical solution will exist global in
time under certain conditions. In this case we will show that the approximated
solution will exist global in time as well.

Note that in (21) the term
$-\nu\langle\boldsymbol{\Delta}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}$
remains positive since

$\begin{split}-\nu\langle\boldsymbol{\Delta}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}&=-\sum_{\ell=1}^{\infty}\frac{1}{\lambda_{\ell}\widehat{\phi}_{\ell}}\sum_{k=1}^{2\ell+1}({\boldsymbol{\theta}}_{h})_{\ell,k}^{div}\int_{{\mathbb{S}}^{2}}{\boldsymbol{\Delta}}^{*}{\boldsymbol{\theta}}_{h}({{\boldsymbol{x}}})\cdot{{\boldsymbol{y}}}_{\ell,k}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}})\\\
&=\sum_{\ell=1}^{\infty}\frac{1}{\widehat{\phi}_{\ell}}\sum_{k=1}^{2\ell+1}({\boldsymbol{\theta}}_{h})_{\ell,k}^{div}\int_{{\mathbb{S}}^{2}}{\boldsymbol{\theta}}_{h}({{\boldsymbol{x}}})\cdot\frac{-{\boldsymbol{\Delta}}^{*}{{\boldsymbol{y}}}_{\ell,k}({{\boldsymbol{x}}})}{\lambda_{\ell}}dS({{\boldsymbol{x}}})\\\
&=\sum_{\ell=1}^{\infty}\frac{1}{\widehat{\phi}_{\ell}}\sum_{k=1}^{2\ell+1}\left(\int_{{\mathbb{S}}^{2}}{\boldsymbol{\theta}}_{h}({{\boldsymbol{x}}})\cdot{{\boldsymbol{y}}}_{\ell,k}({{\boldsymbol{x}}})dS({{\boldsymbol{x}}})\right)^{2}>0.\end{split}$

Unfortunately, the term on the right hand side is none of our native space
norms of the kernel $\phi$. However, it can define a norm which is equivalent
to the ${{\boldsymbol{H}}}^{\sigma+1}_{div}$ norm. Hence, we can find a
constant $C>0$ such that

(22)
$C^{-1}\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}^{2}\leq-\langle\boldsymbol{\Delta}^{*}{{\boldsymbol{v}}},{{\boldsymbol{v}}}\rangle_{\Phi_{div}}\leq
C\|{{\boldsymbol{v}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}^{2}$

for all ${{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma+1}_{div}$.

######  Lemma 5.1.

Let $\sigma>1$, ${{\boldsymbol{f}}}\in
L^{\infty}(0,T_{h};{{\boldsymbol{H}}}^{\sigma})$ and ${{\boldsymbol{u}}}\in
C^{1}([0,T_{h}),{{\boldsymbol{H}}}^{\sigma+1})$. For every test function
$\chi\in\mathcal{F}_{\Phi_{div}}(\mathcal{X})$, the equation

$\begin{split}\langle\partial_{t}{\boldsymbol{\theta}}_{h}-\nu\boldsymbol{\Delta}^{*}{\boldsymbol{\theta}}_{h},\chi\rangle_{\Phi_{div}}=&\langle-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}+\mathcal{P}{{\boldsymbol{f}}}-\mathcal{P}_{\Phi}{{\boldsymbol{f}}},\chi\rangle_{\Phi_{div}}\\\
&-\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}),\chi\rangle_{\Phi_{div}}\\\
&-\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}_{h}),\chi\rangle_{\Phi_{div}}\end{split}$

holds almost everywhere on $[0,T_{h})$.

###### Proof.

Using the property (14) of the discretised Ritz projection, we note that

$-\boldsymbol{\Delta}^{*}\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}=-\boldsymbol{\Delta}^{*}{{\boldsymbol{u}}}+{{\boldsymbol{u}}}-\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}=-\boldsymbol{\Delta}^{*}{{\boldsymbol{u}}}+{\boldsymbol{\rho}}_{h}$

on $\mathcal{X}$. Using (19), we can rewrite the left hand to

$\begin{split}\partial_{t}{\boldsymbol{\theta}}_{h}-\nu\boldsymbol{\Delta}^{*}{\boldsymbol{\theta}}_{h}&=\\!\partial_{t}\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}-\nu\boldsymbol{\Delta}^{*}\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}-\mathcal{P}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right)\\\
&=\\!\partial_{t}\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}-\nu\boldsymbol{\Delta}^{*}{{\boldsymbol{u}}}+\nu{\boldsymbol{\rho}}_{h}-\mathcal{P}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right).\end{split}$

Solving (15) for $-\nu\boldsymbol{\Delta}^{*}{{\boldsymbol{u}}}$, and
inserting it into the equality above yields

$\begin{split}\partial_{t}{\boldsymbol{\theta}}_{h}-\nu\Delta^{*}{\boldsymbol{\theta}}_{h}=&-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}+\mathcal{P}{{\boldsymbol{f}}}-\mathcal{P}_{\Phi}{{\boldsymbol{f}}}\\\
&-\mathcal{P}(\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}}))+\mathcal{P}_{\Phi}(\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}))-\mathcal{P}(\mathcal{C}({{\boldsymbol{u}}}))+\mathcal{P}_{\Phi}(\mathcal{C}({{\boldsymbol{u}}}_{h}))\end{split}$

pointwise on $\mathcal{X}$. Hence, for every
$\chi\in\mathcal{F}_{\Phi_{div}}(\mathcal{X})$ we have

$\begin{split}\langle\partial_{t}{\boldsymbol{\theta}}_{h}-\nu\boldsymbol{\Delta}^{*}{{\boldsymbol{u}}},\chi\rangle_{\Phi_{div}}=&\langle-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}+\mathcal{P}{{\boldsymbol{f}}}-\mathcal{P}_{\Phi}{{\boldsymbol{f}}},\chi\rangle_{\Phi_{div}}\\\
&-\langle\mathcal{P}(\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}}))-\mathcal{P}_{\Phi}(\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})),\chi\rangle_{\Phi_{div}}\\\
&-\langle\mathcal{P}(\mathcal{C}({{\boldsymbol{u}}}))-\mathcal{P}_{\Phi}(\mathcal{C}({{\boldsymbol{u}}}_{h})),\chi\rangle_{\Phi_{div}}.\end{split}$

The fact that we have
$\langle\mathcal{P}{{\boldsymbol{g}}},\chi\rangle_{\Phi_{div}}=\langle{{\boldsymbol{g}}},\chi\rangle_{\Phi_{div}}$
for all ${{\boldsymbol{g}}}\in{{\boldsymbol{H}}}^{\sigma}$ and
$\chi\in\mathcal{F}_{\Phi_{div}}(\mathcal{X})$ together with
$\mathcal{P}_{\Phi}=\mathcal{P}\mathcal{I}_{\Phi}=\mathcal{I}_{\Phi_{div}}$
completes the proof. ∎

As we have now proved the equality in (21), we can go on in finding estimates
for the terms on the right hand side. We will start with the second one, which
addresses the convective part. Here, we will start with finding an estimate
under the condition that the approximated solution is sufficiently good.

######  Lemma 5.2.

Let $\sigma>1$, $t_{0}\in(0,T)$ and
${{\boldsymbol{u}}}(\cdot,t_{0})\in{{\boldsymbol{H}}}^{2\sigma+1}$. Let $M\geq
1$ be a positive number such that
$\|{{\boldsymbol{u}}}(\cdot,t_{0})-{{\boldsymbol{u}}}_{h}(\cdot,t_{0})\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
M$. Then, for $t=t_{0}$, there exists a $C=C({{\boldsymbol{u}}})>0$,
independent of $\mathcal{X}$ and $M$ but linear in
$\|{{\boldsymbol{u}}}\|_{C([0,T];H^{2\sigma}_{div})}$, such that

$\begin{split}\big{|}\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}&\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\big{|}\\\
\leq&Ch_{\mathcal{X}}^{\sigma}M\left(1+\|{{\boldsymbol{u}}}\|_{H^{2\sigma+1}_{div}}\right)\left(\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}+\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\right)\\\
&+C\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}+CM\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}.\end{split}$

###### Proof.

Firstly, we will split up the left hand side into

$\begin{split}\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
&+\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
&+\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}.\end{split}$

We will estimate the three terms separately. With Theorem 3.1 and Lemma 4.2,
the first term can be estimated by

$\begin{split}|\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}|&\leq\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})\|_{\Phi_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
Ch^{\sigma}_{\mathcal{X}}\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})\|_{H^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
Ch^{\sigma}_{\mathcal{X}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma+1}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}.\end{split}$

Using the boundedness of the interpolation operator
$\|\mathcal{I}_{\Phi_{div}}{{\boldsymbol{v}}}\|_{\Phi_{div}}\leq\|{{\boldsymbol{v}}}\|_{\Phi_{div}}$
for all ${{\boldsymbol{v}}}\in{{\boldsymbol{H}}}^{\sigma}_{div}$, we have for
the second term

$\begin{split}|\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}|&\leq\|\mathcal{B}({{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}})\|_{\Phi_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
C\|{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}},\end{split}$

where, using the decomposition
${{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}={\boldsymbol{\rho}}_{h}+{\boldsymbol{\theta}}_{h}$,
the first norm can be estimated with Theorem 3.1 by

(23)
$\|{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq\|{\boldsymbol{\rho}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
Ch^{\sigma}_{\mathcal{X}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}+\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}.$

For the third term we need a little more effort. We begin by splitting it up
by

(24)
$\begin{split}\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
&+\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\theta}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}.\end{split}$

Using the self-adjointness of the interpolation operator
$\mathcal{I}_{\Phi_{div}}$ to split up the first part yields

$\begin{split}\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\langle\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),\mathcal{I}_{\Phi_{div}}{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
=&\langle\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
&+\langle\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),\mathcal{I}_{\Phi_{div}}{\boldsymbol{\theta}}_{h}-{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}},\end{split}$

where we have on the one hand the estimate

$\begin{split}|\langle\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}|&\leq
C\|{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma-1}_{div}}\|{\boldsymbol{\rho}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\\\
&\leq
Ch_{\mathcal{X}}^{\sigma}(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma-1}_{div}}+M)\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\end{split}$

and on the other hand

$\begin{split}|\langle\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\rho}}_{h}),\mathcal{I}_{\Phi_{div}}{\boldsymbol{\theta}}_{h}-{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}|&\leq
C\|{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{\boldsymbol{\rho}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\|\mathcal{I}_{\Phi_{div}}{\boldsymbol{\theta}}_{h}\\!-\\!{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\\
&\leq
C(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\!\\!+\\!M)h_{\mathcal{X}}^{\sigma-1}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}h_{\mathcal{X}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\\\
&\leq
Ch_{\mathcal{X}}^{\sigma}(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+M)\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}.\end{split}$

Finally, using the Cauchy-Schwarz inequality once again, the second part of
(24) can be bounded by

$\begin{split}|\langle\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{\boldsymbol{\theta}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}|&\leq
C\|{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
C(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+M)\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}.\end{split}$

Summing up all the inequalities from above finishes the proof. ∎

Now we will give an estimate for the third term of the left hand side of (21).

######  Lemma 5.3.

Let $\sigma>1$, $t\in(0,T)$ and
${{\boldsymbol{u}}}(\cdot,t)\in{{\boldsymbol{H}}}^{\sigma+1}$. There exists a
constant $C>0$ independent of ${{\boldsymbol{u}}}$, $\mathcal{X}$ and $t$,
such that

$\big{|}\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\big{|}\leq
Ch_{\mathcal{X}}^{\sigma}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}+C\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}.$

###### Proof.

First, we split up the left hand side into

$\begin{split}\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
&+\langle\mathcal{I}_{\Phi_{div}}\left(\mathcal{C}({{\boldsymbol{u}}})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}.\end{split}$

Using the Cauchy-Schwarz inequality, Theorem 3.1 and the fact that the
operator $\mathcal{C}$ is bounded and linear, the first term can easily be
estimated by

$\begin{split}\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}&\leq\|\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
Ch^{\sigma}_{\mathcal{X}}\|\mathcal{C}({{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&\leq
Ch^{\sigma}_{\mathcal{X}}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}.\end{split}$

Again, Cauchy-Schwarz and the boundedness of the operator $\mathcal{C}$ yield

$\begin{split}\langle\mathcal{I}_{\Phi_{div}}\left(\mathcal{C}({{\boldsymbol{u}}})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right),{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\leq&C\|{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}},\end{split}$

where the first norm on the right hand side can be estimated by (23). ∎

We will now prove a local error estimate for our semi-discrete solution. Note
that we will write ${{\boldsymbol{u}}}(t)={{\boldsymbol{u}}}(\cdot,t)$ for the
sake of readability.

######  Lemma 5.4.

Let $\sigma>1$, ${{\boldsymbol{f}}}\in
L^{\infty}(0,T_{h};{{\boldsymbol{H}}}^{\sigma})$ and ${{\boldsymbol{u}}}(t)\in
C^{1}([0,T_{h}),{{\boldsymbol{H}}}^{2\sigma+1})$. Let $M\geq 1$ such that

$T_{M}:=\sup\\{t\in[0,T_{h})\;|\;\max_{s\in[0,t]}\|{{\boldsymbol{u}}}(\cdot,s)-{{\boldsymbol{u}}}_{h}(\cdot,s)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
M\\}$

exists. Then, there exist two constants
$C_{1}=C_{1}({{\boldsymbol{u}}},{{\boldsymbol{f}}},\nu),C_{2}=C_{2}({{\boldsymbol{u}}},{{\boldsymbol{f}}},\nu)>0$,
independent of $\mathcal{X}$ and $M$, such that

$\sup_{t\in[0,T_{M})}\\!\left(\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\!\\!+\\!\nu
h_{\mathcal{X}}\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\right)\\!\leq\\!C_{1}M^{2}e^{C_{2}M^{2}T_{M}}h_{\mathcal{X}}^{\sigma}.$

###### Proof.

Using Lemma 5.1 with $\chi={\boldsymbol{\theta}}_{h}$ yields

$\begin{split}\frac{1}{2}\partial_{t}\|{\boldsymbol{\theta}}_{h}\|^{2}_{\Phi_{div}}-\nu\langle{\boldsymbol{\Delta}}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}=&\left\langle-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}+{{\boldsymbol{f}}}-\mathcal{I}_{\Phi_{div}}f,{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}}\\\
&-\left\langle\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}}\\\
&-\left\langle\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi_{div}}\mathcal{C}({{\boldsymbol{u}}}_{h}),{\boldsymbol{\theta}}_{h}\right\rangle_{\Phi_{div}},\end{split}$

where the first term can be bounded by

$\begin{split}\langle-\partial_{t}{\boldsymbol{\rho}}_{h}+\nu{\boldsymbol{\rho}}_{h}&+{{\boldsymbol{f}}}-\mathcal{I}_{\Phi_{div}}{{\boldsymbol{f}}},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
\leq&\left(\|\partial_{t}{\boldsymbol{\rho}}_{h}\|_{\Phi_{div}}+\nu\|{\boldsymbol{\rho}}_{h}\|_{\Phi_{div}}+\|{{\boldsymbol{f}}}-\mathcal{I}_{\Phi_{div}}{{\boldsymbol{f}}}\|_{\Phi_{div}}\right)\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}.\end{split}$

The first two norms can be estimated by Theorem 3.2, while the third one can
be estimated by Theorem 3.1. Together with Lemma 5.2 and Lemma 5.3 we can
derive the estimate

$\begin{split}\frac{1}{2}\partial_{t}\|{\boldsymbol{\theta}}_{h}\|^{2}_{\Phi_{div}}-\nu\langle{\boldsymbol{\Delta}}^{*}&{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\\
\leq&Ch_{\mathcal{X}}^{\sigma}\left(\|\partial_{t}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}+\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}+\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}\right)\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\\\
&+Ch_{\mathcal{X}}^{\sigma}M\left(1+\|{{\boldsymbol{u}}}\|_{H^{2\sigma+1}_{div}}\right)\left(\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}+\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\right)\\\
&+C\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}+CM\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\\\
&+Ch_{\mathcal{X}}^{\sigma}\|{{\boldsymbol{u}}}\|_{H^{2\sigma}_{div}}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}+C\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}\end{split}$

We will use the Cauchy-Schwarz inequality to separate the terms including
${\boldsymbol{\theta}}_{h}$. Especially, using (22), we will choose a constant
$c_{\nu}>0$, depending on $\nu$, such that

$c_{\nu}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}^{2}\\!<\\!-\frac{\nu}{2}\langle{\boldsymbol{\Delta}}^{*}{\boldsymbol{\theta}}_{h},\\!{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}.$

Collecting all terms including
$\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}$ and
using the Cauchy-Schwarz inequality we can derive

$\begin{split}Ch_{\mathcal{X}}^{\sigma}M&\left(1+\|{{\boldsymbol{u}}}\|_{H^{2\sigma+1}_{div}}\right)\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}+CM\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}\\\
\leq&2C^{2}h_{\mathcal{X}}^{2\sigma}M^{2}c_{\nu}^{-1}\left(1+\|{{\boldsymbol{u}}}\|_{H^{2\sigma+1}_{div}}\right)^{2}+C^{2}M^{2}c_{\nu}^{-1}\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}+c_{\nu}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}^{2}\end{split}$

Overall, we arrive at the estimate

$\begin{split}\frac{1}{2}\partial_{t}\|{\boldsymbol{\theta}}_{h}\|^{2}_{\Phi_{div}}\\!\\!-\\!\nu\langle{\boldsymbol{\Delta}}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\\!\leq&CM^{2}h_{\mathcal{X}}^{2\sigma}\\!\left(\|\partial_{t}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}+\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}\right)\\\
&+CM^{2}h_{\mathcal{X}}^{2\sigma}\\!\left(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma+1}_{div}}^{2}+\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}\right)\\\
&+C(1+M)\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}+c_{\nu}\|{\boldsymbol{\theta}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}},\end{split}$

Rearranging this inequality yields

(25)
$\begin{split}\frac{1}{2}\partial_{t}\|{\boldsymbol{\theta}}_{h}\|^{2}_{\Phi_{div}}-\frac{\nu}{2}\langle{\boldsymbol{\Delta}}^{*}{\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h}\rangle_{\Phi_{div}}\leq&CM^{2}h_{\mathcal{X}}^{2\sigma}\left(\|\partial_{t}{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}+\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}\right)\\\
&+CM^{2}h_{\mathcal{X}}^{2\sigma}\left(\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma+1}_{div}}^{2}+\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}}^{2}\right)\\\
&+C(1+M)\|{\boldsymbol{\theta}}_{h}\|_{\Phi_{div}}^{2}.\end{split}$

Hence, for every $t\in[0,T_{M})$, Gronwall’s inequality implies

$\begin{split}\|{\boldsymbol{\theta}}_{h}(t)\|^{2}_{\Phi_{div}}\leq&\|{\boldsymbol{\theta}}_{h}(0)\|^{2}_{\Phi_{div}}+CM^{2}e^{CM^{2}t}h_{\mathcal{X}}^{2\sigma}\|{{\boldsymbol{u}}}\|_{H^{1}(0,T_{M};{{\boldsymbol{H}}}^{2\sigma}_{div})}^{2}\\\
&+CM^{2}e^{CM^{2}t}h_{\mathcal{X}}^{2\sigma}\left(\|{{\boldsymbol{u}}}\|_{L^{2}(0,T_{M};{{\boldsymbol{H}}}^{2\sigma+1}_{div})}^{2}+\|{{\boldsymbol{f}}}\|_{L^{2}(0,T_{M};{{\boldsymbol{H}}}^{2\sigma}_{div})}^{2}\right).\end{split}$

The initial error can be bounded by

$\|{\boldsymbol{\theta}}_{h}(0)\|^{2}_{\Phi_{div}}\leq\|\mathcal{R}_{\Phi_{div}}{{\boldsymbol{u}}}_{0}-{{\boldsymbol{u}}}_{0}\|^{2}_{\Phi_{div}}+\|{{\boldsymbol{u}}}_{0}-\mathcal{I}_{\Phi_{div}}{{\boldsymbol{u}}}_{0}\|^{2}_{\Phi_{div}}\leq
Ch_{\mathcal{X}}^{\sigma}\|{{\boldsymbol{u}}}_{0}\|_{{{\boldsymbol{H}}}^{2\sigma}_{div}},$

which finishes the proof for the first part. For the second part we can use
(22) again and insert the just derived estimation in (25) and integrate with
respect to the time variable. This yields the same error bound as above for
$\|{\boldsymbol{\theta}}_{h}\|_{L^{2}(0,T_{M};{{\boldsymbol{H}}}^{\sigma+1}_{div})}$.
For the projection error ${\boldsymbol{\rho}}_{h}$, the approximation order is
unfortunately reduced by one because the error is measured in the
${{\boldsymbol{H}}}^{\sigma+1}_{div}$ norm. ∎

Now, we will prove our main result of this paper. As the local in time
solution has been given in Lemma 5.4, it remains to show that the existence of
the semidiscrete solution only depends on the existence interval of the
solution of the projected Navier-Stokes equations. If the conditions of Lemma
4.1 are satisfied and this solution exists global-in-time, our approximation
will exist global in time as well.

######  Theorem 5.5.

Let $\sigma>3$. Suppose the initial velocity satisfies
${{\boldsymbol{u}}}_{0}\in{{\boldsymbol{H}}}^{2\sigma+1}_{div}$ and the right
hand side satisfies ${{\boldsymbol{f}}}\in
C([0,\infty);{{\boldsymbol{H}}}^{2\sigma})$. Let ${{\boldsymbol{u}}}$ be the
global solution of the projected Navier-Stokes equations (15) - (16) on a
maximum time interval $[0,T]$ for a $T<\infty$ or on $[0,\infty)$. Moreover,
let $\mathcal{X}\subset{\mathbb{S}}^{d}$ be a finite set of pair-wise distinct
points with fill distance $h_{\mathcal{X}}\leq h_{0}$ for some $h_{0}>0$.

Then, the semi-discrete solution ${{\boldsymbol{u}}}_{h}\in
C^{1}([0,T_{h});\mathcal{F}_{\Phi_{div}}(\mathcal{X}))$ of the discretised
Navier-Stokes equations exists in all of $[0,T]$ or $[0,\infty)$,
respectively, and there exist constants $C_{1},C_{2}>0$, depending on
${{\boldsymbol{f}}}$, $\nu$, ${{\boldsymbol{u}}}_{0}$ and
${{\boldsymbol{u}}}$, but independent of $\mathcal{X}$, such that

$\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+\nu
h_{\mathcal{X}}\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma+1}_{div}}<C_{1}e^{C_{2}t}h_{\mathcal{X}}^{\sigma}.$

for all $t\in[0,T]$ or $t\in[0,\infty)$, respectively.

###### Proof.

Since the initial error is just the interpolation error of
${{\boldsymbol{u}}}_{0}$, it can be bounded by $Ch_{\mathcal{X}}^{\sigma}$.
Now choose $h_{1}$ small enough such that
$\|{{\boldsymbol{u}}}(0)-{{\boldsymbol{u}}}_{h}(0)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}<1$
for all $\mathcal{X}\subset{\mathbb{S}}^{2}$ with $h_{\mathcal{X}}<h_{1}$ and
fix $M:=2$. Since ${{\boldsymbol{u}}}$ and ${{\boldsymbol{u}}}_{h}$ are
continuous in time, we have that

$T_{M}=\sup\left\\{t\in[0,T_{h})\;|\;\max_{s\in[0,t]}\|{{\boldsymbol{u}}}(s)-{{\boldsymbol{u}}}_{h}(s)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
2\right\\}$

exists with $T_{M}>0$. Applying Lemma 5.4 yields that there exists a constant
$C>0$ such that

$\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
Ch_{\mathcal{X}}^{\sigma}$

for all $t\in[0,T_{M})$. Hence, we find a $0<h_{0}<h_{1}$ such that
$\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
1$ for all $t\in[0,T_{M})$ and $h_{\mathcal{X}}<h_{1}$.

To complete the proof, it remains to show that $T_{M}=T$ or $T_{M}=\infty$,
respectively. We will concentrate on the second case. The proof for $T_{M}=T$
works analogously.

First, suppose that $T_{M}<T_{h}$. Since ${{\boldsymbol{u}}}$ and
${{\boldsymbol{u}}}_{h}$ are continuous in time, we know that
$\|{{\boldsymbol{u}}}(T_{M})-{{\boldsymbol{u}}}_{h}(T_{M})\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
1$. Moreover, there exists an $\epsilon>0$ such that
$T_{M}+\epsilon\in[0,T_{h})$ and

$\|{{\boldsymbol{u}}}(T_{M})-{{\boldsymbol{u}}}_{h}(T_{M})-({{\boldsymbol{u}}}(T_{M}+t)-{{\boldsymbol{u}}}_{h}(T_{M}+t))\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
1$

for all $t\in[0,\epsilon)$. But then we have
$\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\leq
2$ for all $t\in[0,T_{M}+\varepsilon)$, which is a contradiction to the
definition of $T_{M}$. Hence, $T_{M}\geq T_{h}$.

Moreover, ${{\boldsymbol{u}}}_{h}$ is bounded on
${\mathbb{S}}^{2}\times[0,T_{h})$ since we have

$|u_{h}({{\boldsymbol{x}}},t)|\leq\|{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\!\leq\|{{\boldsymbol{u}}}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\!+\|{{\boldsymbol{u}}}(t)-{{\boldsymbol{u}}}_{h}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}\\!\leq\|{{\boldsymbol{u}}}(t)\|_{{{\boldsymbol{H}}}^{\sigma}_{div}}+M$

for all ${{\boldsymbol{x}}}\in{\mathbb{S}}^{2}$ and $t\in[0,T_{h})$. But
according to the continuity criterion for ODEs, this would mean that our
numerical solution would exists globally in time. Hence, we have $T_{M}\geq
T_{h}=\infty$. ∎

###  5.2. Calculation of the pressure

To obtain an approximation for the pressure $p$ at least up to a constant, we
have to solve the curl-free projection of the Navier-Stokes equations (17).
The approximated version of this equation is given by

(26)
$\nabla_{*}p_{h}=(\operatorname{id}-\mathcal{P})\mathcal{I}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right),\quad\text{in
}{\mathbb{S}}^{2}\times(0,T).$

Fortunately, we can obtain this approximation without much effort. Since we
already need to calculate the right hand side of (19), we already have time-
dependent coefficients
${\boldsymbol{\beta}}_{j}(t)={\boldsymbol{\beta}}_{j}\in
T_{{{\boldsymbol{x}}}_{j}}({\mathbb{S}}^{2})$, $1\leq j\leq N$, such that

$\mathcal{P}\mathcal{I}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right)=\sum_{j=1}^{N}\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}.$

Just using these coefficients, we can give an approximation to the pressure by
simply exchanging the kernel function by its curl-free part and using (11) to
arrive at

$\begin{split}(\operatorname{id}-\mathcal{P})\mathcal{I}_{\Phi}\left({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h})\right)&=\sum_{j=1}^{N}\Phi_{curl}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}\\\
&=\nabla^{*}\sum_{j=1}^{N}(\nabla^{*})^{T}\phi(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}.\end{split}$

Hence an approximation of the pressure $p$, up to a constant, is simply given
by

$p_{h}=\sum_{j=1}^{N}(\nabla^{*})^{T}\phi(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}.$

Since the integral over the surface-gradient of
$\phi(\cdot,{{\boldsymbol{x}}}_{j})\in H^{\sigma}$ vanishes for every $1\leq
j\leq N$, see [10], the mean integral of our approximation $p_{h}$ vanishes,
too. Therefore, we assume that the integral over the true pressure p also
vanishes. The following error estimate for the approximation of the pressure
can be established.

######  Theorem 5.6.

Under the assumptions of Theorem 5.5 and that the meant integral of $p$
vanishes, there exists a constant $C>0$ such that

$\|p-p_{h}\|_{L^{2}(0,T;\\!H^{\sigma+1})}\leq Ch_{\mathcal{X}}^{\sigma-1}$

holds for all data sites $\mathcal{X}\in{\mathbb{S}}^{2}$ with fill distance
$h_{\mathcal{X}}\leq h_{0}$ and all $t\in[0,T)$.

###### Proof.

Since both $p$ and $p_{h}$ have zero mean integral, we have

$\|p-p_{h}\|_{H^{\sigma+1}}\leq\sqrt{2}\|\nabla^{*}p-\nabla^{*}p_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}.$

Inserting equation (17) and (26) for $\nabla^{*}p$ and $\nabla^{*}p_{h}$,
respectively, and splitting up the singe terms yields

$\begin{split}\|p-p_{h}\|_{H^{\sigma+1}}\leq&\|{{\boldsymbol{f}}}-\mathcal{I}_{\Phi}{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}+\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
&+\|\mathcal{C}({{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{C}({{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}.\end{split}$

With Theorem 3.1, the first term can be estimated by

$\|{{\boldsymbol{f}}}-\mathcal{I}_{\Phi}{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\leq
ch^{\sigma}_{\mathcal{X}}\|{{\boldsymbol{f}}}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}.$

The second term can be splitted up to

$\begin{split}\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\leq&\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
&+\|\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}},\end{split}$

where the first part can simply be bounded by

$\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{I}_{\Phi}\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\leq
ch_{\mathcal{X}}^{\sigma}\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{2\sigma}_{curl}}\leq
ch_{\mathcal{X}}^{\sigma}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{curl}}\|{{\boldsymbol{u}}}\|_{H^{2\sigma+1}_{curl}}.$

For the second part, we use that the operator $\mathcal{I}_{\Phi}$ is bounded
in ${{\boldsymbol{H}}}^{\sigma}_{curl}$. Hence, we have

$\begin{split}\|\mathcal{I}_{\Phi}(\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h}))\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\leq&\|\mathcal{B}({{\boldsymbol{u}}},{{\boldsymbol{u}}})-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
\leq&\|\mathcal{B}({{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
&+\|\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
\leq&C\|{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma+1}_{curl}}.\end{split}$

The last term can be split and estimated like the previous one. We obtain

$\begin{split}\|\mathcal{C}({{\boldsymbol{u}}})\\!-\mathcal{I}_{\Phi}\mathcal{C}({{\boldsymbol{u}}}_{h})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\leq&\|\mathcal{C}({{\boldsymbol{u}}})\\!-\mathcal{I}_{\Phi}\mathcal{C}({{\boldsymbol{u}}})\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\!\\!+\\!\|\mathcal{I}_{\Phi}\left(\mathcal{C}({{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h})\right)\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\\\
\leq&C\left(h_{\mathcal{X}}^{\sigma}\|{{\boldsymbol{u}}}\|_{{{\boldsymbol{H}}}^{2\sigma}_{curl}}+\|{{\boldsymbol{u}}}-{{\boldsymbol{u}}}_{h}\|_{{{\boldsymbol{H}}}^{\sigma}_{curl}}\right).\end{split}$

Inserting the results from Theorem 5.5, squaring and integrating over time
finishes the proof. ∎

##  6\. Numerical considerations

In this section we will rewrite the spatial-discretized Navier-Stokes
equations (19) - (20) into a time-depended linear system, which can then be
solved using standard ODE solvers. Moreover, we take a look at the numerical
effort to show that the given method can be evaluated efficiently. In the end,
we will evaluate the new method on a standard benchmark problem, see [6, 9,
16], to show the stability of the new method.

###  6.1. Algorithm and Cost

Using the representation of the approximated solution (18), we can write the
discretised Navier-Stokes equations (19) - (20) in matrix-vector form given by

(27) $\displaystyle A_{\Phi_{div}}\dot{{\boldsymbol{\alpha}}}$
$\displaystyle=\nu
A_{\Delta^{*}\Phi_{div}}{\boldsymbol{\alpha}}+{{\boldsymbol{g}}}({\boldsymbol{\alpha}})\quad\text{on
}(0,T),$ (28) $\displaystyle A_{\Phi_{div}}{\boldsymbol{\alpha}}(0)$
$\displaystyle=({{\boldsymbol{u}}}_{0}({{\boldsymbol{x}}}_{k}))_{1\leq k\leq
N},$

where the kernel-based $3N\times 3N$-matrices are given by

$A_{\Phi_{div}}=(\Phi_{div}({{\boldsymbol{x}}}_{k},{{\boldsymbol{x}}}_{j}))_{1\leq
j,k\leq N},\quad
A_{\Delta^{*}\Phi_{div}}=(\Delta^{*}\Phi_{div}({{\boldsymbol{x}}}_{k},{{\boldsymbol{x}}}_{j}))_{1\leq
j,k\leq N}.$

Note that both matrices $A_{\Phi_{div}}$ and $-A_{\Delta^{*}\Phi_{div}}$ are
positive definite on a $2N$-dimensional subspace of ${\mathbb{R}}^{3N}$ since
the kernel function $\Phi_{div}$ as well as $-\Delta^{*}\Phi_{div}$ are
positive definite on the tangent space of the sphere. To solve the $3N\times
3N$ linear systems (27) and (28) in the implementation, we can reduce it to a
$2N\times 2N$ system by introducing coordinates and bases for the various
tangent planes of the sphere. The corresponding matrices of the reduced system
are classically positive definite. For details and the exact implementation,
see [15].

The right hand side term
${{\boldsymbol{g}}}:{\mathbb{R}}^{3N}\to{\mathbb{R}}^{3N}$ is given by

(29)
${{\boldsymbol{g}}}({\boldsymbol{\alpha}}):=\left(\mathcal{P}_{\Phi}({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h}))({{\boldsymbol{x}}}_{k})\right)_{1\leq
k\leq N},$

where ${{\boldsymbol{u}}}_{h}$ is given by (18). To calculate the value of
${{\boldsymbol{g}}}$ for a given ${\boldsymbol{\alpha}}$, we have to solve the
interpolation problem

$\sum_{j=1}^{N}\Phi({{\boldsymbol{x}}}_{k},{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}=({{\boldsymbol{f}}}-\mathcal{B}({{\boldsymbol{u}}}_{h},{{\boldsymbol{u}}}_{h})-\mathcal{C}({{\boldsymbol{u}}}_{h}))({{\boldsymbol{x}}}_{k}),\quad
1\leq k\leq N.$

The function ${{\boldsymbol{g}}}$ is then given by
${{\boldsymbol{g}}}({\boldsymbol{\alpha}})=\sum_{j=1}^{N}\Phi_{div}(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\beta}}_{j}$.

To solve the received ODE, standard explicit or semi-implicit methods can be
used. However, using explicit solvers, one has to be aware that the time step
$\tau$ is probably coupled to the filling distance $h_{\mathcal{X}}$ and has
to satisfy some kind of CFL condition for the system to remain stable.
however, for a investigation of the numerical cost of the new method, we will
neglect the fact and restrict ourselves to two simple examples for the time
discretization, the explicit and the semi-implicit Euler algorithm. Let
$\tau>0$ be the time discretization parameter. Then, a time step for an
explicit Euler scheme is given by

(30) ${\boldsymbol{\alpha}}_{n+1}={\boldsymbol{\alpha}}_{n}+\tau
A_{\Phi_{div}}^{-1}\left(\nu
A_{\Delta^{*}\Phi_{div}}{\boldsymbol{\alpha}}_{n}+{{\boldsymbol{g}}}({\boldsymbol{\alpha}}_{n})\right)$

and a time step of the semi-implicit Euler scheme is given by

(31)
$\begin{split}{\boldsymbol{\alpha}}_{n+1}&=\tau\left(A_{\Phi_{div}}-\nu\tau
A_{\Delta^{*}\Phi_{div}}\right)^{-1}\left(A_{\Phi_{div}}{\boldsymbol{\alpha}}_{n}+{{\boldsymbol{g}}}({\boldsymbol{\alpha}}_{n})\right)\\\
&={\boldsymbol{\alpha}}_{n}+\tau\left(A_{\Phi_{div}}-\nu\tau
A_{\Delta^{*}\Phi_{div}}\right)^{-1}\left(\nu
A_{\Delta^{*}\Phi_{div}}{\boldsymbol{\alpha}}_{n}+{{\boldsymbol{g}}}({\boldsymbol{\alpha}}_{n}).\right)\end{split}$

The costs of the given schemes depend highly on the ability to solve the two
linear systems given in (29) and in (30) or (31), respectively. This would
have a cost of $\mathcal{O}(N^{3})$ for each linear system using standard
solvers in each time-step. However these, costs can be reduced by various
strategies. Since the matrices are constant in time, one can take a Cholesky
decomposition on the start of the cost of $\mathcal{O}(N^{3})$. Then, the cost
of the solution of the linear system in each time-step reduces to
$\mathcal{O}(N^{2})$. Instead of a Cholesky decomposition, it is also possible
to use Lagrange basis functions, with which it is theoretically possible to
reduce the effort of one timestep to $\mathcal{O}(N^{2})$. However, as in the
case of Cholesky decomposition, it we need to calculate the Lagrange basis
function at the beginning of the cost of $\mathcal{O}(N^{3})$. We will give a
little preview about the Lagrange basis in the conclusion. Note that both, the
Cholesky decomposition as well as the Lagrange functions have to be calculated
once for a given point set, which means that we can use the same decomposition
for different initial conditions or external forces ${{\boldsymbol{f}}}$.

Another possibility is to use the fact that the kernel functions have compact
support. By using scaled kernel functions of the form
$\phi_{\epsilon}=\phi(\cdot/\epsilon)$ for an $\epsilon>0$, the matrices
appearing in the schemes above can be made sparse for small $\epsilon$. In
this case, the effort can be reduced to $\mathcal{O}(N^{2}\log(N))$ by using
special sparse solvers. However, it is to be expected that the theoretical
convergence rate will depend on the scaling parameter $\epsilon$, which means
that $\epsilon$ has to be coupled to the fill distance $h_{\mathcal{X}}$ in a
certain way in order to still have convergence for $h_{\mathcal{X}}\to 0$ and
$\epsilon\to 0$.

Of course, any other time discretization scheme would be possible. We refer to
general IMEX (Implicit-Explicit) schemes, see [2], which have better stability
properties than explicit methods. Furthermore, it is also possible to give
fully implicit schemes, see [21]. Unfortunately, those schemes are highly
inefficient since in every time-step a non-linear equation system has to be
solved.

###  6.2. Numerical Validation

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=10

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=20

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=30

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=40

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=50

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=60

Figure 1. Velocity field ${{\boldsymbol{u}}}_{h}$ at various times.

We will now demonstrate the stability of the new method on a test case which
is motivated by the numerical demonstration in [9], see also [6, 16]. This
simulation should only be seen as a numerical test to validate that the new
method is stable. A proper numerical validation of the results in chapter 3 is
unfortunately not possible, since the author is not aware of any analytical
solution of the Navier-Stokes equations with which a numerical solution of the
method could be compared.

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=10

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=20

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=30

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=40

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=50

Velocity field ${{\boldsymbol{u}}}_{h}$ at time t=60

Figure 2. Pressure $p_{h}$ at various times.

The initial velocity ${{\boldsymbol{u}}}_{0}$, which is chosen in analogy to
[9], is given by ${{\boldsymbol{u}}}_{0}={{\boldsymbol{L}}}^{*}p$, where
$p({{\boldsymbol{x}}})=\sum_{|{{\boldsymbol{n}}}|<20}a_{{{\boldsymbol{n}}}}{{\boldsymbol{x}}}^{{\boldsymbol{n}}}$
is a polynomial of degree $19$ in ${\mathbb{R}}^{3}$ such that
$\|{{\boldsymbol{u}}}_{0}\|_{{{\boldsymbol{L}}}^{2}}=1$. The coefficients
$a_{{{\boldsymbol{n}}}}$ were generated by a random function such that
$a_{{{\boldsymbol{n}}}}\sim|{{\boldsymbol{n}}}|^{-1}$. This construction
ensures that the initial velocity is a tangential and divergence-free velocity
field. The time-depended external force $f$ is given by

$f(t,{{\boldsymbol{x}}})=\gamma(t){{\boldsymbol{y}}}_{3,0}({{\boldsymbol{x}}})$

with

$\gamma(t)=\begin{cases}1&\text{for }t\in[0,10],\\\
\cos\left(\frac{\pi}{5}t\right)e^{(10-t)/5}&\text{for
}t\in]10,60].\end{cases}$

Note that the initial condition as well as the forcing term both satisfy the
conditions of Lemma 4.1, such that an analytical solution of the problem
exists global in time. The viscosity coefficient and the angular velocity are
set to $\nu=10^{-4}$ and $\Omega=1$. The set of data points $\mathcal{X}$ is
given by $N=2500$ points which are distributed by minimizing the potential
energy of the points, see [30]. For the basis kernel function, we used the
unscaled Wendland function $\phi_{4}$ from Table 1.

As the time integration scheme, we used a third order implicit-explicit Runge-
Kutta method which can be found in [2, Section 2.4] with a timestep
$\tau=10^{-2}$. We used the explicit solver in the IMEX scheme to handle the
term ${{\boldsymbol{g}}}$ in (27) and the implicit solver for the Laplace
term.

The approximated velocity field ${{\boldsymbol{u}}}_{h}$ as well as the
approximated pressure $p_{h}$ can be seen in Figure 1 and Figure 2,
respectively, for $t=10,20,30,40,50,60$. The approximated
${{\boldsymbol{L}}}^{2}$-norms of the velocity and the pressure are given by

$e_{{\boldsymbol{u}}}(t)=\sqrt{\frac{4\pi}{N}\sum_{j=1}^{N}\|{{\boldsymbol{u}}}_{h}({{\boldsymbol{x}}}_{j},t)\|_{2}^{2}},\quad
e_{p}(t)=\sqrt{\frac{4\pi}{N}\sum_{j=1}^{N}|p_{h}({{\boldsymbol{x}}}_{j},t)|^{2}},\quad
t\in[0,60].$

The time evolution of both norms can be seen in Figure 3.

Figure 3. Time evolution of $e_{{\boldsymbol{u}}}$ (left, red) and $e_{p}$
(right, red) in comparison with the norm of the external force
${{\boldsymbol{f}}}$ (blue).

As one can see, the solution remains stable in the whole time interval
$[0,60]$. As in [9] and [16], the initial random velocity evolves into a flow
with large structures. Also the pressure develops large, smooth structures
over the sphere . In Figure 3, one can see that $e_{{\boldsymbol{u}}}$ as well
as $e_{p}$ are increasing in the time interval $[0,10]$, where the norm of the
forcing term is constant at its maximum. As the external force decreases for
$t>10$, the damping effect of the diffusion term becomes dominant and
$e_{{\boldsymbol{u}}}$ as well as $e_{p}$ decay slowly in time, as it can also
be seen in [16]. The results shown are thus in qualitative agreement with
those of [9] and [16].

##  7\. Conclusion and Remarks

In this paper, we developed a new kernel-based collocation method for the
incompressible Navier-Stokes equations on the two dimensional sphere. We
proved convergence of the method based on the fill distance of the collocation
points and the smoothness of the initial data for both the velocity and the
pressure of the system.

The advantages of this method are various. The new method is a meshfree
method, which means that the positions of the data points $\mathcal{X}$ are
independent of any underlying grid on the sphere. Moreover, it needs no
quadrature formula for the computation, which makes it less vulnerable to
additional errors. The method leads to high-order analytically divergence-free
approximations of the velocity field, which can be computed by any smoothness
depending on the choice of the kernel function. By its kernel-based character,
it is simple to give high-order approximations of the pressure by simply
exchanging the kernel function. Finally, it is at least as efficient as
comparable methods like the one of Fengler and Freeden [9] or Ganesh, Le Gia
and Sloan [16], as the complexity can be reduced to $\sim N^{2}$ calculations
per time step, where $N$ is is number of data points.

This can be achieved by using the Lagrangian basis. The Lagrange functions of
the basis $\Phi(\cdot,{{\boldsymbol{x}}}_{j})\in{\mathbb{R}}^{3\times 3}$ is
given by the functions

$\Psi_{k}=\sum_{j=1}^{N}\Phi(\cdot,{{\boldsymbol{x}}}_{j}){\boldsymbol{\gamma}}_{j}^{k},\quad
1\leq k\leq N,$

such that we have the interpolation condition
$\Psi_{k}({{\boldsymbol{x}}}_{j})=\delta_{k,j}I_{3\times 3}$ for every $1\leq
k,j\leq N$. Given a Lagrange basis, the interpolation of a tangential vector-
field ${{\boldsymbol{u}}}$ is just given by

${{\boldsymbol{u}}}_{h}({{\boldsymbol{x}}})=\sum_{j=1}^{N}\Psi(\cdot,{{\boldsymbol{x}}}_{j}){{\boldsymbol{u}}}({{\boldsymbol{x}}}_{j}),$

which means among other things that the linear systems in (30) and (31) can be
given without any further calculation. However, the evaluation of the
approximation ${{\boldsymbol{u}}}_{h}$ becomes more expensive, since the cost
of an evaluation of the function $\Psi_{k}$ is linear in $N$. Hence, the cost
of the whole method will be $\sim N^{2}$. The use of local Lagrange functions,
see for example [4], which reduces the cost of an evaluation of $\Psi$ at the
cost of an additional error, is part of current research.

However, there remains much space to improve and investigate the current
method. The work mainly involves the analysis of discretisation in space. An
exact analysis of the time discretisation, the choice of the time
discritisation scheme as well as the the time step $\tau$ and the associated
stability of the method need to be researched further.

Moreover, the numerical behaviour of the method has to be further
investigated. We only validated the method on a benchmark test case with
unknown analytical solution and on a nice node set. The method has to be
tested on some examples with known analytical solution to determine the
numerical order of convergence and to compare them with the results from
Theorem 5.5 and Theorem 5.6.

Moreover, in the presented test case, only unscaled kernels where used.
However, kernel functions of the form $\phi_{\epsilon}(x)=\phi(x/\epsilon)$
should improve the efficiency of the method since the occurring matrices will
be sparser. However, this will be achieved at the cost of accuracy, since it
is expected that the interploation error behaves at least like
$\sim\epsilon^{-\sigma}$. An investigation of the numerical error depending on
$h_{\mathcal{X}}$ and $\epsilon$ has still to be done for the presented
method.

## References

  * [1] M. S. Agranovič, _Elliptic singular integro-differential operators_ , Uspehi Mat. Nauk 20 (1965), no. 5 (125), 3–120. MR 0198017
  * [2] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri, _Implicit-explicit runge-kutta methods for time-dependent partial differential equations_ , Applied Numerical Mathematics 25 (1997), no. 2, 151 – 167, Special Issue on Time Integration.
  * [3] M. N. Benbourhim and A. Bouhamidi, _Meshless pseudo-polyharmonic divergence-free and curl-free vector fields approximation_ , SIAM J. Math. Anal. 42 (2010), no. 3, 1218–1245. MR 2644920
  * [4] M. D. Buhmann, _Multivariate cardinal interpolation with radial-basis functions_ , Constr. Approx. 6 (1990), no. 3, 225–255. MR 1054754
  * [5] Chongsheng Cao, Mohammad A. Rammaha, and Edriss S. Titi, _The Navier-Stokes equations on the rotating $2$-D sphere: Gevrey regularity and asymptotic degrees of freedom_, Z. Angew. Math. Phys. 50 (1999), no. 3, 341–360. MR 1697711
  * [6] A. Debussche, T. Dubois, and R. Temam, _The nonlinear Galerkin method: A multi-scale method applied to the simulation of homogeneous turbulent flows_ , Theoret. Comput. Fluid Dynamics 7 (1995), no. 2, 279–315.
  * [7] F. Dodu and C. Rabut, _Vectorial interpolation using radial-basis-like functions_ , vol. 43, 2002, Radial basis functions and partial differential equations, pp. 393–411. MR 1883575
  * [8] P. Farrell, K. Gillow, and H. Wendland, _Multilevel interpolation of divergence-free vector fields_ , IMA J. Numer. Anal. 37 (2017), no. 1, 332–353. MR 3614888
  * [9] M. J. Fengler and 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. Comput. 27 (2005), no. 3, 967–994. MR 2199916
  * [10] W. Freeden, T. Gervens, and M. Schreiner, _Constructive approximation on the sphere with applications to geomathematics_ , Numerical Mathematics and Scientific Computation, The Clarendon Press, Oxford University Press, New York, 1998, With applications to geomathematics. MR 1694466
  * [11] E. J. Fuselier, _Improved stability estimates and a characterization of the native space for matrix-valued RBFs_ , Adv. Comput. Math. 29 (2008), no. 3, 269–290. MR 2438345
  * [12] E. J. Fuselier, V. Shankar, and G. B. Wright, _A high-order radial basis function (RBF) Leray projection method for the solution of the incompressible unsteady Stokes equations_ , Comput. & Fluids 128 (2016), 41–52. MR 3461854
  * [13] E. J. Fuselier and G. B. Wright, _Stability and error estimates for vector field interpolation and decomposition on the sphere with RBFs_ , SIAM J. Numer. Anal. 47 (2009), no. 5, 3213–3239. MR 2551192
  * [14] Edward J. Fuselier, _Sobolev-type approximation rates for divergence-free and curl-free RBF interpolants_ , Math. Comp. 77 (2008), no. 263, 1407–1423. MR 2398774
  * [15] Edward J. Fuselier, Francis J. Narcowich, Joseph D. Ward, and Grady B. Wright, _Error and stability estimates for surface-divergence free RBF interpolants on the sphere_ , Math. Comp. 78 (2009), no. 268, 2157–2186. MR 2521283
  * [16] M. Ganesh, Q. T. Le Gia, and I. H. Sloan, _A pseudospectral quadrature method for Navier-Stokes equations on rotating spheres_ , 2011, pp. 1397–1430. MR 2785463
  * [17] A. A. Il’in, _Navier-Stokes and Euler equations on two-dimensional closed manifolds_ , Mat. Sb. 181 (1990), no. 4, 521–539. MR 1055527
  * [18] by same author, _Partially dissipative semigroups generated by the Navier-Stokes system on two-dimensional manifolds, and their attractors_ , Mat. Sb. 184 (1993), no. 1, 55–88. MR 1211366
  * [19] A. A. Il’in and A. N. Filatov, _Unique solvability of the Navier-Stokes equations on a two-dimensional sphere_ , Dokl. Akad. Nauk SSSR 301 (1988), no. 1, 18–22. MR 953596
  * [20] C. Keim, _Collocation methods for the navier-stokes equations_ , Ph.D. thesis, University of Bayreuth, Bayreuth, 2016.
  * [21] C. Keim and H. Wendland, _A high-order, analytically divergence-free approximation method for the time-dependent Stokes problem_ , SIAM J. Numer. Anal. 54 (2016), no. 2, 1288–1312. MR 3490500
  * [22] S. Lowitzsch, _Error estimates for matrix-valued radial basis function interpolation_ , J. Approx. Theory 137 (2005), no. 2, 238–249. MR 2186949
  * [23] by same author, _Matrix-valued radial basis functions: stability estimates and applications_ , Adv. Comput. Math. 23 (2005), no. 3, 299–315. MR 2136542
  * [24] C. Müller, _Spherical harmonics_ , Lecture Notes in Mathematics, vol. 17, Springer-Verlag, Berlin-New York, 1966. MR 0199449
  * [25] F. Narcowich, J. Ward, and G. Wright, _Divergence-free rbfs on surfaces_ , Journal of Fourier Analysis and Applications 13 (2007), 643–663.
  * [26] Francis J. Narcowich and Joseph D. Ward, _Scattered data interpolation on spheres: error estimates and locally supported basis functions_ , SIAM J. Math. Anal. 33 (2002), no. 6, 1393–1410. MR 1920637
  * [27] D. Schräder and H. Wendland, _A high-order, analytically divergence-free discretization method for Darcy’s problem_ , Math. Comp. 80 (2011), no. 273, 263–277. MR 2728979
  * [28] by same author, _An extended error analysis for a meshfree discretization method of Darcy’s problem_ , SIAM J. Numer. Anal. 50 (2012), no. 2, 838–857. MR 2914288
  * [29] Y. N. Skiba, _Mathematical problems of the dynamics of incompressible fluid on a rotating sphere_ , Springer, Cham, 2017. MR 3645447
  * [30] I. H. Sloan and R. S. Womersley, _Extremal systems of points and numerical integration on the sphere_ , Adv. Comput. Math. 21 (2004), no. 1-2, 107–125. MR 2065291
  * [31] R. Temam, _Navier-Stokes equations. Theory and numerical analysis_ , North-Holland Publishing Co., Amsterdam-New York-Oxford, 1977, Studies in Mathematics and its Applications, Vol. 2. MR 0609732
  * [32] H. Wendland, _Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree_ , Advances in Computational Mathematics 4 (1995), no. 1, 389–396.
  * [33] by same author, _Scattered data approximation_ , Cambridge Monographs on Applied and Computational Mathematics, vol. 17, Cambridge University Press, Cambridge, Cambridge, 2005. MR 2131724
  * [34] by same author, _Divergence-free kernel methods for approximating the Stokes problem_ , SIAM J. Numer. Anal. 47 (2009), no. 4, 3158–3179. MR 2551162