
\documentclass[reqno]{amsart}
\usepackage{amssymb,graphicx}

\AtBeginDocument{{\noindent\small
{\em Electronic Journal of Differential Equations},
Vol. 2004(2004), No. 42, pp. 1--13.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or
http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2004 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}

\title[\hfilneg EJDE-2004/42\hfil Degenerate elliptic equations]
{A note on a degenerate elliptic equation
 with applications for lakes and seas}

\author[D. Bresch, J. Lemoine, \& F. Guillen-G.\hfil EJDE-2004/42\hfilneg]
{Didier  Bresch, J\'er\^ome  Lemoine, \& Francisco Gu\'{\i}llen-Gonzalez}
  % in alphabetical order


\address{Didier  Bresch\hfill\break
 Laboratoire de Mod\'elisation et de Calcul - Imag (CNRS UMR 5523),
 Universit\'e Joseph Fourier,
 38041 Grenoble cedex, France}
\email{Didier.Bresch@imag.fr}

\address{J\'er\^ome  Lemoine \hfill\break
Laboratoire de Math\'ematiques  Appliqu\'ees,
Universit\'e Blaise Pascal et C.N.R.S.,
63177 Aubi\`ere cedex, France}
\email{Jerome.lemoine@math.univ-bpclermont.fr}

\address{Francisco  Guillen-Gonzalez \hfill\break
Dpto. de Ecuaciones Diferenciales y An\'alisis Num\'erico,
Universidad de Sevilla,
Aptdo. 1160, 41080 Sevilla, Spain}
\email{guillen@us.es}

\date{}
\thanks{Submitted June 15, 2002. Published March 23, 2004.}
\subjclass[2000]{35Q30, 35B40, 76D05}
\keywords{Regularity result, degenerate elliptic equation, geophysics,
\hfill\break\indent weighted Sobolev spaces, splitting projection method}

\begin{abstract}
   In this paper, we give an intermediate regularity result
   on a degenerate elliptic equation with a weight blowing up
   on the boundary.  This kind of equations is encountoured
   when  modelling some phenomena linked to seas or lakes.
   We give some examples where such regularity is useful.
\end{abstract}

\maketitle
\newtheorem{theorem}{Theorem}[section] % theorems numbered with section #
\newtheorem{lemma}[theorem]{Lemma}

\section{Introduction} \label{sec1}

   This paper is devoted to a degenerate elliptic equation that we can
find in several models in oceanography when we consider a domain with
a depth vanishing on the shore.
   A lot of mathematical studies in oceanography assume a domain with a
strictly positive depth in order to prevent the study in weighted spaces.
   Few papers have been devoted to coefficients with degenerated behavior.

   Regularity results in weighted Sobolev spaces on degenerate elliptic
equations have been studied for instance in \cite{2,3} with a
vanishing weight on the boundary that implies no boundary
condition on the unknown.
   Here we study a degenerate elliptic equation with a weight with a
blowing up comportment on the shore.
   A $H^2$ regularity in weighted spaces is proved allowing to consider
general weights.
   We will obtain such regularity by a careful study of the weight adapting the
standard method of translation, see \cite{8}.
   For example, we use and adapt some results on weighted Sobolev
spaces that have been studied in \cite{7}.

    Section  \ref{sec2} is devoted to the regularity result related to
the degenerate elliptic equation.
    Then in Section \ref{sec3}, we  explain why this kind of equation
is important in oceanography.
    In the last section we describe precise examples where
such regularity is used.
    At first we give some examples where such regularity result is used
for existence result or error estimates that means
the planetary-geostrophic equations and the vertical geostrophic
equations.
   Then we give an example where it is used in a splitting projection method.
  We also mention that such equation is obtained from the Great Lake equation.
  We remark that this degenerate elliptic equation may be found
in an other field such as in electromagnetism with the maxwell's system,
see \cite{18}.
   Similar degenerate elliptic equation may also be encountoured for a
problem related to Saint-Venant's equations if we want to apply
Babuska-Brezzi's Inf-Sup Lemma in weighted spaces, see \cite{1}.


\section{The degenerate elliptic equation} \label{sec2}

Let $\mathcal{O}$ be a two-dimensional domain.
This section is devoted to a regularity result on the following degenerate
elliptic equation: Given $h:\mathcal{O} \to \mathbb{R}$ and $g: \mathcal{
\theta} \to \mathbb{R}$,
the problem is to find $\Psi:\mathcal{O}\to \mathbb{R}$ such that
\begin{equation}
-\nabla_{\rm x}\cdot (\frac 1h\nabla_{\rm x} \Psi) =  g
\text{in } {\mathcal{O}},
\quad \Psi= 0  \text{ on } \partial\mathcal{O}. \label{e1}
\end{equation}
Here, $h$ is a function data satisfying
\begin{gather}
 h\in W^{1,\infty}(\mathcal{O}), \quad h>0 \text{ in } \mathcal{O}, \label{e2}\\
 h({\rm x}) = \varphi(\delta({\rm x})) \text{in a neightbourhood of }
 \partial\mathcal{O} \label{e3}
\end{gather}
with $\delta({\rm x})=\mathop{\rm dist}({\rm x},\mathcal{O})$.
    Moreover we assume that
\begin{gather}
\varphi \text{ is a non decreasing Lipschitz function},  \quad \varphi(0)
= 0,
\label{e4} \\
 \exists \> c>0 \text{ such that: }\ \forall s>0,\
\Bigl|\frac{\varphi'(s)}{\varphi(s)}\Bigr|\le \frac cs \label{e5}
\end{gather}
and for all $c_1,c_2>0$, there exist $\alpha_1,\alpha_2>0$ such
that
\begin{equation}
 \forall s,r >0,\;    c_1\le\frac sr\le c_2\Longrightarrow
      \alpha_1\le \frac{\varphi(s)}{\varphi(r)}\le\alpha_2. \label{e6}
\end{equation}

\noindent{\textbf{Remark} Note that $\varphi(s)= c\,s^\alpha$, $0<\alpha<1$ and $c>0$,
satisfies the previous hypothesis. \smallskip


Now, we define the function space
\begin{equation}
H(\mathcal{O}) = \{\Psi\in L^2(\mathcal{O}) :
\frac{\nabla_{\rm x} \Psi}{h^{1/2}}\in (L^2(\mathcal{O}))^2, \Psi = 0
\text{ on } \partial\mathcal{O}\} \label{e7}
\end{equation}
endowed with the norm
$$
\|\Psi\|_{H(\mathcal{O})} =\|\frac{\nabla_{\rm x} \Psi }{ h^{1/2}}\|_{(L^2(\mathcal{O}))^2}.
$$
where $h$ is defined by \eqref{e2}--\eqref{e3}.
  We remark that $\|\cdot\|_{H(\mathcal{O})}$ is a norm since $\frac 1h \ge c >0$
and $\Psi = 0$ on $\partial\mathcal{O}$ implies that there exists $c>0$
 such that
for all $\Psi \in {H(\mathcal{O})}$,
$$
\|\Psi\|_{L^2(\mathcal{O})}
 \le c \|\frac{\nabla_{\rm x} \Psi }{ h^{1/2}}\|_{(L^2(\mathcal{O}))^2}.
$$


\begin{lemma} \label{lm1}
Let $H(\mathcal{O})$ be defined by \eqref{e7} with $h$ satisfying
\eqref{e2}--\eqref{e6}. Then
$\mathcal{D}(\mathcal{O})$ is dense in ${H(\mathcal{O})}$.
\end{lemma}

The proof of this lemma is similar to the proof of \cite[Theorem 11.2]{7}.
Therefore, we omit it.
  The main result of this paper is the following

\begin{theorem} \label{thm2}
Let $\mathcal{O}$ be a (two-dimensional) bounded domain of class
$\mathcal{C^3}$.  Let $g$ be such that $\delta h^{1/2} g\in L^2(\mathcal{
O})$
with $h$ satisfying \eqref{e2}--\eqref{e6}.
There exists a unique solution $\Psi$ of \eqref{e1} such that
$\Psi \in {H(\mathcal{O})}$ and
   $$ \|\Psi\|_{H(\mathcal{O})} \le c \|\delta h^{1/2} g\|_{L^2(\mathcal{
O})}.
   $$
  Moreover, if
\begin{equation}
h^{1/2} g \in L^2(\mathcal{O}) \label{e8}
\end{equation}
then
\begin{gather}
h^{1/2} \nabla_{\rm x}(\frac 1h\nabla_{\rm x} \Psi) \in (L^2(\mathcal{O})
)^4,
     \label{e9}\\
\|h^{1/2} \nabla_{\rm x}(\frac 1h\nabla_{\rm x} \Psi)\|_{(L^2(\mathcal{O}
))^4}
     \le c\|h^{1/2} g\|_{L^2(\mathcal{O})}
\end{gather}
with $c$ a constant depending only on the domain.
\end{theorem}

\begin{proof} \textbf{Weak solutions:}\quad
The existence and uniqueness of weak solutions of \eqref{e1} follows
from the Lax-Milgram theorem, since
\begin{align*}
     \|\Psi\|^2_{H(\mathcal{O})} = \int_{\mathcal{O}} g \Psi
  &  \le \|\delta h^{1/2} g \|_{L^2(\mathcal{O})}
       \|\frac{\Psi }{ \delta h^{1/2}}\|_{L^2(\mathcal{O})} \\
  &  \le c\| \delta h^{1/2} g\|_{L^2(\mathcal{O})}
     \|\frac{\nabla_{\rm x} \Psi }{ h^{1/2}}\|_{(L^2(\mathcal{O}))^2}.
\end{align*}
   In the previous estimate we have used Hardy's inequality in
weighted space which will be proved in Lemma \ref{lm3}.

\noindent\textbf{Regularity:}
We use the usual difference quotients (cf. Brezis \cite{8}).
The interior regularity is well known since $h \ge c(\omega) >0$ in
each $\omega \Subset \mathcal{O}$.
To obtain the regularity result up to the boundary, we define
a local diffeomorphism $T$ which preserves the normal direction.
   More precisely we define the local diffeomorphism $T: Q \to V$ by
$T(x^*,r) = (x^*,\alpha(x^*)) + r\, n(x^*,\alpha(x^*))$ for all $ (x^*,
r) \in Q$
where $\partial \mathcal{O}$ is locally the graph of a $\mathcal{C}^3$
function $\alpha$
(see Figure 1).
   This property, combined with the hypothesis \eqref{e3} of $h$, will be
strongly used in the sequel.
   We also define a cut-off function
$\theta$ such that
\begin{gather*}
 \theta \in {\mathcal{C}}^2 \text{ in } V, \cr
 \theta = 0 \text{ on }\mathbb{R}^2\backslash  V, \quad
   \theta = 1 \text{ in } \mathcal{V}_1, \quad  \cr
 \frac{\partial\theta}{ \partial \widetilde n} = 0 \text{ on a neightbourhood of }
    \partial\mathcal{O}\cap V, \cr
\end{gather*}
where $\mathcal{V}_1\Subset V$ and the extension $\widetilde n$ of the
normal $n$ is
defined for all $(x,y)\in V$ by
   $$ \widetilde n(x,y)= n(T(x^*,0))$$
where $(x,y)= T(x^*,r)$, $(x^*,r)\in Q$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width= 0.6\textwidth]{fig1.eps}
\end{center}
\caption{The local diffeomorphism $T$}
\end{figure}

\noindent\textbf{remark}
Due to the $\mathcal{C}^3$ regularity of $\partial\mathcal{O}$, $T$ is a
$\mathcal{C}^2$ diffeormorphism from $Q$ onto a neightbourhood of $(x_0,y_0)$
denoted by $V$. \smallskip

  Multiplying \eqref{e1} by $\theta$, and denoting $\xi = \theta\Psi$,
it follows that $\xi$ is the (unique) solution in
$H(V\cap \mathcal{O})$ of
\begin{equation}
\begin{gathered}
 -\nabla_{\rm x}\cdot(\frac 1h\nabla_{\rm x}\xi ) = f \text{ in }
   V \cap \mathcal{O}, \\
 \xi\vert_{\partial(V \cap\mathcal{O})} = 0, \\
\end{gathered} \label{e10}
\end{equation}
where $f= \theta g + {\nabla_{\rm x} h\cdot \nabla_{\rm x} \theta\over
h^2}\Psi -\frac 1h \Delta_{\rm x}\theta \Psi -{2\over h}\nabla_{\rm x}\theta
\cdot\nabla_{\rm x}\Psi$ (since $h^{1/2}f\in L^2(V\cap\mathcal{O})$).
   For this, it suffices to check that
   $$
{\nabla_{\rm x} h\cdot \nabla_{\rm x} \theta\over h^{3/2}}\Psi \in L^2(V\cap\mathcal{O}).
$$
  Indeed, since ${\partial\theta/\partial\widetilde n} = 0$ on a neightbourhood of
$\partial\mathcal{O}\cap V$, then
   $$
\nabla_{\rm x} h\cdot \nabla_{\rm x} \theta
={\partial h \over \partial \widetilde \tau} {\partial \theta \over \partial \widetilde \tau}
= 0
$$
on a neightbourhood of $\partial\mathcal{O}\cap V$.
   We recall that $h$ doesn't depend on $\widetilde \tau$ (where
$\widetilde\tau$ is
defined as $\widetilde n$) using that $h$ is given by \eqref{e3} and using the
definition of $T$.

   Now we use the difference quotient technic on \eqref{e10} to deduce
the weight regularity announced in the theorem.
   Let $\varphi:\mathcal{O}\to\mathbb{R}^n$ ($n\ge 1$) be a function.
   We denote, for all $(x^*,r) \in Q_+$,
\begin{gather*}
\widetilde \varphi (x^*,r) = \varphi(T(x^*,r)),\\
a_{kl} = \sum_{j= 1}^2 \partial_j T^{-1}_k (T(x^*,r))\partial_j T^{-1}_l
        (T(x^*,r)) |\mathop{\rm Jac} T(x^*,r)|, \quad k,l= 1,2 ,\\
\hat k (x^*,r) = \widetilde f(x^*,r) |\mathop{\rm Jac} T(x^*,r)|.
\end{gather*}
Then we get that $\widetilde\xi$ is the unique solution in
$\widetilde H(Q_+)$ of:
\begin{equation}
\int_{Q_+} \sum_{k,l} {a_{kl}\over \widetilde h} \,\partial_k\widetilde \xi
      \,\partial_l\widetilde\varphi dx^*dr = \int_{Q_+} \hat k \widetilde
\varphi d x^*dr      \label{e11}
\end{equation}
for all $\widetilde \varphi\in \widetilde H(Q_+)$ where
$$
\widetilde H(Q_+) = \{\widetilde \varphi\in L^2(Q_+) :
\widetilde h^{-1/2}\nabla_{\rm x} \widetilde \varphi\in (L^2(Q_+))^2,
\widetilde \varphi = 0 \text{ on } \partial Q_+\}.
$$
We choose $\widetilde \varphi = D_{-\tau}(D_\tau \widetilde \xi)$ with
$\tau= |\tau| e_1$,
$D_{\tau}\widetilde \xi= \bigl(\widetilde \xi(x+\tau)-\widetilde
\xi(x)\bigr)/|\tau|$
and
$|\tau|$ small enough in order to obtain
$\widetilde \varphi \in H(Q_+)$.
   Using
$$
\|{1\over \widetilde h^{1/2}} D_{-\tau}(D_\tau \widetilde \xi) \|_{L^2(Q_+)}
      \le c \|{1\over \widetilde h^{1/2}} \nabla_{\rm x} (D_\tau
\widetilde \xi) \|_{L^2(Q_+)}
$$
and
   $$ \|\widetilde h^{1/2} \hat k \|_{L^2(Q_+)}
      \le c \| h^{1/2}g\|_{L^2(\mathcal{O})}
$$
we get
\begin{equation}
 \|\hat k D_{-\tau}(D_\tau \widetilde \xi) \|_{L^1(Q_+)}
 \le c\| h^{1/2}g\|_{L^2(\mathcal{O})}
 \|{1\over \widetilde h^{1/2}} \nabla_{\rm x} (D_\tau \widetilde\xi) \|_{
L^2(Q_+)}.
     \label{e12}
\end{equation}
  Moreover, denoting
$$
I = \sum_{k,l} \int_{Q_+} D_\tau ({1\over \widetilde h} a_{kl}
\partial_k \widetilde \xi) \partial_l(D_\tau\widetilde\xi)
$$
since $D_\tau (a_{kl}/\widetilde h)= D_\tau (a_{kl})/ \widetilde
h$, (recall that $\widetilde h$ does not depend on $\tau$) and $T\in
\mathcal{C}^2$, we have
\begin{equation}
\begin{aligned}
I& \ge c\|{1\over\widetilde h^{1/2}}\nabla_{\rm x}
(D_\tau\widetilde\xi)\|_{(L^2(Q_+))^2}^2
      - c\|{1\over\widetilde h^{1/2}} \nabla_{\rm x}  \widetilde\xi\|_{L^2}
      \|{1\over \widetilde h^{1/2}}\nabla_{\rm x} (D_\tau\widetilde\xi)
\|_{(L^2(Q_+))^2}
\\
& \ge c\|{1\over\widetilde h^{1/2}}\nabla_{\rm x}
(D_\tau\widetilde\xi)\|_{(L^2(Q_+))^2}^2
      -  c\| h^{1/2}g\|_{L^2(\mathcal{O})}\|
        {1\over \widetilde h^{1/2}}\nabla_{\rm x} (D_\tau\widetilde\xi)
      \|_{(L^2(Q_+))^2}
\end{aligned} \label{e13}
\end{equation}
   Using the variational formulation satisfied by $\widetilde
\xi$, \eqref{e12} and \eqref{e13} we get
\begin{equation}
\|{1\over \widetilde h^{1/2}}\nabla_{\rm x} (D_\tau\widetilde\xi)
\|_{(L^2(Q_+))^2}
\le c\|h^{1/2}g\|_{L^2(\mathcal{O})}. \label{e14}
\end{equation}
Thus, by classical arguments,
\begin{equation}
{\partial_1^2\widetilde \xi \over \widetilde h^{1/2}}\in L^2(Q_+)
 \quad\text{and}\quad
    {\partial_2\partial_1\widetilde \xi \over \widetilde h^{1/2}}\in L^2(
Q_+),
     \label{e15}
\end{equation}
and their respective norms are bounded by
$c\|h^{1/2} g\|_{L^2(\mathcal{O})}$.
  In particular,
\begin{equation}
\widetilde h^{1/2} \partial_1({\widetilde h}^{-1} \partial_1\widetilde \xi)
   \in L^2(Q_+), \quad
   \widetilde h^{1/2} \partial_1({\widetilde h}^{-1} \partial_2\widetilde\xi)
   \in L^2(Q_+) \label{e16}
\end{equation}
and their respective norms depend continuously on $h^{1/2} g$ in
$L^2(\mathcal{O})$.

   We remark that contrary to the homogeneous case, that means the
standard Laplacian operator, we have not yet the regularity
$\widetilde h^{1/2}\partial_2({\widetilde h}^{-1}\partial_1\widetilde \xi)
   \in L^2(Q_+)$.
   We will obtain such regularity using the hypothesis \eqref{e5} on $h$.
   Indeed, in the distribution sense,
\begin{equation}
  \widetilde h^{1/2}\partial_2({1\over \widetilde  h}\partial_1\widetilde\xi )
    = {-\partial_2\widetilde  h\over \widetilde  h^{3/2}}\partial_1\widetilde \xi
      +{1\over \widetilde  h^{1/2}}\partial_2\partial_1\widetilde \xi. \label{e17}
\end{equation}
  Since $\partial_1\widetilde\xi= 0$ on $\partial Q_+$ then
\eqref{e15} yields $\partial_1\widetilde\xi\in H(Q_+)$.
    Using the Hardy's inequality \eqref{e20} and Hypothesis \eqref{e5}, we
get
   $$
\int_{Q_+}\Bigl|{\partial_2\widetilde  h\over \widetilde
h^{3/2}}\partial_1\widetilde \xi\Bigr|^2
      \le c\int_{Q_+}{|\partial_1\widetilde \xi|^2\over \widetilde\delta^2\widetilde h}
      \le c\int_{Q_+}{|\partial_2\partial_1\widetilde\xi|^2\over \widetilde h}.
$$
   Thus, using the regularity \eqref{e15}, we get from \eqref{e17}
\begin{gather}
\widetilde h^{1/2}\partial_2({1\over \widetilde  h}\partial_1\widetilde
\xi)\in L^2(Q_+),   \label{e18} \\
\|\widetilde h^{1/2}\partial_2(\partial_1\widetilde \xi/ \widetilde  h)
\|_{(L^2(Q_+))^2}
      \le c\|h^{1/2}g\|_{L^2(\mathcal{O})}. \nonumber
\end{gather}
   Now we use the variational formulation \eqref{e11} satisfied by
$\widetilde \xi$
to obtain the regularity on
$\widetilde h^{1/2}\partial_2(\partial_2\widetilde \xi/\widetilde  h)$.
  We have
  $$ \bigl|\int_{Q_+} {a_{22}\over \widetilde h}
     \partial_2 \widetilde  \xi \partial_2 \Phi \Bigr|
     \le c \|h^{1/2}g\|_{L^2(\mathcal{O})}\|{1\over \widetilde
      h^{1/2}}\Phi\|_{L^2(Q_+)}$$
for all $\Phi \in \mathcal{D}(Q_+)$.
   Using now the weak regularity of $\widetilde \xi$ and
$a_{22} \ge c >0$ in $Q_+$, this gives
\begin{equation}
\widetilde h^{1/2}\partial_2({1\over \widetilde  h}\partial_2\widetilde \xi )
    \in L^2(Q_+). \label{e19}
\end{equation}
   Therefore, \eqref{e16}, \eqref{e18} and \eqref{e19} give the regularity
\eqref{e9}.
\end{proof}


\begin{lemma}[Hardy's inequality in weighted spaces] \label{lm3}
  Let $h$ satisfy \eqref{e2}--\eqref{e6} and let $\Psi\in H(\mathcal{O})$. Then
\begin{equation}
\|{\Psi\over \delta h^{1/2}}\|_{L^2(\mathcal{O})}
\le c \| {\nabla_{\rm x} \Psi\over h^{1/2}}\|_{L^2(\mathcal{O})^2}\label{e20}
\end{equation}
where $c$  depends only on $\mathcal{O}$.
\end{lemma}

\begin{figure}[t]
\includegraphics[width= 0.6\textwidth]{fig2.eps}
\caption{The local coordinates}
\end{figure}

\begin{proof}
The proof of this lemma is similar to the proof of the classical Hardy's
inequality (see for instance \cite{14}) introducing the corresponding
weight.
   By density, it suffices to consider $\Psi\in \mathcal{D}(\mathcal{O})$.
   The interior estimate is obvious.
In the local coordinates (see Figure 2), we write
\begin{align*}
&\int_{\alpha(x)}^z {|\Psi(x,y)|^2 dy \over |y-\alpha(x)|^2 \varphi(y-\alpha(x))}  \\
&  \le 2\int_{\alpha(x)}^z\Bigl(\int_y^{+\infty}
        {dt\over |t-\alpha(x)|^2\varphi(t-\alpha(x))}\Bigr)
        \Psi(x,y)\partial_y\Psi(x,y) dy\\
& \le 2\int_{\alpha(x)}^z {|\Psi(x,y)||\partial_y\Psi(x,y)|\over
|y-\alpha(x)|\varphi(y-\alpha(x))}dy \\
& \le 2\Bigl(\int_{\alpha(x)}^z {|\Psi(x,y)|^2\over
|y-\alpha(x)|^2\varphi(y-\alpha(x))}dy\Bigr)^{1/2}
       \Bigl(\int_{\alpha(x)}^z {|\partial_y \Psi(x,y)|^2\over
\varphi(y-\alpha(x))}dy\Bigr)^{1/2}.
\end{align*}
Therefore,
   $$
\int_{\alpha(x)}^z {|\Psi(x,y)|^2 \over |y-\alpha(x)|^2\varphi(y-\alpha(x))} \,dy
\le 4 \int_{\alpha(x)}^z {|\partial_y \Psi(x,y)|^2 \over \varphi(y-\alpha(x))}\,dy.
$$
   Thus, integrating with respect to $x$, we get
$$  \int_{V\cap\mathcal{O}}{|\Psi(x,z)|^2\over
 |z-\alpha(x)|^2\varphi(z-\alpha(x))}\, dz\,dx
 \le 4\int_{V\cap\mathcal{O}} {|\partial_y \Psi|^2\over
 \varphi(\xi-\alpha(x))}\,d\xi\, dx.
$$
   Since $\alpha$ is smooth enough, there exists $c>1$ such that, for all
$(x,z)\in V\cap\mathcal{O}$,
   $$\delta(x,z)\le |z-\alpha(x)|\le c\delta(x,z).
$$
Therefore, using \eqref{e5}--\eqref{e6}, we get
   $$  \int_{V\cap\mathcal{O}}{|\Psi|^2\over \delta^2\varphi(\delta)}
       \le c\int_{V\cap\mathcal{O}} {|\partial_y \Psi|^2\over \varphi(\delta)}
$$
and the result follows. \end{proof}

\section{Importance of this degenerate equation} \label{sec3}

Let us introduce the three-dimensional oceanographic domain
$$
\Omega =\{({\rm x},z)\in \mathbb{R} ^3 : {\rm x}= (x,y)\in \mathcal{O},
    -h({\rm x}) <z<0 \}
$$
with $\mathcal{O} \subset \mathbb{R} ^2$ the surface domain and
$h:\overline{\mathcal{O}} \to \mathbb{R}$ with $h>0$
in $\mathcal{O}$, the bottom function. Moreover,
$\Gamma_s =\overline{\mathcal{O}} \times \{0\}$ is the  surface boundary
and $\Gamma_b =\partial \Omega \setminus\Gamma_s$ the bottom.

   We denote $\nabla=(\nabla_{\rm x}, \partial_z) $ the three dimensional gradient
vector (with $\nabla_{\rm x}=(\partial_x , \partial_y)$ the vectorial
horizontal part) and  $\Delta$ is the Laplace operator.
  We explain, in this section, why such degenerate equation naturally
appears in different models issued from oceanography when hydrostatic
pressure is assumed.
  In all these equations, the field $u=(v,w)$ and the pressure $p$
satisfy the equation
\begin{equation}
\begin{aligned}
 Lv + \nabla_{\rm x}\, p = f, \quad \partial_z p = 0 \quad\text{ in }\Omega, \\
 \partial_z w = - \mathop{\rm div}_x v \quad\text{ in } \Omega,
\end{aligned} \label{e21}
\end{equation}
and at least one of the the boundary conditions
\begin{equation}
(v,w)\cdot n_{\partial\Omega}= 0, \quad
\overline v \cdot n_{\partial \mathcal{O}} = 0 \label{e22}
\end{equation}
where we use the notation
$\overline v = \int_{-h}^0 v\, dz$.
   We note that $L$ is a certain operator (algebraic or differential), see
\eqref{e29}, \eqref{e31} or \eqref{e36}  for some examples.

\noindent\textbf{Remark}
  Of course Boundary conditions \eqref{e22} are not necessary or
sufficient to solve System \eqref{e21}.
  We have to choose other boundary conditions following the choice of
the operator $L$. \smallskip


  Integrating the divergence free equation with respect to the vertical
coordinate and using the boundary condition \eqref{e22}, part 1, we
obtain
\begin{equation}
 \nabla_{\rm x}\cdot \overline v = 0 \text{ in } \mathcal{O}.  \label{e23}
\end{equation}
If the domain is simply connected then, using \eqref{e23},
there exists a stream function $\Psi$ such that
\begin{equation}
 \overline v = \nabla_{\rm x}^\bot \Psi \text{ in } \mathcal{O},    \label{e24}
\end{equation}
where $\nabla_x^\bot$ is the $2D$ curl operator,  i.e., $(-\partial_y,\partial_x)$.
The boundary condition \eqref{e22}, part 2, gives
\begin{equation}
\Psi = 0 \quad\text{on} \partial \mathcal{O}. \label{e25}
\end{equation}
   We assume that $v$ maybe formally written as
\begin{equation}
 v = A\nabla_{\rm x} p + g_1\quad \text{ in } \Omega. \label{e26}
\end{equation}
where $A$ is a matrix function (see the examples below).
   The purpose is to obtain some regularity result on $v$.
   Integrating \eqref{e26} with respect to $z$ (taking into account that
$\partial_z p =0$ in $\Omega$), we obtain
$$
\overline v = \overline A\nabla_{\rm x} p + \overline g_1,
$$
where $\overline A= \int_{-h}^0 A$ and $\overline g_1= \int_{-h}^0 g_1=
$.
Therefore, using \eqref{e24}, we obtain
$$
\nabla^\bot_{\rm x} \Psi = \overline A \nabla_{\rm x} p + \overline g_1
$$
and thus, assuming $\overline A$ invertible
\begin{equation}
\nabla_{\rm x} p = B(\nabla_{\rm x}^\bot \Psi - \overline g_1)   \label{e27}
\end{equation}
where $B=(\overline A)^{-1}$.
  Taking the horizontal curl operator of \eqref{e27}, using that
$\nabla_{\rm x}^\bot\cdot \nabla_{\rm x} = 0$, we get
\begin{equation}
\nabla_{\rm x}^\bot \cdot (B(\nabla_{\rm x}^\bot \Psi - \overline g_1)) =
0 \text{ in } \mathcal{O}, \quad \Psi= 0 \text{ on } \partial \mathcal{=
O}.
\label{e28}
\end{equation}
  On the other-hand, \eqref{e26} and \eqref{e27} yield
   $$
\nabla_{\rm x} v = \nabla_{\rm x} \Bigl( A\bigl(B(\nabla_{\rm x}^\bot \Psi - \overline g_1)\bigr)\Bigr) + \nabla_{\rm x} g_1.
$$
   Thus the regularity of $\nabla_{\rm x} v$ depends on the regularity of
$\Psi$ and $g_1$.
   Theorem \ref{thm2} may be extended easily to more general
degenerate elliptic equations including for instance \eqref{e28}.

Now assume that $A= \mathop{\rm Id}$ then we get $B= 1/h \mathop{\rm Id}$
and therefore
$A\bigl(B(\nabla_{\rm x}^\bot \Psi )\bigr) = \nabla^\bot \Psi/h$.
   Then the regularity of $\nabla_{\rm x} v$ in $(L^2(\Omega))^4$ is
given by the regularity of $\Psi$
$$
h^{1/2}\nabla_{\rm x}(h^{-1}\nabla_{\rm x}\Psi)\in (L^2(\mathcal{O}))^4
$$
deduced from Theorem \ref{thm2}.\smallskip

  Let us give now some applications of such regularity results on the
stream function.


\section{Some applications for lakes and seas} \label{sec4}

  We consider again, in the three first examples, the three-dimensional
domain
   $$ \Omega =\{({\rm x},z)\in \mathbb{R}^3 : {\rm x}= (x,y)\in \mathcal{O} ,
    \ -h({\rm x}) <z<0 \}.$$


\subsection{Planetary geostrophic equation} % subsec4.1
   Let us consider the hyperviscous parame\-trization for the planetary
equations given by
\begin{equation}
\begin{gathered}
 \varepsilon_H v + f v^\bot + \nabla_{\rm x} p = 0, \quad \partial_z p =
+ T = 0, \\
 \mathop{\rm div}_{\rm x} v + \partial_z w = 0 , \quad
     (v,w)\cdot n_{\partial\Omega} = 0, \quad
     \overline v \cdot n_{\partial \mathcal{O}} = 0, \\
 \partial_t T - K_h \Delta_{\rm x} T - K_v \partial_z^2 T + v\cdot \nabla_{\rm x} T
     + w\partial_z T  + \gamma \,\mathcal{D} T = Q,
\end{gathered}\label{e29}
\end{equation}
with $f=(1+\beta y)$ where $\beta$ is a constant, with suitable bounda=
ry
and initial conditions and $\mathcal{D}$ a suitable fourth order
differential operator.
   This system has been studied in \cite{6} using some
relations between $T$ and $u= (v,w)$ jointly with the regularity result
proved in Theorem \ref{thm2}.
   Here, we only recall how to get the estimate on ${\rm div}_{\rm x} v$
which is necessary to obtain the existence result for \eqref{e29} (the
regularity result proved in Theorem \ref{thm2} is used there).

  Let $(v,p,T)$ be an approximate sequence of
\eqref{e29} built for instance using a Galerkin method.
   Then we get the following estimate.

\begin{theorem} \label{thm4}
   Let $\mathcal{O}$ be a bounded simply connected domain of class $\mathcal{C}^3$
and $h$ satisfied the hypothesis of Theorem 2.
   Then
$$ \|{\rm div}_{\rm x}v\|_{L^2(\Omega)} \le
    c\Bigl( \|\mathcal{D}^{1/2} T\|_{L^2(\Omega)}
       + \|\nabla_{\rm x} T\|_{(L^2(\Omega))^2}\Bigr).
$$
where $c$ depends only on $\Omega$.
\end{theorem}

   This estimate allows to pass to the limit on the approximate
solutions and to get an existence result of weak solutions for
\eqref{e29}.

\begin{proof}
    The hydrostatic equation $\partial_z p + T= 0$ reads
$p= p_s + \int_z^0 T$ with $p_s= p_s({\rm x})$.
   Then,
\begin{equation}
 \widetilde M v + \nabla_{\rm x} p_s = \nabla_{\rm x} \int_z^0 T \label{e30}
\end{equation}
with $\widetilde M=\begin{pmatrix}
  \varepsilon_H & - f       \\
  f      & \varepsilon_H  \end{pmatrix}$.
   Using the same procedure than in Section \ref{sec3}, we prove
that there exists $\Psi$ such that $\overline v = \nabla_{\rm x}^\bot
\Psi$.
   The stream function $\Psi$ satisfies a degenerate elliptic
equations similar to \eqref{e1} and therefore it is possible
to prove that
  $$
h^{1/2} \nabla_{\rm x}(\frac 1h\nabla_{\rm x} \Psi) \in (L^2(\mathcal{O}))^4
$$
using Theorem \ref{thm2} and the weak regularity satisfied by
$T$.
   This regularity result implies the required estimate on
$\mathop{\rm div}_{\rm x} v$ using \eqref{e30} and the relation
   $$
-\nabla_{\rm x} p_s =   {\widetilde M\over h}\nabla_{\rm x}^\bot \Psi
        + \frac 1h \int_{-h}^0 \nabla_{\rm x}(\int_z^0 T)
$$
coming from the hydrostatic approximation.
  The suitable fourth order differential operator is obtained
using the relation between $v$ and $T$.
\end{proof}


\subsection{Vertical-geostrophic equations} %subsec4.2
    We recall here the equations studied in \cite{5} in domains with
vanishing depth.
These equations may also be obtained in lubrication and thin film
theory.  We consider that the flow satisfies the vertical-geostrophic
equations
\begin{equation}
\begin{gathered}
 - \partial_z^2 v + k v^\bot + \nabla_{\rm x} p = 0, \quad \partial_z p = 0, \\
 \mathop{\rm div}_x v + \partial_z w = 0,  \\
 \partial_z v = g,  \quad w = 0 \text{ on } \Gamma_s, \\
 (v,w) = 0 \text{ on } \Gamma_b, \quad
  \overline v \cdot n_{\partial \mathcal{O}} = 0.
\end{gathered}\label{e31}
\end{equation}
  If $k= 0$, integrating two times with respect to $z$, we prove that
$(v,w)$ is given by
\begin{equation}
v = {1 \over 2} (z^2-h^2)\nabla_{\rm x} p + (z+h) g, \quad
    w = - \int_0^z {\rm div}_{\rm x} v. \label{e32}
\end{equation}
Similar formula may be obtained if $k\not = 0$,  see for instance
\cite{5}.
  Using the same lines than in Section \ref{sec3}, there exists
$\Psi$ such that
$$
\nabla_{\rm x}^\bot \Psi = {h^3\over 3} \nabla_{\rm x} p + {h^2\over 2}
g,
$$
and therefore we get, from \eqref{e32},
$$
v = {3 \over 2\,h^3} (z^2-h^2) \nabla_{\rm x}^\bot \Psi + (z+h) g
       - {3\over 4 h}(z^2-h^2)g, \quad
    w = - \int_0^z {\rm div}_{\rm x} v.
$$
  Using the regularity of the stream function $\Psi$ associated to
$\overline v$,
that implies
$$
h^{3/2}\nabla_{\rm x} ({1\over h^3}\nabla_{\rm x} \Psi) \in (L^2(\mathcal
{O}))^4.
$$
  This gives a regularity result of $(v,w)$ in
$(H^1(\Omega))^3$ assuming that $g$ is smooth enough, $h$ satisfying
the hypothesis of Theorem \ref{thm2} and
$h\in W^{2,\infty}$.
   We note that ${\rm div}_{\rm x} v$ depends only on the first
derivative of $\Psi$.

\subsection{Splitting-projection methods for the hydrostatic
Navier-Stokes equations} %subsec4.3

   In many geophysical fluids it is natural to assume hydrostatic
pressure and the ``rigid lid'' hypothesis (see \cite{15}), hence the
3D Navier-Stokes equations gives to the hydrostatic Navier-Stokes
equations, see \cite{12,13}.
   For instance, these equations model the general
circulation in lakes and oceans.
   For simplicity, we impose constant density, cartesian coordinates and
temperature and salinity effects decoupled of the dynamic flow. This
gives the hydrostatic Navier-Stokes equations:
\begin{equation}
\begin{gathered}
  \partial _t v \,+\, ( u\cdot \nabla ) v\,
- \nu \, \Delta v + f v^\bot\,+\,
 \nabla_{\rm x}\, p_{s} =\, {\rm F}  \text{ in $ \Omega\times (0,T)$,}
  \\
 w (t;{\rm x},z)  = \int _z^0 \mathop{\rm div}{}_{\rm x} v
  (t;{\rm x},s)\, ds, \quad
  \mathop{\rm div}{}_{\rm x} \overline v  = 0  \text{ in $ \mathcal{O}
\times
(0,T)$,} \\
 v|_{ \Gamma_b} = 0, \quad  \nu \partial_z v|_{ \Gamma_s}    = {\rm g
}_s ,\quad
  v_{\vert t= 0}  = v_0   \text{ in } \Omega.
\end{gathered}\label{e33}
\end{equation}
   The unknowns of this problem are:
$u= (v,w):\Omega\times  (0,T) \to \mathbb{R}^3$ the 3D velocity field
(with  $v= (v_1 , v_2 )$ the corresponding horizontal velocity field) and
$p_s: \mathcal{O}\times (0,T) \to \mathbb{R}$ a potential function, defined only
on the surface $\mathcal{O}$.
  The term $f v^\bot$ represents the Coriolis forces, being
$f $ a constant (depending on the  angular velocity of the earth and  the
latitude), ${\rm F} : \Omega \times (0,T) \to \mathbb{R} ^2$ are the
horizontal external forces, ${\rm g}_s:\Gamma_s \times (0,T) \to \mathbb{R} ^2$
models the friction effects on the surface and $v_0: \Omega \to \mathbb{R}^2$ is the initial velocity.

\subsection*{(i) Description of the splitting projection scheme,\cite{9,10}}

    We consider a regular partition of the time interval $[0,T]$ by
$M$ subintervals of length $k=T/M$, hence we have the nodes $\{ t_m= m\,
k\}_{m= 0,\dots,M}$.

\noindent {\bf Initialization:} Given $v^0 \simeq v_0 $, to
compute  $w^0=\int_z^0 {\rm div}_{\rm x} v^0 $  in $\Omega$.

\noindent {\bf Time Step $m$:}

{\bf Sub-step 1:}
  Given $u^{m-1}=(v^{m-1}, w^{m-1}) $, to find
$\widetilde {v^{m}}:\Omega \to\mathbb{R}^2$ such that
\begin{equation}
\begin{gathered}
  \frac 1k ({v^{m}}- v^{m-1})+ (u^{m-1} \cdot
\nabla )\widetilde {v^{m}}
  -\nu\Delta \widetilde {v^{m}} + f(\widetilde {v^{m}})^\bot =
{\rm f}^{m}\quad  \text{ in  } \Omega ,  \\
  \nu  \partial_z \widetilde {v^{m}}|_{\Gamma_s}={\rm g}_s^{m}, \quad
\widetilde {v^{m}}|_{\Gamma_b}  = 0 .  \end{gathered}\label{e34}
\end{equation}

{\bf Sub-step 2}: Given $\widetilde {v^{m}}$, to find
  $v^{m}:\Omega \to\mathbb{R}^2 $ and $p_s^{m}:\mathcal{O} \to\mathbb{R}$
 such that
\begin{equation}
\begin{gathered}
 \frac 1k (v^{m}-{\widetilde v}^{m})
  + \nabla_{\rm x}\, p_s^{m}  =  0 \quad \text{ in } \Omega , \\
 \mathop{\rm div}_{\rm x} \overline {v^{m}} =0\quad
\text{ in } \mathcal{O},\quad
\overline {v^{m}} \cdot  n_{\partial \mathcal{O}  } = 0
  \text{ on } \partial \mathcal{O}.
\end{gathered}\label{e35}
\end{equation}

  {\bf Sub-step 3}: Given $v^{m}$,
to compute $w^{m} ({\rm x},z)  =  \int _z^0 {\rm div}_{\rm
x} v^{m}({\rm x},s) \, ds$. \medskip

\noindent   Note that sub-step 2 is equivalent to an elliptic-Neuman prob=
lem in
$\mathcal{O}$ for $p_s^m$ (which degenerates when $h= 0$ on $\partial\mathcal{O}$)
and the explicit relation $v_m = \widetilde v_m - k \nabla_x p_s^m$.
   Indeed, integrating with respect top $z$ from $-h$ to $0$, we get
\begin{equation}
\overline{v^m} - \overline{{\widetilde v}^m}
       + k h \nabla_x p_s^m = 0 \text{ in } \mathcal{O}.
      \label{e36}
\end{equation}
  Hence, taking ${\rm div}_x$ and multiplying by $n_{\partial\mathcal{O}}$,
$$
k\mathop{\rm div}_x(h\nabla_x p_s^m) = \mathop{\rm div}_x \overline{{\widetilde v}^m}
  \text{ in } \mathcal{O}, \quad
     h \nabla_x p_s^m\cdot n_{\partial\mathcal{O}} = 0.
$$

\subsection*{(ii) A Coriolis correction scheme}
  In \cite{9}, a variant of this scheme is considered, where one has
Coriolis in an explicit way in the first sub-step (changing
$({\widetilde v}^{m})^\bot$ by
$(v^{m-1})^\bot$ in \eqref{e34}, but with a correction of the Coriolis
terms in the second sub-step, changing \eqref{e35} by
\begin{equation}
\begin{gathered}
 \frac 1k (v^{m}-{\widetilde{v}}^{m})
+ f(v^{m} - v^{m-1})^\bot + \nabla_{\rm x}\, p_s^{m}  =  0 \quad \text{in }
\Omega , \\
 \nabla_{\rm x} \cdot\overline{v^{m}} =0 \quad
\text{in } \mathcal{O},\quad
  \overline{v^{m}} \cdot  n_{\partial \mathcal{O}  } = 0
\quad \text{on  }\partial \mathcal{O}.
\end{gathered}\label{e37}
\end{equation}

\subsection*{Stability of these schemes}
  It is well known in the Navier-Stokes framework (see for instance
\cite{16}, \cite{17}), that splitting-projection schemes are
stables in the
$H^1$-norm for both velocities $\widetilde v^{m}$ and $v^m$.
  In order to obtain stability in $H^1(\Omega)$ for $v^m$ in the
hydrotatic Navier-Stokes equations (that is  fundamental to prove the
convergence of the scheme) it is necessary to get the $H^2(\Omega)$
stability for the pressure
$p_s^m$.
   Now, in the Primitive Equations framework,  taking into account that
$p_s^m$ does not depend on $z$, this $H^2(\Omega)$ regularity is equivale=
nt
to a weight $H^2(\mathcal{O})$ regularity depending on the bathymetry $h$.
   This $H^2$ weighted regularity can be deduced from Theorem
\ref{thm2}, see \cite{9} following the steps described in
Section \ref{sec3}.

\subsection{Great-Lake equations} %subsec4.4
   This example is not an illustration of Section \ref{sec3} but, as we
will see,  we obtain the same degenerate equation on the stream function.
   Let us give here an equation studied in \cite{11} in a domain
with a depth with sidewalls.
   This system will be studied in \cite{4} with a depth vanishing on
the shore.

  We assume that the mean flow $(u,p)$ satisfies in a
simply connected two-dimen\-sional space domain $\mathcal{O}$, the following system
\begin{gather*}
 \partial_t u + u\cdot \nabla_{\rm x} u + \nabla_{\rm x} p =0 \text{ in }
\mathcal{O}, \\
 \nabla_{\rm x} \cdot (hu) = 0 \text{ in } \mathcal{O}, \\
 hu\cdot n = 0 \text{ on } \partial\mathcal{O}, \\
 u(0)= u_0.
\end{gather*}
  Denoting $\omega = \mathop{\rm curl}_{\rm x} u / h$, and using the equation
$\nabla_{\rm x}\cdot (hu) = 0$, we can easily prove that the system is equivalent
to
\begin{equation}
\begin{gathered}
 \partial_t \omega + u\cdot \nabla_{\rm x} \omega  = 0, \quad
    u = {\nabla_{\rm x}^\bot\Psi}/h \text{ in } \mathcal{O}, \\
 \nabla_{\rm x}\cdot (\frac 1h \nabla_{\rm x}\Psi) = h\omega
  \text{ in } \mathcal{O}, \quad
   \Psi = 0 \text{ on } \partial\mathcal{O}, \\
 u(0)= u_0.
\end{gathered}\label{e38}
\end{equation}
  We obtain again the degenerate equation on $\Psi$ that we have
studied in Section \ref{sec2}.
  Now, the "$H^2$ regularity" obtained in Section \ref{sec2} is not
enough to follow the study  of Youdovitch made on the standard Euler
equations. An "$L^r$ regularity" in weighted spaces is required.
   This result remains as an interesting open question.


\subsection*{Acknowledgements}
   The first author has been partially supported by the IDOPT project in
Grenoble.
   The third author has been partially supported by project
 BFM2000-1317.

\begin{thebibliography}{00}

\bibitem{1} M. Amara, D. Capatina.  Title ???,
In preparation


\bibitem{2}
P. Bolley, J. Camus.
Sur une classe d'op\'erateurs elliptiques et d\'eg\'en\'er\'es \`a
plusieurs variables. Contributions \`a l'analyse fonctionnelle, pp.
55--140.
\textit{Bull. Soc. Math. France},  Mem. No. 34, Soc. Math. France,
Paris, (1973).

\bibitem{3}
P. Bolley, J. Camus, G. M\'etivier.
Estimations de Schauder et r\'egularit\'e hold\'e\-rienne pour une
classe de probl\`emes aux limites singuliers.
\textit{Comm. Partial Differential Equations},  {\bf 11}, (1986), no.
11, 1135--1203.

\bibitem{4}
D. Bresch, A. Kazhikhov, M. Sy.
On the rigid-lid shallow-water equation with vanishing depth on the shore.
In preparation.

\bibitem{5}
D. Bresch, J. Lemoine, J. Simon.
 A geostrophic model with vertical diffusion.
\textit{Nonlinear Analysis, Theory, Methods and Applications},
 {\bf 43}/4, (2001), {\rm 449--470}.

\bibitem{6}
D. Bresch, M. Sy.
Porous convection in rotating thin domains: the planetary
geostrophic equations, used in geophysical fluid dynamics,
revisited,
\it Cont. Mech. Thermodyn. \rm 15 (2003) 3, 247-263.


\bibitem{7}
A. Kufner.
\textit{Weighted Sobolev Spaces}.
 Teubner-Texte zur Mathematik, Band 31, (1980).

\bibitem{8}
H. Br\'ezis.
\textit{Analyse fonctionnelle, th\'eorie et applications},
 Masson (1987).

\bibitem{9}
F. Guill\'en-Gonz\'alez, M. V. Redondo-Neble.
Numerical Analysis for some  time fractional-step projection
methods for the Pimitive Equations.
 In preparation.

\bibitem{10}
F. Guill\'en-Gonz\'alez, M. V. Redondo-Neble, J. R. Rodr\'{\i}-guez-Galv\'en.
An\'alisis Num\'erico y resoluci\'on efectiva de las Ecuaciones
Primitivas con esquemas de tipo proyecci\'on.
In \textit{Proc. XVII CEDYA, VII CMA}, Universidad de Salamanca (2001).

\bibitem{11}
C. D. Levermore, M. Oliver, E. S. Titi.
 Global well-posedness for models of shallow water in a basin with
a varying bottom.
\textit{Indiana Univ. Math.}  J. 45 (1996), no. 2, 479--510.

\bibitem{12}
J.--L. Lions, R. Teman, S. Wang.
 New formulations of the primitives equations of the  atmosphere
and applications.
\textit{Nonlinearity}, 5 (1992), 237-288.

\bibitem{13}
J.--L. Lions, R. Teman, S. Wang.
On the equations of the large scale Ocean.
\textit{Nonlinearity},  5 (1992), 1007-1053.

\bibitem{14}
J. Necas
\textit{Les m\'ethodes directes en th\'eorie des \'equations elliptiques}.
\rm Masson, (1967).

\bibitem{15}
J. Pedlosky.
\textit{Geophysical fluid dynamics}.  Springer-Verlag, (1987).

\bibitem{16}
J. Shen.
On error estimates of projection methods for Navier-Stokes equations:
first-order schemes.
\textit{SIAM Journal Num. Anal.}  Vol. 29 (1992),  57-77.

\bibitem{17}
R. Temam.
\textit{Navier-Stokes Equations}.
AMS Chelsea publishing, (2001).

\bibitem{18}
H.-M. Yin.
Optimal regularity of solution to a degenerate elliptic system arising in
electromagnetic fields.
\textit{Comm. Pure and Appl. Anal.}
 {\bf 1}, 1, (2002), 127--134.
\end{thebibliography}

\end{document}

