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 @ 11544

Last change on this file since 11544 was 11544, checked in by nicolasmartin, 3 years ago

Missing modifications from previous commit

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