\chapter{Incompressibility}
\label{app:incompressibility}


The numerical procedure for the incompressible correction is summarized  
in the following:
\bn

\item Guess the trial surface pressure $\tilde{p}_{s}^{n}$ from the 
previous time step ($\tilde{p}_{s}^{n} =p_{s}^{n-1}$) and integrate the 
hydrostatic equation~(\ref{eq:Hydrostatic}) to get the intermediate 
pressure field $\tilde{p}^{n}$ over the whole domain.

\item Update the trial integral average $\tilde{u}^{n+1}$,$\tilde{v}^{n+1}$
in the control volume using equation~(\ref{eq:numerical_discritization}) 
and the analogous control volume equation for $v$.
   
\item Interpolate $\tilde{u}^{n+1},\tilde{v}^{n+1}$ to 
$\tilde{U}_{i,j,k}^{n+1},\tilde{V}_{i,j,k}^{n+1}$
at the cell face using $4^{th}$ order interpolation from
equation~(\ref{eq:interp}).

\item Adjust the rigid-lid pressure such that the divergence
of the vertically integrated horizontal velocity matches the specified 
value(near zero, but small precipitation/evaporation and river sources are
included to give non-zero $W$ at the rigid lid which is slightly porous
to fresh water).
This procedure uses mass conservation and the detailed procedure is 
listed as follows.
\bi
\item Set the final pressure $p^{n}=\tilde{p}^{n}+\Delta \tilde{p}$
where $\Delta \tilde{p}$ is due to the change of rigid-lid pressure
and thus independent of depth. $\tilde{p}^{n}$ is derived  
hydrostatically from trial rigid-lid pressure in the first step. Then 
the final velocities can be written as
\begin{align}
U^{n+1}&=\tilde{U}^{n+1}+\Delta \tilde{U} \label{eq:U_final} \\
V^{n+1}&=\tilde{V}^{n+1}+\Delta \tilde{V},
\label{eq:V_final}
\end{align}
where
$\Delta \tilde{U}=-\Delta t \fc{\pd \Delta \tilde{p}}{\pd x}$ and
$\Delta \tilde{V}=-\Delta t \fc{\pd \Delta \tilde{p}}{\pd y}$.
\item The final velocity field satisfies mass conservation from equation
(\ref{eq:3Dmass}) so that the vertically-integrated horizontal velocity   
can be determined by vertical velocity at the top and bottom:
\be
\int_{0}^{D}(\fc{\pd U^{n+1}}{\pd x}+\fc{\pd V^{n+1}}{\pd y})dz=   
W^{n+1}(0)-W^{n+1}(D)=0
\ee

\item Substitute the final velocity integration in terms of trial 
velocity and
the correction in equations (\ref{eq:U_final}) and (\ref{eq:V_final}),
thus deriving a 2-D Poisson equation for rigid-lid pressure adjustment. 
Solve the Poisson equation (\ref{eq:2d_poisson}) using an efficient EVP 
elliptic solver~\cite{Roache95}.
\be
\int_{0}^{D}(\fc{\pd \tilde{U}^{n+1}}{\pd x}+\fc{\pd 
\tilde{V}^{n+1}}{\pd y})dz=
\int_{0}^{D}-(\fc{\pd \Delta \tilde{U}}{\pd x}+\fc{\pd \Delta 
\tilde{V}}{\pd y})dz
\label{eq:2d_poisson}
\ee
\ei

\item The velocity $U^{n+1},V^{n+1}$ is corrected from the pressure
adjustment $\Delta \tilde{p}$ in the last step.

\item Finally, interpolate cell average changes
$\Delta \tilde{u},\Delta \tilde{v}$
from $\Delta \tilde{U}$,$\Delta \tilde{V}$ using $4^{th}$ order
interpolation, then  
$u^{n+1}=\tilde{u}^{n+1}+\Delta \tilde{u}$ and
$v^{n+1}=\tilde{v}^{n+1}+\Delta \tilde{v}$.
\en

\newpage
