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	% $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