Discrete total energy conservation : vector invariant form

Total energy conservation

The discrete form of the total energy conservation, (C.3), is given by:

$\displaystyle \partial_t \left( \sum\limits_{i,j,k} \biggl\{ \frac{u^2}{2}  b_u + \frac{v^2}{2}  b_v + \rho   g   z_t  b_t \biggr\} \right)$ $\displaystyle =0$    

which in vector invariant forms, it leads to:

\begin{displaymath}\begin{split}\sum\limits_{i,j,k} \biggl\{ u  \partial_t u \;...
... \biggl\{ \rho  g \partial_t (z_t)  b_t \biggr\} \end{split}\end{displaymath} (C.7)

Substituting the discrete expression of the time derivative of the velocity either in vector invariant, leads to the discrete equivalent of the four equations (C.6).

Vorticity term (coriolis + vorticity part of the advection)

Let $ q$, located at $ f$-points, be either the relative ( $ q=\zeta / e_{3f}$), or the planetary ( $ q=f/e_{3f}$), or the total potential vorticity ( $ q=(\zeta +f) /e_{3f}$). Two discretisation of the vorticity term (ENE and EEN) allows the conservation of the kinetic energy.

Vorticity Term with ENE scheme (ln_dynvor_ene=.true.)

For the ENE scheme, the two components of the vorticity term are given by :

$\displaystyle - e_3   q \;{\textbf{k}}\times {\textbf {U}}_h \equiv \left( {{ ...
...( e_{2u} e_{3u} u \right)}^{ j+1/2}} ^{ i} \hfill  \end{array}} } \right)$    

This formulation does not conserve the enstrophy but it does conserve the total kinetic energy. Indeed, the kinetic energy tendency associated to the vorticity term and averaged over the ocean domain can be transformed as follows:

  $\displaystyle \int\limits_D - \left( e_3   q \;\textbf{k} \times \textbf{U}_h \right) \cdot \textbf{U}_h \; dv$        
  $\displaystyle \qquad \qquad {\begin{array}{*{20}l} &\equiv \sum\limits_{i,j,k} ...
... U }^{ j+1/2}\; \overline{ V }^{ i+1/2} \biggr\} \quad \equiv 0 \end{array} }$    

In other words, the domain averaged kinetic energy does not change due to the vorticity term.

Vorticity Term with EEN scheme (ln_dynvor_een=.true.)

With the EEN scheme, the vorticity terms are represented as:

\begin{equation*}\left\{ { \begin{aligned}+q e_3   v &\equiv +\frac{1}{e_{1u} ...
...e_{3u}  u \right)^{i+i_p}_{j+j_p-1/2}  \end{aligned} } \right.\end{equation*}

where the indices $ i_p$ and $ j_p$ take the following value: $ i_p = -1/2$ or $ 1/2$ and $ j_p = -1/2$ or $ 1/2$, and the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$, defined at $ T$-point, are given by:

$\displaystyle _i^j \mathbb{Q}^{i_p}_{j_p} = \frac{1}{12}  \left( q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p} \right)$ (C.9)

This formulation does conserve the total kinetic energy. Indeed,

  $\displaystyle \int\limits_D - \textbf{U}_h \cdot \left( \zeta \;\textbf{k} \times \textbf{U}_h \right) \; dv$        
$\displaystyle \equiv \sum\limits_{i,j,k}$ $\displaystyle \biggl\{ \left[ \sum_{\substack{i_p, k_p}} {^{i+1/2-i_p}_j}\math...
...}\mathbb{Q}^{i_p}_{j_p} \; U^{i+i_p}_{j+1/2-j_p} \right] V^{i}_{j+1/2} \biggr\}$        
$\displaystyle \equiv \sum\limits_{i,j,k}$ $\displaystyle \sum_{\substack{i_p, k_p}} \biggl\{   {^{i+1/2-i_p}_j}\mathbb{...
...}\mathbb{Q}^{i_p}_{j_p} \; U^{i+i_p}_{j+1/2-j_p}   V^{i}_{j+1/2}  \; \biggr\}$        

Expending the summation on $ i_p$ and $ k_p$, it becomes:

$\displaystyle \allowdisplaybreaks\equiv \sum\limits_{i,j,k}$ $\displaystyle \biggl\{   {^{i+1}_j }\mathbb{Q}^{-1/2}_{+1/2} \;V^{i+1}_{j+1/2...
...{j} - {^i_{j}\quad}\mathbb{Q}^{-1/2}_{+1/2} \; U^{i-1/2}_{j} \; V^{ i}_{j+1/2}$    
  $\displaystyle + {^{i+1}_j }\mathbb{Q}^{-1/2}_{-1/2} \; V^{i+1}_{j-1/2} \; U^{ ...
...i_{j+1} }\mathbb{Q}^{-1/2}_{-1/2} \; U^{i-1/2}_{j+1} \; V^{ i}_{j+1/2} \biggr.$        
  $\displaystyle + {^{i}_j\quad}\mathbb{Q}^{+1/2}_{+1/2} \; V^{i}_{j+1/2} \; U^{ ...
...i_{j}\quad}\mathbb{Q}^{+1/2}_{+1/2} \; U^{i+1/2}_{j} \; V^{ i}_{j+1/2} \biggr.$        
  $\displaystyle + {^{i}_j\quad}\mathbb{Q}^{+1/2}_{-1/2} \; V^{i}_{j-1/2} \; U^{ ...
...1} }\mathbb{Q}^{+1/2}_{-1/2} \; U^{i+1/2}_{j+1}\; V^{ i}_{j+1/2}  \; \biggr\}$        

The summation is done over all $ i$ and $ j$ indices, it is therefore possible to introduce a shift of $ -1$ either in $ i$ or $ j$ direction in some of the term of the summation (first term of the first and second lines, second term of the second and fourth lines). By doning so, we can regroup all the terms of the summation by triad at a ($ i$,$ j$) point. In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($ i$,$ j$) indices. It becomes:

$\displaystyle \allowdisplaybreaks\allowdisplaybreaks \equiv \sum\limits_{i,j,k}$ $\displaystyle \biggl\{   {^{i}_j}\mathbb{Q}^{-1/2}_{+1/2} \left[ V^{i}_{j+1/2}  U^{ i-1/2}_{j} - U^{i-1/2}_{j}   V^{ i}_{j+1/2} \right]$    
  $\displaystyle + {^{i}_j}\mathbb{Q}^{-1/2}_{-1/2} \left[ V^{i}_{j-1/2}   U^{ i-1/2}_{j} - U^{i-1/2}_{j}   V^{ i}_{j-1/2} \right] \biggr.$        
  $\displaystyle + {^{i}_j}\mathbb{Q}^{+1/2}_{+1/2} \left[ V^{i}_{j+1/2}   U^{ i+1/2}_{j} - U^{i+1/2}_{j}   V^{ i}_{j+1/2} \right] \biggr.$        
  $\displaystyle + {^{i}_j}\mathbb{Q}^{+1/2}_{-1/2} \left[ V^{i}_{j-1/2}   U^{ i...
...j} - U^{i+1/2}_{j-1}   V^{ i}_{j-1/2} \right]  \; \biggr\} \qquad \equiv  0$        

Gradient of Kinetic Energy / Vertical Advection

The change of Kinetic Energy (KE) due to the vertical advection is exactly balanced by the change of KE due to the horizontal gradient of KE :

$\displaystyle \int_D \textbf{U}_h \cdot \frac{1}{e_3 } \omega \partial_k \textb...
... \frac{1}{2} \int_D { \frac{{\textbf{U}_h}^2}{e_3} \partial_t ( e_3) \;dv }  $    

Indeed, using successively (4.11) ($ i.e.$ the skew symmetry property of the $ \delta$ operator) and the continuity equation, then (4.11) again, then the commutativity of operators $ \overline { \cdot  }$ and $ \delta$, and finally (4.12) ($ i.e.$ the symmetry property of the $ \overline { \cdot  }$ operator) applied in the horizontal and vertical directions, it becomes:

  $\displaystyle - \int_D \textbf{U}_h \cdot$   KEG$\displaystyle \;dv = - \int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv$    
$\displaystyle \equiv$ $\displaystyle - \sum\limits_{i,j,k} \frac{1}{2} \biggl\{ \frac{1} {e_{1u}} \del...
...} \left[ \overline {u^2}^{ i} + \overline {v^2}^{ j} \right] v  b_v \biggr\}$          
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} \frac{1}{2} \left( \overline {u^2}^{ i} + ...
...)\; \biggl\{ \delta_{i} \left[ U \right] + \delta_{j} \left[ V \right] \biggr\}$          
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle - \sum\limits_{i,j,k} \frac{1}{2} \left( \overline {u^2}^{ i} + ...
...l\{ \frac{b_t}{e_{3t}} \partial_t (e_{3t}) + \delta_k \left[ W \right] \biggr\}$    
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle + \sum\limits_{i,j,k} \frac{1}{2} \delta_{k+1/2} \left[ \overline...
...} \left( \overline {u^2}^{ i} + \overline {v^2}^{ j} \right) \;\partial_t b_t$          
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle + \sum\limits_{i,j,k} \frac{1} {2} \left( \overline{\delta_{k+1/2...
...b_t}^{ {i+1/2}} + \frac{v^2}{2} \partial_t \overline{b_t}^{ {j+1/2}} \right)$          

Assuming that $ b_u= \overline{b_t}^{ i+1/2}$ and $ b_v= \overline{b_t}^{ j+1/2}$, or at least that the time derivative of these two equations is satisfied, it becomes:

$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle \sum\limits_{i,j,k} \frac{1} {2} \biggl\{ \; \overline{W}^{ i+1/...
...k} \left( \frac{u^2}{2} \partial_t b_u + \frac{v^2}{2} \partial_t b_v \right)$    
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle \sum\limits_{i,j,k} \biggl\{ \; \overline{W}^{ i+1/2}\; \overlin...
...k} \left( \frac{u^2}{2} \partial_t b_u + \frac{v^2}{2} \partial_t b_v \right)$          
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle \sum\limits_{i,j,k} \biggl\{ \; \frac{1} {b_u } \; \overline { \o...
...k} \left( \frac{u^2}{2} \partial_t b_u + \frac{v^2}{2} \partial_t b_v \right)$          

The first term provides the discrete expression for the vertical advection of momentum (ZAD), while the second term corresponds exactly to (C.7), therefore:

$\displaystyle \equiv$ $\displaystyle \int\limits_D \textbf{U}_h \cdot$   ZAD$\displaystyle \;dv + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t (e_3) \;dv }$    
$\displaystyle \equiv$ $\displaystyle \int\limits_D \textbf{U}_h \cdot w \partial_k \textbf{U}_h \;dv + \frac{1}{2} \int_D { {\textbf{U}_h}^2 \frac{1}{e_3} \partial_t (e_3) \;dv }$    

There is two main points here. First, the satisfaction of this property links the choice of the discrete formulation of the vertical advection and of the horizontal gradient of KE. Choosing one imposes the other. For example KE can also be discretized as $ 1/2 ({\overline u^{ i}}^2 + {\overline v^{ j}}^2)$. This leads to the following expression for the vertical advection:

$\displaystyle \frac{1} {e_3 }\; \omega\; \partial_k \textbf{U}_h \equiv \left( ...
...t[ \overline v^{ j+1/2} \right]}}^{ j+1/2,k} \hfill  \end{array}} } \right)$    

a formulation that requires an additional horizontal mean in contrast with the one used in NEMO. Nine velocity points have to be used instead of 3. This is the reason why it has not been chosen.

Second, as soon as the chosen $ s$-coordinate depends on time, an extra constraint arises on the time derivative of the volume at $ u$- and $ v$-points:

$\displaystyle e_{1u} e_{2u} \partial_t (e_{3u}) =\overline{ e_{1t} e_{2t}\;\partial_t (e_{3t}) }^{ i+1/2}$    
$\displaystyle e_{1v} e_{2v} \partial_t (e_{3v}) =\overline{ e_{1t} e_{2t}\;\partial_t (e_{3t}) }^{ j+1/2}$    

which is (over-)satified by defining the vertical scale factor as follows:

$\displaystyle e_{3u} = \frac{1}{e_{1u} e_{2u}}\;\overline{ e_{1t}^{ } e_{2t}^{ } e_{3t}^{ } }^{ i+1/2}$ (C.10)
$\displaystyle e_{3v} = \frac{1}{e_{1v} e_{2v}}\;\overline{ e_{1t}^{ } e_{2t}^{ } e_{3t}^{ } }^{ j+1/2}$ (C.11)

Blah blah required on the the step representation of bottom topography.....

Pressure Gradient Term

When the equation of state is linear ($ i.e.$ when an advection-diffusion equation for density can be derived from those of temperature and salinity) the change of KE due to the work of pressure forces is balanced by the change of potential energy due to buoyancy forces:

$\displaystyle - \int_D \left. \nabla p \right\vert _z \cdot \textbf{U}_h \;dv =...
...\rho  \textbf {U} \right)  g z \;dv + \int_D g  \rho \; \partial_t (z) \;dv$    

This property can be satisfied in a discrete sense for both $ z$- and $ s$-coordinates. Indeed, defining the depth of a $ T$-point, $ z_t$, as the sum of the vertical scale factors at $ w$-points starting from the surface, the work of pressure forces can be written as:

  $\displaystyle - \int_D \left. \nabla p \right\vert _z \cdot \textbf{U}_h \;dv \...
...1/2} [p_t] - g\;\overline \rho^{ i+1/2}\;\delta_{i+1/2} [z_t] \Bigr) \; u\;b_u$        
  $\displaystyle \qquad \qquad \qquad \qquad \qquad \quad    - \frac{1} {e_{2v}}...
... - g\;\overline \rho^{ j+1/2}\delta_{j+1/2} [z_t] \Bigr) \; v\;b_v \; \biggr\}$        

Using successively (4.11), $ i.e.$ the skew symmetry property of the $ \delta$ operator, (6.4), the continuity equation, (6.20), the hydrostatic equation in the $ s$-coordinate, and $ \delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} $, which comes from the definition of $ z_t$, it becomes:

$\displaystyle \allowdisplaybreaks\allowdisplaybreaks \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \biggl\{ \overline\rho^{ i+1/2} U \del...
..._{j+1/2}[z_t] +\Bigl( \delta_i[U] + \delta_j [V] \Bigr)\;\frac{p_t}{g} \biggr\}$    
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \biggl\{ \overline\rho^{ i+1/2} U \del...
... \partial_t (e_{3t}) + \delta_k \left[ W \right] \right) \frac{p_t}{g} \biggr\}$    
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \biggl\{ \overline\rho^{ i+1/2} U \del...
...t] + \frac{W}{g}\;\delta_{k+1/2} [p_t] - \frac{p_t}{g} \partial_t b_t \biggr\}$    
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \biggl\{ \overline\rho^{ i+1/2} U \del...
...] - W\;e_{3w} \overline \rho^{ k+1/2} - \frac{p_t}{g} \partial_t b_t \biggr\}$    
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \biggl\{ \overline\rho^{ i+1/2} U \del...
...e \rho^{ k+1/2}\;\delta_{k+1/2} [z_t] - \frac{p_t}{g} \partial_t b_t \biggr\}$    
$\displaystyle \allowdisplaybreaks \equiv$ $\displaystyle - \sum\limits_{i,j,k} g \; z_t \biggl\{ \delta_i \left[ U\; \over...
...2} \right] \biggr\} - \sum\limits_{i,j,k} \biggl\{ p_t\;\partial_t b_t \biggr\}$    
$\displaystyle \equiv$ $\displaystyle + \sum\limits_{i,j,k} g \; z_t \biggl\{ \partial_t ( e_{3t}  \rho) \biggr\} \; b_t - \sum\limits_{i,j,k} \biggl\{ p_t\;\partial_t b_t \biggr\}$    

The first term is exactly the first term of the right-hand-side of (C.7). It remains to demonstrate that the last term, which is obviously a discrete analogue of $ \int_D \frac{p}{e_3} \partial_t (e_3)\;dv$ is equal to the last term of (C.7). In other words, the following property must be satisfied:

$\displaystyle \sum\limits_{i,j,k} \biggl\{ p_t\;\partial_t b_t \biggr\} \equiv \sum\limits_{i,j,k} \biggl\{ \rho  g \partial_t (z_t)  b_t \biggr\}$    

Let introduce $ p_w$ the pressure at $ w$-point such that $ \delta_k [p_w] = - \rho  g e_{3t}$. The right-hand-side of the above equation can be transformed as follows:

$\displaystyle \sum\limits_{i,j,k} \biggl\{ \rho  g \partial_t (z_t)  b_t \biggr\}$ $\displaystyle \equiv - \sum\limits_{i,j,k} \biggl\{ \delta_k [p_w] \partial_t (z_t)  e_{1t} e_{2t} \biggr\}$    
  $\displaystyle \equiv + \sum\limits_{i,j,k} \biggl\{ p_w  \delta_{k+1/2} [\part...
...sum\limits_{i,j,k} \biggl\{ p_w  \partial_t (e_{3w})  e_{1t} e_{2t} \biggr\}$    
  $\displaystyle \equiv + \sum\limits_{i,j,k} \biggl\{ p_w  \partial_t (b_w) \biggr\} %

therefore, the balance to be satisfied is:

$\displaystyle \sum\limits_{i,j,k} \biggl\{ p_t\;\partial_t (b_t) \biggr\} \equiv \sum\limits_{i,j,k} \biggl\{ p_w  \partial_t (b_w) \biggr\}$    

which is a purely vertical balance:

$\displaystyle \sum\limits_{k} \biggl\{ p_t\;\partial_t (e_{3t}) \biggr\} \equiv \sum\limits_{k} \biggl\{ p_w  \partial_t (e_{3w}) \biggr\}$    

Defining $ p_w = \overline{p_t}^{ k+1/2}$

Note that this property strongly constrains the discrete expression of both the depth of $ T-$points and of the term added to the pressure gradient in the $ s$-coordinate. Nevertheless, it is almost never satisfied since a linear equation of state is rarely used.

Gurvan Madec and the NEMO Team
NEMO European Consortium2017-02-17