s_algorithm/text/nonlin_frsurf.tex

% $Header: /u/gcmpack/mitgcmdoc/part2/nonlin_frsurf.tex,v 1.3 2001/09/24 19:30:40 jmc Exp $
% $Name:  $


\subsection{Non-linear free surface}

Recently, two options have been added to the model (and 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$) could be considered
to be a constant ($b_s=$ constant)
\marginpar{add a reference to part.1 here}
but is in fact dependent on the position ($x,y,r$)
since we linearize:
$$\Phi_{surf}=\int_{R_o}^{R_o+\eta} b dr \simeq b_s \eta
~\mathrm{with}~ b_s = b(\theta,S,r)_{r=R_o} 
\simeq b_s(\theta_{ref}(R_o),S_{ref}(R_o),R_o)$$
Note that, for convenience, the effect on $b_s$ of the local surface
$\theta,S$ are not considered here, but are incorporated in to
$\Phi'_{hyd}$.

For the ocean, $b_s = g \rho_{surf} / \rho_o \simeq g$ is a very good
approximation since the relative difference in surface density are
usually small and only due to local $\theta,S$ gradients (because the
upper surface, $R_o = 0$, is essentially flat). Therefore, they can
easily be incorporated in $\Phi'_{hyd}$.

For the atmosphere, however, because of topographic effects, the
reference surface pressure $R_o$ has large spacial variations that
are responsible for significant $b_s$ variations (from 0.8 to 1.2
$[m^3/kg]$). For this reason, we use a non-uniform linear coefficient
$b_s$.

In practice, in an oceanic configuration or when the default value
(TRUE) of the parameter {\bf uniformLin\_PhiSurf} is used, then $b_s$
is simply set to $g$ for the ocean and $1.$ for the atmosphere.
Turning {\bf uniformLin\_PhiSurf} to "FALSE", tells the code to
evaluate $b_s$ from the reference vertical profile $\theta_{ref}$
({\it S/R INI\_LINEAR\_PHISURF}) according to the reference surface
pressure $P_o$ ($=R_o$): $b_s = 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_{fixed} =
\eta + R_o - R_{fixed}$ In the linear free surface approximation
(detailed before), only the fixed part of it ($R_o - R_{fixed})$ 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.  Here we discuss sections the barotropic part. In
sections \ref{sect:freesurf-tracer-advection} and
\ref{sect:freesurf-momentum-advection} we consider the baroclinic
component.


The continuous form of the model equations remains unchanged, except
for the 2D continuity equation (\ref{eq-tCsC-eta}) which is now
integrated from $R_{fixed}(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}_h \cdot \int_{R_{fixed}}^{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_h \Phi_{surf}$), this adds a non-linear term to the free
surface equation. Several options for the time discretization of this
non-linear part have been tested.

If the column thickness is evaluated at time step $n$, and with
implicit treatment of the surface potential gradient, equations
(\ref{eq-solve2D} and \ref{eq-solve2D_rhs}) becomes:
\begin{eqnarray*}
\epsilon_{fs} {\eta}^{n+1} -
{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed})
{\bf \nabla}_h b_s {\eta}^{n+1}
= {\eta}^*
\end{eqnarray*}
where
\begin{eqnarray*}
{\eta}^* = \epsilon_{fs} \: {\eta}^{n} -
\Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr
\: + \: \epsilon_{fw} \Delta_t (P-E)^{n}
\end{eqnarray*} 
This method requires us 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}_h \cdot \Delta t^2 (R_o-R_{fixed})
{\bf \nabla}_h b_s {\eta}^{n+1}
= {\eta}^*
+{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n})
{\bf \nabla}_h b_s {\eta}^{n}
\end{eqnarray*} 
This formulation allows one to keep the initial solver matrix
unchanged though throughout the integration, since the non-linear free
surface only affects the RHS.

Finally, another option is a "linearized" formulation where the total
column thickness appears only in the integral term of the RHS
(\ref{eq-solve2D_rhs}) but not directly in the equation
(\ref{eq-solve2D}).


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

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

To illustrate this, 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 above (section
\ref{sect:??}, we have:
\begin{eqnarray*}
h^{n+1} = h^{n} - \Delta_t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) \\
\end{eqnarray*}
The discretized form of the tracer equation must adopt the same
``form'' in the computation of tracer fluxes, that is, the same value
of $h$, as used in the continuity equation:
\begin{eqnarray*}
h^{n+1} \, \theta^{n+1} = h^n \, \theta^n 
        - \Delta_t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1})
\end{eqnarray*}

For Adams-Bashforth time-stepping, we implement this scheme slightly
differently from the linear free-surface method, using two steps: the
variation of the water column thickness (from $h^n$ to $h^{n+1}$) is
not incorporated directly into the tracer equation.  Instead, the
model uses the $G_\theta$ terms (first step) as in the linear free
surface formulation (with the "{\it surface correction}" turned "on",
see tracer section):
$$
G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) 
         - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n
$$
Then, in a second step, the thickness variation (expansion/reduction)
is taken into account:
$$
\theta^{n+1} = \theta^n + \Delta_t \frac{h^n}{h^{n+1}} G_\theta^{(n+1/2)} 
$$
Note that with a simple forward time step (no Adams-Bashforth), 
since
$
\dot{r}_{surf}^{n+1} 
= - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = (h^{n+1} - h^{n})
/ \Delta_t
$
these two formulations are equivalent. 

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