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

Last change on this file since 10442 was 10442, checked in by nicolasmartin, 2 years ago

Front page edition, cleaning in custom LaTeX commands and add index for single subfile compilation

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