Discrete total energy conservation : vector invariant form

Total energy conservation

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

$$\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) =0$$

which in vector invariant forms, it leads to:

$$\begin{split}\sum\limits_{i,j,k} \biggl\{ u \partial_t u \;... ... \biggl\{ \rho g \partial_t (z_t) b_t \biggr\} \end{split}$$ (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 :

$$- 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:

$$\int\limits_D - \left( e_3 q \;\textbf{k} \times \textbf{U}_h \right) \cdot \textbf{U}_h \; dv$$

$$\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:

$$\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.$$

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:

$$_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,

$$\int\limits_D - \textbf{U}_h \cdot \left( \zeta \;\textbf{k} \times \textbf{U}_h \right) \; dv$$

$$\equiv \sum\limits_{i,j,k} \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\}$$

$$\equiv \sum\limits_{i,j,k} \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\}$$

$$i_p k_p$$

$$\allowdisplaybreaks\equiv \sum\limits_{i,j,k} \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}$$

$$+ {^{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.$$

$$+ {^{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.$$

$$+ {^{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\}$$

$$i j -1 i j i j i j$$

$$\allowdisplaybreaks\allowdisplaybreaks \equiv \sum\limits_{i,j,k} \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]$$

$$+ {^{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.$$

$$+ {^{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.$$

$$+ {^{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 :

$$\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:

$$- \int_D \textbf{U}_h \cdot \;dv = - \int_D \textbf{U}_h \cdot \nabla_h \left( \frac{1}{2}\;{\textbf{U}_h}^2 \right)\;dv$$

$$\equiv - \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\}$$

$$\equiv + \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\}$$

$$\allowdisplaybreaks \equiv - \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\}$$

$$\allowdisplaybreaks \equiv + \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$$

$$\allowdisplaybreaks \equiv + \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)$$

$$b_u= \overline{b_t}^{ i+1/2} b_v= \overline{b_t}^{ j+1/2}$$

$$\allowdisplaybreaks \equiv \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)$$

$$\allowdisplaybreaks \equiv \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)$$

$$\allowdisplaybreaks \equiv \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:

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

$$\equiv \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:

$$\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:

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

$$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:

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

$$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:

$$- \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:

$$- \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$$

$$\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\}$$

$$i.e. \delta s \delta_{k+1/2} \left[ z_t \right] \equiv e_{3w} z_t$$

$$\allowdisplaybreaks\allowdisplaybreaks \equiv + \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\}$$

$$\equiv + \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\}$$

$$\equiv + \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\}$$

$$\equiv + \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\}$$

$$\equiv + \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\}$$

$$\allowdisplaybreaks \equiv - \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\}$$

$$\equiv + \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:

$$\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:

$$\sum\limits_{i,j,k} \biggl\{ \rho g \partial_t (z_t) b_t \biggr\} \equiv - \sum\limits_{i,j,k} \biggl\{ \delta_k [p_w] \partial_t (z_t) e_{1t} e_{2t} \biggr\}$$

$$\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\}$$

$$\equiv + \sum\limits_{i,j,k} \biggl\{ p_w \partial_t (b_w) \biggr\} %$$

therefore, the balance to be satisfied is:

$$\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:

$$\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.