New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
chap_conservation.tex in NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO/subfiles – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO/subfiles/chap_conservation.tex @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

File size: 17.3 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5\chapter{Invariants of the Primitive Equations}
6\label{chap:CONS}
7
8\thispagestyle{plain}
9
10\chaptertoc
11
12\paragraph{Changes record} ~\\
13
14{\footnotesize
15  \begin{tabularx}{\textwidth}{l||X|X}
16    Release & Author(s) & Modifications \\
17    \hline
18    {\em   4.0} & {\em ...} & {\em ...} \\
19    {\em   3.6} & {\em ...} & {\em ...} \\
20    {\em   3.4} & {\em ...} & {\em ...} \\
21    {\em <=3.4} & {\em ...} & {\em ...}
22  \end{tabularx}
23}
24
25\clearpage
26
27The continuous equations of motion have many analytic properties.
28Many quantities (total mass, energy, enstrophy, etc.) are strictly conserved in the inviscid and unforced limit,
29while ocean physics conserve the total quantities on which they act (momentum, temperature, salinity) but
30dissipate their total variance (energy, enstrophy, etc.).
31Unfortunately, the finite difference form of these equations is not guaranteed to
32retain all these important properties.
33In constructing the finite differencing schemes, we wish to ensure that
34certain integral constraints will be maintained.
35In particular, it is desirable to construct the finite difference equations so that
36horizontal kinetic energy and/or potential enstrophy of horizontally non-divergent flow,
37and variance of temperature and salinity will be conserved in the absence of dissipative effects and forcing.
38\citet{arakawa_JCP66} has first pointed out the advantage of this approach.
39He showed that if integral constraints on energy are maintained,
40the computation will be free of the troublesome "non linear" instability originally pointed out by
41\citet{phillips_TAMS59}.
42A consistent formulation of the energetic properties is also extremely important in carrying out
43long-term numerical simulations for an oceanographic model.
44Such a formulation avoids systematic errors that accumulate with time \citep{bryan_JCP97}.
45
46The general philosophy of OPA which has led to the discrete formulation presented in {\S}II.2 and II.3 is to
47choose second order non-diffusive scheme for advective terms for both dynamical and tracer equations.
48At this level of complexity, the resulting schemes are dispersive schemes.
49Therefore, they require the addition of a diffusive operator to be stable.
50The alternative is to use diffusive schemes such as upstream or flux corrected schemes.
51This last option was rejected because we prefer a complete handling of the model diffusion,
52\ie\ of the model physics rather than letting the advective scheme produces its own implicit diffusion without
53controlling the space and time structure of this implicit diffusion.
54Note that in some very specific cases as passive tracer studies, the positivity of the advective scheme is required.
55In that case, and in that case only, the advective scheme used for passive tracer is a flux correction scheme
56\citep{Marti1992?, Levy1996?, Levy1998?}.
57
58%% =================================================================================================
59\section{Conservation properties on ocean dynamics}
60\label{sec:CONS_Invariant_dyn}
61
62The non linear term of the momentum equations has been split into a vorticity term,
63a gradient of horizontal kinetic energy and a vertical advection term.
64Three schemes are available for the former (see {\S}~II.2) according to the CPP variable defined
65(default option\textbf{?}or \textbf{key{\_}vorenergy} or \textbf{key{\_}vorcombined} defined).
66They differ in their conservative properties (energy or enstrophy conserving scheme).
67The two latter terms preserve the total kinetic energy:
68the large scale kinetic energy is also preserved in practice.
69The remaining non-diffusive terms of the momentum equation
70(namely the hydrostatic and surface pressure gradient terms) also preserve the total kinetic energy and
71have no effect on the vorticity of the flow.
72
73\textbf{* relative, planetary and total vorticity term:}
74
75Let us define as either the relative, planetary and total potential vorticity, \ie, ?, and ?, respectively.
76The continuous formulation of the vorticity term satisfies following integral constraints:
77\[
78  % \label{eq:CONS_vor_vorticity}
79  \int_D {{\textbf {k}}\cdot \frac{1}{e_3 }\nabla \times \left( {\varsigma
80        \;{\mathrm {\mathbf k}}\times {\textbf {U}}_h } \right)\;dv} =0
81\]
82
83\[
84  % \label{eq:CONS_vor_enstrophy}
85  if\quad \chi =0\quad \quad \int\limits_D {\varsigma \;{\textbf{k}}\cdot
86    \frac{1}{e_3 }\nabla \times \left( {\varsigma {\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} =-\int\limits_D {\frac{1}{2}\varsigma ^2\,\chi \;dv}
87  =0
88\]
89
90\[
91  % \label{eq:CONS_vor_energy}
92  \int_D {{\textbf{U}}_h \times \left( {\varsigma \;{\textbf{k}}\times {\textbf{U}}_h } \right)\;dv} =0
93\]
94where $dv = e_1\, e_2\, e_3\, di\, dj\, dk$ is the volume element.
95(II.4.1a) means that $\varsigma $ is conserved. (II.4.1b) is obtained by an integration by part.
96It means that $\varsigma^2$ is conserved for a horizontally non-divergent flow.
97(II.4.1c) is even satisfied locally since the vorticity term is orthogonal to the horizontal velocity.
98It means that the vorticity term has no contribution to the evolution of the total kinetic energy.
99(II.4.1a) is obviously always satisfied, but (II.4.1b) and (II.4.1c) cannot be satisfied simultaneously with
100a second order scheme.
101Using the symmetry or anti-symmetry properties of the operators (Eqs II.1.10 and 11),
102it can be shown that the scheme (II.2.11) satisfies (II.4.1b) but not (II.4.1c),
103while scheme (II.2.12) satisfies (II.4.1c) but not (II.4.1b) (see appendix C).
104Note that the enstrophy conserving scheme on total vorticity has been chosen as the standard discrete form of
105the vorticity term.
106
107\textbf{* Gradient of kinetic energy / vertical advection}
108
109In continuous formulation, the gradient of horizontal kinetic energy has no contribution to the evolution of
110the vorticity as the curl of a gradient is zero.
111This property is satisfied locally with the discrete form of both the gradient and the curl operator we have made
112(property (II.1.9)~).
113Another continuous property is that the change of horizontal kinetic energy due to
114vertical advection is exactly balanced by the change of horizontal kinetic energy due to
115the horizontal gradient of horizontal kinetic energy:
116
117\begin{equation} \label{eq:CONS_keg_zad}
118\int_D {{\textbf{U}}_h \cdot \nabla _h \left( {1/2\;{\textbf{U}}_h ^2} \right)\;dv} =-\int_D {{\textbf{U}}_h \cdot \frac{w}{e_3 }\;\frac{\partial
119{\textbf{U}}_h }{\partial k}\;dv}
120\end{equation}
121
122Using the discrete form given in {\S}II.2-a and the symmetry or anti-symmetry properties of
123the mean and difference operators, \autoref{eq:CONS_keg_zad} is demonstrated in the Appendix C.
124The main point here is that satisfying \autoref{eq:CONS_keg_zad} links the choice of the discrete forms of
125the vertical advection and of the horizontal gradient of horizontal kinetic energy.
126Choosing one imposes the other.
127The discrete form of the vertical advection given in {\S}II.2-a is a direct consequence of
128formulating the horizontal kinetic energy as $1/2 \left( \overline{u^2}^i + \overline{v^2}^j \right) $ in
129the gradient term.
130
131\textbf{* hydrostatic pressure gradient term}
132
133In continuous formulation, a pressure gradient has no contribution to the evolution of the vorticity as
134the curl of a gradient is zero.
135This properties is satisfied locally with the choice of discretization we have made (property (II.1.9)~).
136In addition, when the equation of state is linear
137(\ie\ when an advective-diffusive equation for density can be derived from those of temperature and salinity)
138the change of horizontal kinetic energy due to the work of pressure forces is balanced by the change of
139potential energy due to buoyancy forces:
140
141\[
142  % \label{eq:CONS_hpg_pe}
143  \int_D {-\frac{1}{\rho_o }\left. {\nabla p^h} \right|_z \cdot {\textbf {U}}_h \;dv} \;=\;\int_D {\nabla .\left( {\rho \,{\textbf{U}}} \right)\;g\;z\;\;dv}
144\]
145
146Using the discrete form given in {\S}~II.2-a and the symmetry or anti-symmetry properties of
147the mean and difference operators, (II.4.3) is demonstrated in the Appendix C.
148The main point here is that satisfying (II.4.3) strongly constraints the discrete expression of the depth of
149$T$-points and of the term added to the pressure gradient in $s-$coordinates: the depth of a $T$-point, $z_T$,
150is defined as the sum the vertical scale factors at $w$-points starting from the surface.
151
152\textbf{* surface pressure gradient term}
153
154In continuous formulation, the surface pressure gradient has no contribution to the evolution of vorticity.
155This properties is trivially satisfied locally as (II.2.3)
156(the equation verified by $\psi$ has been derived from the discrete formulation of the momentum equations,
157vertical sum and curl).
158Nevertheless, the $\psi$-equation is solved numerically by an iterative solver (see {\S}~III.5),
159thus the property is only satisfied with the accuracy required on the solver.
160In addition, with the rigid-lid approximation, the change of horizontal kinetic energy due to the work of
161surface pressure forces is exactly zero:
162\[
163  % \label{eq:CONS_spg}
164  \int_D {-\frac{1}{\rho_o }\nabla _h } \left( {p_s } \right)\cdot {\textbf{U}}_h \;dv=0
165\]
166
167(II.4.4) is satisfied in discrete form only if
168the discrete barotropic streamfunction time evolution equation is given by (II.2.3) (see appendix C).
169This shows that (II.2.3) is the only way to compute the streamfunction,
170otherwise there is no guarantee that the surface pressure force work vanishes.
171
172%% =================================================================================================
173\section{Conservation properties on ocean thermodynamics}
174\label{sec:CONS_Invariant_tra}
175
176In continuous formulation, the advective terms of the tracer equations conserve the tracer content and
177the quadratic form of the tracer, \ie
178\[
179  % \label{eq:CONS_tra_tra2}
180  \int_D {\nabla .\left( {T\;{\textbf{U}}} \right)\;dv} =0
181  \;\text{and}
182  \int_D {T\;\nabla .\left( {T\;{\textbf{U}}} \right)\;dv} =0
183\]
184
185The numerical scheme used ({\S}II.2-b) (equations in flux form, second order centred finite differences) satisfies
186(II.4.5) (see appendix C).
187Note that in both continuous and discrete formulations, there is generally no strict conservation of mass,
188since the equation of state is non linear with respect to $T$ and $S$.
189In practice, the mass is conserved with a very good accuracy.
190
191%% =================================================================================================
192\subsection{Conservation properties on momentum physics}
193\label{subsec:CONS_Invariant_dyn_physics}
194
195\textbf{* lateral momentum diffusion term}
196
197The continuous formulation of the horizontal diffusion of momentum satisfies the following integral constraints~:
198\[
199  % \label{eq:CONS_dynldf_dyn}
200  \int\limits_D {\frac{1}{e_3 }{\mathrm {\mathbf k}}\cdot \nabla \times \left[ {\nabla
201        _h \left( {A^{lm}\;\chi } \right)-\nabla _h \times \left( {A^{lm}\;\zeta
202            \;{\mathrm {\mathbf k}}} \right)} \right]\;dv} =0
203\]
204
205\[
206  % \label{eq:CONS_dynldf_div}
207  \int\limits_D {\nabla _h \cdot \left[ {\nabla _h \left( {A^{lm}\;\chi }
208        \right)-\nabla _h \times \left( {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)}
209    \right]\;dv} =0
210\]
211
212\[
213  % \label{eq:CONS_dynldf_curl}
214  \int_D {{\mathrm {\mathbf U}}_h \cdot \left[ {\nabla _h \left( {A^{lm}\;\chi }
215        \right)-\nabla _h \times \left( {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)}
216    \right]\;dv} \leqslant 0
217\]
218
219\[
220  % \label{eq:CONS_dynldf_curl2}
221  \mbox{if}\quad A^{lm}=cste\quad \quad \int_D {\zeta \;{\mathrm {\mathbf k}}\cdot
222    \nabla \times \left[ {\nabla _h \left( {A^{lm}\;\chi } \right)-\nabla _h
223        \times \left( {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)} \right]\;dv}
224  \leqslant 0
225\]
226
227\[
228  % \label{eq:CONS_dynldf_div2}
229  \mbox{if}\quad A^{lm}=cste\quad \quad \int_D {\chi \;\nabla _h \cdot \left[
230      {\nabla _h \left( {A^{lm}\;\chi } \right)-\nabla _h \times \left(
231          {A^{lm}\;\zeta \;{\mathrm {\mathbf k}}} \right)} \right]\;dv} \leqslant 0
232\]
233
234(II.4.6a) and (II.4.6b) means that the horizontal diffusion of momentum conserve both the potential vorticity and
235the divergence of the flow, while Eqs (II.4.6c) to (II.4.6e) mean that it dissipates the energy, the enstrophy and
236the square of the divergence.
237The two latter properties are only satisfied when the eddy coefficients are horizontally uniform.
238
239Using (II.1.8) and (II.1.9), and the symmetry or anti-symmetry properties of the mean and difference operators,
240it is shown that the discrete form of the lateral momentum diffusion given in
241{\S}II.2-c satisfies all the integral constraints (II.4.6) (see appendix C).
242In particular, when the eddy coefficients are horizontally uniform,
243a complete separation of vorticity and horizontal divergence fields is ensured,
244so that diffusion (dissipation) of vorticity (enstrophy) does not generate horizontal divergence
245(variance of the horizontal divergence) and \textit{vice versa}.
246When the vertical curl of the horizontal diffusion of momentum (discrete sense) is taken,
247the term associated to the horizontal gradient of the divergence is zero locally.
248When the horizontal divergence of the horizontal diffusion of momentum (discrete sense) is taken,
249the term associated to the vertical curl of the vorticity is zero locally.
250The resulting term conserves $\chi$ and dissipates $\chi^2$ when the eddy coefficient is horizontally uniform.
251
252\textbf{* vertical momentum diffusion term}
253
254As for the lateral momentum physics, the continuous form of the vertical diffusion of
255momentum satisfies following integral constraints~:
256
257conservation of momentum, dissipation of horizontal kinetic energy
258
259\[
260  % \label{eq:CONS_dynzdf_dyn}
261  \begin{aligned}
262    & \int_D {\frac{1}{e_3 }}  \frac{\partial }{\partial k}\left( \frac{A^{vm}}{e_3 }\frac{\partial {\textbf{U}}_h }{\partial k} \right) \;dv = \overrightarrow{\textbf{0}} \\
263    & \int_D \textbf{U}_h \cdot \frac{1}{e_3} \frac{\partial}{\partial k} \left( {\frac{A^{vm}}{e_3 }}{\frac{\partial \textbf{U}_h }{\partial k}} \right) \;dv \leq 0 \\
264  \end{aligned}
265\]
266conservation of vorticity, dissipation of enstrophy
267\[
268  % \label{eq:CONS_dynzdf_vor}
269  \begin{aligned}
270    & \int_D {\frac{1}{e_3 }{\mathrm {\mathbf k}}\cdot \nabla \times \left( {\frac{1}{e_3
271          }\frac{\partial }{\partial k}\left( {\frac{A^{vm}}{e_3 }\frac{\partial {\mathrm
272                  {\mathbf U}}_h }{\partial k}} \right)} \right)\;dv} =0 \\
273    & \int_D {\zeta \,{\mathrm {\mathbf k}}\cdot \nabla \times \left( {\frac{1}{e_3
274          }\frac{\partial }{\partial k}\left( {\frac{A^{vm}}{e_3 }\frac{\partial {\mathrm
275                  {\mathbf U}}_h }{\partial k}} \right)} \right)\;dv} \leq 0 \\
276  \end{aligned}
277\]
278conservation of horizontal divergence, dissipation of square of the horizontal divergence
279\[
280  % \label{eq:CONS_dynzdf_div}
281  \begin{aligned}
282    &\int_D {\nabla \cdot \left( {\frac{1}{e_3 }\frac{\partial }{\partial
283            k}\left( {\frac{A^{vm}}{e_3 }\frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}}
284          \right)} \right)\;dv} =0 \\
285    & \int_D {\chi \;\nabla \cdot \left( {\frac{1}{e_3 }\frac{\partial }{\partial
286            k}\left( {\frac{A^{vm}}{e_3 }\frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}}
287          \right)} \right)\;dv} \leq 0 \\
288  \end{aligned}
289\]
290
291In discrete form, all these properties are satisfied in $z$-coordinate (see Appendix C).
292In $s$-coordinates, only first order properties can be demonstrated,
293\ie\ the vertical momentum physics conserve momentum, potential vorticity, and horizontal divergence.
294
295%% =================================================================================================
296\subsection{Conservation properties on tracer physics}
297\label{subsec:CONS_Invariant_tra_physics}
298
299The numerical schemes used for tracer subgridscale physics are written in such a way that
300the heat and salt contents are conserved (equations in flux form, second order centred finite differences).
301As a form flux is used to compute the temperature and salinity,
302the quadratic form of these quantities (\ie\ their variance) globally tends to diminish.
303As for the advective term, there is generally no strict conservation of mass even if,
304in practice, the mass is conserved with a very good accuracy.
305
306\textbf{* lateral physics: }conservation of tracer, dissipation of tracer
307variance, i.e.
308
309\[
310  % \label{eq:CONS_traldf_t_t2}
311  \begin{aligned}
312    &\int_D \nabla\, \cdot\, \left( A^{lT} \,\Re \,\nabla \,T \right)\;dv = 0 \\
313    &\int_D \,T\, \nabla\, \cdot\, \left( A^{lT} \,\Re \,\nabla \,T \right)\;dv \leq 0 \\
314  \end{aligned}
315\]
316
317\textbf{* vertical physics: }conservation of tracer, dissipation of tracer variance, \ie
318
319\[
320  % \label{eq:CONS_trazdf_t_t2}
321  \begin{aligned}
322    & \int_D \frac{1}{e_3 } \frac{\partial }{\partial k}\left( \frac{A^{vT}}{e_3 }  \frac{\partial T}{\partial k}  \right)\;dv = 0 \\
323    & \int_D \,T \frac{1}{e_3 } \frac{\partial }{\partial k}\left( \frac{A^{vT}}{e_3 }  \frac{\partial T}{\partial k}  \right)\;dv \leq 0 \\
324  \end{aligned}
325\]
326
327Using the symmetry or anti-symmetry properties of the mean and difference operators,
328it is shown that the discrete form of tracer physics given in {\S}~II.2-c satisfies all the integral constraints
329(II.4.8) and (II.4.9) except the dissipation of the square of the tracer when non-geopotential diffusion is used
330(see appendix C).
331A discrete form of the lateral tracer physics can be derived which satisfies these last properties.
332Nevertheless, it requires a horizontal averaging of the vertical component of the lateral physics that
333prevents the use of implicit resolution in the vertical.
334It has not been implemented.
335
336\subinc{\input{../../global/epilogue}}
337
338\end{document}
Note: See TracBrowser for help on using the repository browser.