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/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_conservation.tex @ 10354

Last change on this file since 10354 was 10354, checked in by nicolasmartin, 5 years ago

Vast edition of LaTeX subfiles to improve the readability by cutting sentences in a more suitable way
Every sentence begins in a new line and if necessary is splitted around 110 characters lenght for side-by-side visualisation,
this setting may not be adequate for everyone but something has to be set.
The punctuation was the primer trigger for the cutting process, otherwise subordinators and coordinators, in order to mostly keep a meaning for each line

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