s_algorithm/text/nonlin_frsurf.tex

% $Header:  $
% $Name:  $

%\section{Time Integration}

\subsection{Non-linear free surface}

Recently, 2 options have added to the model
(and therefore, have not yet been extensively tested)
that concern the free surface formulation.

%------------------------------------------
\subsubsection{Non-uniform linear-relation for the surface potential}

The linear relation between 
surface pressure / geo- potential ($\Phi_{surf}$)
and surface displacement ($\eta$)
has been considered as uniform ($b_{surf} = Constant$)
but is in fact
dependent on the position 
since we have 
$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_{surf} \eta$

For the ocean, $b_{surf} = g \rho_{surf} / \rho_o \simeq g$
is a fairly good approximation since the relative difference
in surface density are usually small.
For the atmosphere, because of topographic effects,
the reference surface pressure $R_o$ has large spacial differences
that are responsible for significant $b_surf$ variations
(from 0.8 to 1.2 $[]$).

Practically, in an oceanic configuration or
when the default value (TRUE) of the parameter {\bf uniformLin\_PhiSurf}
is used ) 
then $b_{surf}$ is simply set to $g$ for the ocean
and $1.$ for the atmosphere.\\
Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to
evaluate $b_{surf}$ from the reference Theta,Q vertical profile
({\it S/R INI\_LINEAR\_PHISURF})
according to the reference surface pressure $P_o$ ($=R_o$):
$b_{surf} = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$

%------------------------------------------
\subsubsection{Free surface effect on column total thickness
(Non-linear free surface)}

The total thickness of the fluid column is
$r_{surf} - R_{min} = \eta + R_o - R_{min}$
In the linear free surface approximation
(detailed before), only the fixed part of
it ($R_o - R_{min})$ is considered when we integrate the 
continuity equation or compute tracer and momentum advection term.

This approximation is dropped when using 
the non linear free surface formulation. 
Details follow here after for the barotropic part
and in the 2 next subsection for the baroclinic
part.

%------------------------------------------
% Non Linear Barotropic part

The continuous form of the model equations remains 
unchanged, except for the 2D continuity equation
(\ref{eq-cont-2D}) that is now integrated 
from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ :

\begin{displaymath}
\epsilon_{fs} \partial_t \eta =
\left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
- {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v} dr
+ \epsilon_{fw} (P-E)
\end{displaymath}

Since $\eta$ has a direct effect on the horizontal
velocity (through $\nabla_r \Phi_{surf}$), this
add a non linear term to the free surface equation.

Regarding the time discretization of this non linear part,
several option are now being tested:

With the column thickness evaluated at time step $n$,
and the surface potential gradient still implicit,
equation (\ref{solve_2d} \& \ref{solve_2d_rhs})
become:
\begin{eqnarray*}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n}+R_o-R_{min})
{\bf \nabla}_r B_o {\eta}^{n+1}
= {\eta}^*
%\label{solve_2d}
\end{eqnarray*}
where
\begin{eqnarray*}
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
\Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v}^* dr
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
%\label{solve_2d_rhs}
\end{eqnarray*} 
This method requires to update the solver matrix at each time step.

Alternatively, the non-linear contribution can be evaluated fully
explicitly:
\begin{eqnarray*}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min})
{\bf \nabla}_r B_o {\eta}^{n+1}
= {\eta}^*
+{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n})
{\bf \nabla}_r B_o {\eta}^{n}
\end{eqnarray*} 
This formulation allows to keep the initial solver matrix
since the non-linear free surface only affects the RHS.

An other option is a "linearized" formulation where the 
total column thickness appears only in the integral term of 
the RHS (\ref{solve_2d_rhs}) but not directly in 
the equation (\ref{solve_2d}).

%------------------------------------------
\subsubsection{Free surface effect on the surface level thickness
(Non-linear free surface): Tracer advection}

To ensure a global tracer conservation (i.e., the total amount)
as well as the local one (see tracer section for more details), 
the change in the surface level thickness must be consistent with
the way the continuity equation is integrated, both in 
in the barotropic part (to find $\eta$) and baroclinic part
(to find $\dot{r}$).

To illustrate this, let's consider the shallow water model,
with uniform cartesian horizontal grid:
$$
\partial_t h + \nabla \cdot h \vec{\bf v} = 0
$$
where $h$ is the total thickness of the water column.
To conserve the tracer $\theta$ we have to discretize:
$$
\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
$$
Using the implicit (non-linear) free surface described before, we have:
\begin{eqnarray*}
h^{n+1} = h^{n} - \Delta_t \nabla (h^n \, \vec{\bf v}^{n+1} ) \\
\end{eqnarray*}
The discretized form of the tracer equation must use the same
"geometry" to compute the tracer fluxes, that is, the same value of
h, as the continuity equation did:
\begin{eqnarray*}
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n 
        - \Delta_t \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1})
\end{eqnarray*}

In order to deal with the Adams-Bashforth time stepping,
we implement this scheme slightly differently:
the variation of the water column thickness (from
$h^n$ to $h^{n+1}$)
is not incorporated directly to the tracer equation.
Instead,
the model continues to evaluate the $G_\theta$ term 
exactly as it did with the Linear free surface formulation,
with the "{\it surface correction}" turned "on" (see tracer section):
$$
G_\theta = \left(- \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1}) 
         + w_{surf}^{n+1} \theta^n \right) / h^n
$$
Then after,
thickness variation (expansion/reduction) is taken into account :
$$
\theta^{n+1} = (\theta^n + \Delta_t G_\theta^{(n+1)}) \frac{h^{n+1}}{h^n}
$$
Note that with a simple forward time step (no Adams-Bashforth), 
since
$
w_{surf} = - \nabla (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
/ \Delta_t
$
those 2 formulations are equivalent. 

The implementation in the MITgcm follows this scheme.
The model "geometry" (here the {\bf hFacC,W,S}) is updated
at the beginning of {\it SOLVE\_FOR\_PRESSURE}
according to the current $\eta$ field.
Then, at the end of the time step, the variables are
advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$.
At the next time step, the tracer tendency ($G_\theta$) is computed 
using the same geometry, that is now consistent with
$\eta^{n-1}$.
Finally, in S/R {\it TIMESTEP\_TRACER}, the expansion/reduction
ratio is applied to the surface level to compute the new tracer field.

%------------------------------------------
\subsubsection{Free surface effect on the surface level thickness
(Non-linear free surface): Momentum advection}     

Regarding momentum advection,
the vector invariant formulation is similar to the
advective form (as opposed to the flux form) and therefore
does not need specific modification to include the 
free surface effect on the surface level thickness.
Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s) 
at one given place (like describe before) is sufficient.

With the flux form formulation, advection of momentum
can be treated exactly as the tracer advection is,
Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
subroutine {\it TIMESTEP}.

1	jmc	1.1	% $Header: $
2			% $Name: $
3
4			%\section{Time Integration}
5
6			\subsection{Non-linear free surface}
7
8			Recently, 2 options have added to the model
9			(and therefore, have not yet been extensively tested)
10			that concern the free surface formulation.
11
12			%------------------------------------------
13			\subsubsection{Non-uniform linear-relation for the surface potential}
14
15			The linear relation between
16			surface pressure / geo- potential ($\Phi_{surf}$)
17			and surface displacement ($\eta$)
18			has been considered as uniform ($b_{surf} = Constant$)
19			but is in fact
20			dependent on the position
21			since we have
22			$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_{surf} \eta$
23
24			For the ocean, $b_{surf} = g \rho_{surf} / \rho_o \simeq g$
25			is a fairly good approximation since the relative difference
26			in surface density are usually small.
27			For the atmosphere, because of topographic effects,
28			the reference surface pressure $R_o$ has large spacial differences
29			that are responsible for significant $b_surf$ variations
30			(from 0.8 to 1.2 $[]$).
31
32			Practically, in an oceanic configuration or
33			when the default value (TRUE) of the parameter {\bf uniformLin\_PhiSurf}
34			is used )
35			then $b_{surf}$ is simply set to $g$ for the ocean
36			and $1.$ for the atmosphere.\\
37			Turning {\bf uniformLin\_PhiSurf} to "FALSE", allows to
38			evaluate $b_{surf}$ from the reference Theta,Q vertical profile
39			({\it S/R INI\_LINEAR\_PHISURF})
40			according to the reference surface pressure $P_o$ ($=R_o$):
41			$b_{surf} = c_p \kappa (P_o / Pc)^{(\kappa - 1)} \theta_{ref}(P_o)$
42
43			%------------------------------------------
44			\subsubsection{Free surface effect on column total thickness
45			(Non-linear free surface)}
46
47			The total thickness of the fluid column is
48			$r_{surf} - R_{min} = \eta + R_o - R_{min}$
49			In the linear free surface approximation
50			(detailed before), only the fixed part of
51			it ($R_o - R_{min})$ is considered when we integrate the
52			continuity equation or compute tracer and momentum advection term.
53
54			This approximation is dropped when using
55			the non linear free surface formulation.
56			Details follow here after for the barotropic part
57			and in the 2 next subsection for the baroclinic
58			part.
59
60			%------------------------------------------
61			% Non Linear Barotropic part
62
63			The continuous form of the model equations remains
64			unchanged, except for the 2D continuity equation
65			(\ref{eq-cont-2D}) that is now integrated
66			from $R_{min}(x,y)$ up to $r_{surf}=R_o+\eta$ :
67
68			\begin{displaymath}
69			\epsilon_{fs} \partial_t \eta =
70			\left. \dot{r} \right\|_{r=r_{surf}} + \epsilon_{fw} (P-E) =
71			- {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v} dr
72			+ \epsilon_{fw} (P-E)
73			\end{displaymath}
74
75			Since $\eta$ has a direct effect on the horizontal
76			velocity (through $\nabla_r \Phi_{surf}$), this
77			add a non linear term to the free surface equation.
78
79			Regarding the time discretization of this non linear part,
80			several option are now being tested:
81
82			With the column thickness evaluated at time step $n$,
83			and the surface potential gradient still implicit,
84			equation (\ref{solve_2d} \& \ref{solve_2d_rhs})
85			become:
86			\begin{eqnarray*}
87			\epsilon_{fs} {\eta}^{n+1} -
88			{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n}+R_o-R_{min})
89			{\bf \nabla}_r B_o {\eta}^{n+1}
90			= {\eta}^*
91			%\label{solve_2d}
92			\end{eqnarray*}
93			where
94			\begin{eqnarray*}
95			{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
96			\Delta t {\bf \nabla}_r \cdot \int_{R_{min}}^{R_o+\eta} \vec{\bf v}^* dr
97			\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
98			%\label{solve_2d_rhs}
99			\end{eqnarray*}
100			This method requires to update the solver matrix at each time step.
101
102			Alternatively, the non-linear contribution can be evaluated fully
103			explicitly:
104			\begin{eqnarray*}
105			\epsilon_{fs} {\eta}^{n+1} -
106			{\bf \nabla}_r \cdot \Delta t^2 (R_o-R_{min})
107			{\bf \nabla}_r B_o {\eta}^{n+1}
108			= {\eta}^*
109			+{\bf \nabla}_r \cdot \Delta t^2 (\eta^{n})
110			{\bf \nabla}_r B_o {\eta}^{n}
111			\end{eqnarray*}
112			This formulation allows to keep the initial solver matrix
113			since the non-linear free surface only affects the RHS.
114
115			An other option is a "linearized" formulation where the
116			total column thickness appears only in the integral term of
117			the RHS (\ref{solve_2d_rhs}) but not directly in
118			the equation (\ref{solve_2d}).
119
120			%------------------------------------------
121			\subsubsection{Free surface effect on the surface level thickness
122			(Non-linear free surface): Tracer advection}
123
124			To ensure a global tracer conservation (i.e., the total amount)
125			as well as the local one (see tracer section for more details),
126			the change in the surface level thickness must be consistent with
127			the way the continuity equation is integrated, both in
128			in the barotropic part (to find $\eta$) and baroclinic part
129			(to find $\dot{r}$).
130
131			To illustrate this, let's consider the shallow water model,
132			with uniform cartesian horizontal grid:
133			$$
134			\partial_t h + \nabla \cdot h \vec{\bf v} = 0
135			$$
136			where $h$ is the total thickness of the water column.
137			To conserve the tracer $\theta$ we have to discretize:
138			$$
139			\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v})= 0
140			$$
141			Using the implicit (non-linear) free surface described before, we have:
142			\begin{eqnarray*}
143			h^{n+1} = h^{n} - \Delta_t \nabla (h^n \, \vec{\bf v}^{n+1} ) \\
144			\end{eqnarray*}
145			The discretized form of the tracer equation must use the same
146			"geometry" to compute the tracer fluxes, that is, the same value of
147			h, as the continuity equation did:
148			\begin{eqnarray*}
149			h^{n+1} \, \theta^{n+1} = h^n \, \theta^n
150			- \Delta_t \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1})
151			\end{eqnarray*}
152
153			In order to deal with the Adams-Bashforth time stepping,
154			we implement this scheme slightly differently:
155			the variation of the water column thickness (from
156			$h^n$ to $h^{n+1}$)
157			is not incorporated directly to the tracer equation.
158			Instead,
159			the model continues to evaluate the $G_\theta$ term
160			exactly as it did with the Linear free surface formulation,
161			with the "{\it surface correction}" turned "on" (see tracer section):
162			$$
163			G_\theta = \left(- \nabla (h^n \, \theta^n \, \vec{\bf v}^{n+1})
164			+ w_{surf}^{n+1} \theta^n \right) / h^n
165			$$
166			Then after,
167			thickness variation (expansion/reduction) is taken into account :
168			$$
169			\theta^{n+1} = (\theta^n + \Delta_t G_\theta^{(n+1)}) \frac{h^{n+1}}{h^n}
170			$$
171			Note that with a simple forward time step (no Adams-Bashforth),
172			since
173			$
174			w_{surf} = - \nabla (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
175			/ \Delta_t
176			$
177			those 2 formulations are equivalent.
178
179			The implementation in the MITgcm follows this scheme.
180			The model "geometry" (here the {\bf hFacC,W,S}) is updated
181			at the beginning of {\it SOLVE\_FOR\_PRESSURE}
182			according to the current $\eta$ field.
183			Then, at the end of the time step, the variables are
184			advanced in time, so that $\eta^n$ becomes $\eta^{n-1}$.
185			At the next time step, the tracer tendency ($G_\theta$) is computed
186			using the same geometry, that is now consistent with
187			$\eta^{n-1}$.
188			Finally, in S/R {\it TIMESTEP\_TRACER}, the expansion/reduction
189			ratio is applied to the surface level to compute the new tracer field.
190
191			%------------------------------------------
192			\subsubsection{Free surface effect on the surface level thickness
193			(Non-linear free surface): Momentum advection}
194
195			Regarding momentum advection,
196			the vector invariant formulation is similar to the
197			advective form (as opposed to the flux form) and therefore
198			does not need specific modification to include the
199			free surface effect on the surface level thickness.
200			Updating the {\bf hFacC,W,S} and the {\bf recip\_hFac}(s)
201			at one given place (like describe before) is sufficient.
202
203			With the flux form formulation, advection of momentum
204			can be treated exactly as the tracer advection is,
205			Here the expansion/reduction factors ($hFacW^{n+1}/hFacW^n$ for $u$
206			and $hFacS^{n+1}/hFacS^n$ for $v$) are simply applied in the
207			subroutine {\it TIMESTEP}.
208