--- manual/s_overview/continuous_eqns.tex 2001/08/08 16:16:18 1.1 +++ manual/s_overview/continuous_eqns.tex 2001/09/27 01:57:17 1.4 @@ -1,18 +1,18 @@ -% $Header: /home/ubuntu/mnt/e9_copy/manual/s_overview/Attic/continuous_eqns.tex,v 1.1 2001/08/08 16:16:18 adcroft Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/manual/s_overview/Attic/continuous_eqns.tex,v 1.4 2001/09/27 01:57:17 cnh Exp $ % $Name: $ \section{Continuous equations in `r' coordinates} To render atmosphere and ocean models from one dynamical core we exploit `isomorphisms' between equation sets that govern the evolution of the -respective fluids - see fig.4% +respective fluids - see fig.4 \marginpar{ Fig.4. Isomorphisms}. One system of hydrodynamical equations is written down and encoded. The model variables have different interpretations depending on whether the atmosphere or ocean is being studied. Thus, for example, the vertical coordinate `$r$' is interpreted as pressure, $p$, if we are modeling the atmosphere and height, $z$, if we are modeling the ocean. A -complete list of the isomorphisms is given in table 1.% +complete list of the isomorphisms is given in table 1. \marginpar{ Table 1. Isomorphisms} @@ -22,24 +22,24 @@ depend on $\theta $, $S$, and $p$. The equations that govern the evolution of these fields, obtained by applying the laws of classical mechanics and thermodynamics to a Boussinesq, Navier-Stokes fluid are, written in terms of -a generic vertical coordinate, $r$, see fig.5% +a generic vertical coordinate, $r$, see fig.5 \marginpar{ Fig.5 The vertical coordinate of model}: \[ -\frac{D\vec{\mathbf{v}_{h}}}{Dt}+\left( 2\vec{\Omega}\times \vec{\mathbf{v}}% -\right) _{h}+\mathbf{\nabla }_{h}\phi =\left( \mathcal{F}_{\vec{\mathbf{v}}}% +\frac{D\vec{\mathbf{v}_{h}}}{Dt}+\left( 2\vec{\Omega}\times \vec{\mathbf{v}} +\right) _{h}+\mathbf{\nabla }_{h}\phi =\left( \mathcal{F}_{\vec{\mathbf{v}}} \mathcal{+D}_{\vec{\mathbf{v}}}\right) _{h}\text{horizontal mtm} \] \[ -\frac{D\dot{r}}{Dt}+\widehat{k}\cdot \left( 2\vec{\Omega}\times \vec{\mathbf{% -v}}\right) +\frac{\partial \phi }{\partial r}+b=\left( \mathcal{F}_{\vec{% +\frac{D\dot{r}}{Dt}+\widehat{k}\cdot \left( 2\vec{\Omega}\times \vec{\mathbf{ +v}}\right) +\frac{\partial \phi }{\partial r}+b=\left( \mathcal{F}_{\vec{ \mathbf{v}}}\mathcal{+D}_{\vec{\mathbf{v}}}\right) _{r}\text{vertical mtm} \] \begin{equation} -\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \dot{r}}{% +\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \dot{r}}{ \partial r}=0\text{ continuity} \label{incompressible} \end{equation} @@ -53,7 +53,7 @@ \] \[ -\frac{DS}{Dt}=\mathcal{F}_{S}\text{ }\mathcal{+D}_{S}\text{ humidity/salinity% +\frac{DS}{Dt}=\mathcal{F}_{S}\text{ }\mathcal{+D}_{S}\text{ humidity/salinity } \] @@ -69,10 +69,10 @@ \] \[ -\mathbf{\nabla }=\mathbf{\nabla }_{h}+\widehat{k}\frac{\partial }{\partial r}% +\mathbf{\nabla }=\mathbf{\nabla }_{h}+\widehat{k}\frac{\partial }{\partial r} \text{ is the `grad' operator} \] -with $\mathbf{\nabla }_{h}$ operating in the horizontal and $\widehat{k}% +with $\mathbf{\nabla }_{h}$ operating in the horizontal and $\widehat{k} \frac{\partial }{\partial r}$ operating in the vertical, where $\widehat{k}$ is a unit vector in the vertical @@ -106,7 +106,7 @@ \] \[ -\mathcal{F}_{\vec{\mathbf{v}}}\text{ and }\mathcal{D}_{\vec{\mathbf{v}}}% +\mathcal{F}_{\vec{\mathbf{v}}}\text{ and }\mathcal{D}_{\vec{\mathbf{v}}} \text{ are forcing and dissipation of }\vec{\mathbf{v}} \] @@ -201,11 +201,13 @@ The boundary conditions at top and bottom are given by: -\begin{eqnarray*} -&&\omega =0~\text{at }r=R_{fixed}\text{ (top of the atmosphere)} \\ +\begin{eqnarray} +&&\omega =0~\text{at }r=R_{fixed} \label{eq:fixed-bc-atmos} +\text{ (top of the atmosphere)} \\ \omega &=&\frac{Dp_{s}}{Dt}\text{; at }r=R_{moving}\text{ (bottom of the atmosphere)} -\end{eqnarray*} +\label{eq:moving-bc-atmos} +\end{eqnarray} Then the (hydrostatic form of) eq(\ref{incompressible}) yields a consistent set of atmospheric equations which, for convenience, are written out in $p$ @@ -216,11 +218,15 @@ In the ocean we interpret: \begin{eqnarray} -r &=&z\text{ is the height} \\ -\dot{r} &=&\frac{Dz}{Dt}=w\text{ is the vertical velocity} \\ -\phi &=&\frac{p}{\rho _{c}}\text{ is the pressure} \\ +r &=&z\text{ is the height} +\label{eq:ocean-z}\\ +\dot{r} &=&\frac{Dz}{Dt}=w\text{ is the vertical velocity} +\label{eq:ocean-w}\\ +\phi &=&\frac{p}{\rho _{c}}\text{ is the pressure} +\label{eq:ocean-p}\\ b(\theta ,S,r) &=&\frac{g}{\rho _{c}}\left( \rho (\theta ,S,r)-\rho _{c}\right) \text{ is the buoyancy} +\label{eq:ocean-b} \end{eqnarray} where $\rho _{c}$ is a fixed reference density of water and $g$ is the acceleration due to gravity.\noindent @@ -231,15 +237,17 @@ The surface of the ocean is given by: $R_{moving}=\eta $ -The position of the resting free surface of the ocean is given by $% +The position of the resting free surface of the ocean is given by $ R_{o}=Z_{o}=0$. Boundary conditions are: -\begin{eqnarray*} -w &=&0~\text{at }r=R_{fixed}\text{ (ocean bottom)} \\ -w &=&\frac{D\eta }{Dt}\text{ at }r=R_{moving}=\eta \text{ (ocean surface)} -\end{eqnarray*} +\begin{eqnarray} +w &=&0~\text{at }r=R_{fixed}\text{ (ocean bottom)} +\label{eq:fixed-bc-ocean}\\ +w &=&\frac{D\eta }{Dt}\text{ at }r=R_{moving}=\eta \text{ (ocean surface) +\label{eq:moving-bc-ocean}} +\end{eqnarray} where $\eta $ is the elevation of the free surface. Then eq(\ref{incompressible}) yields a consistent set of oceanic equations @@ -252,30 +260,30 @@ \begin{equation} \phi (x,y,r)=\phi _{s}(x,y)+\phi _{hyd}(x,y,r)+\phi _{nh}(x,y,r) -\label{pressuresplit} +\label{eq:phi-split} \end{equation} and write eq(\ref{incompressible}a) in the form: \begin{equation} \frac{\partial \vec{\mathbf{v}_{h}}}{\partial t}+\mathbf{\nabla }_{h}\phi _{s}+\mathbf{\nabla }_{h}\phi _{hyd}+\epsilon _{nh}\mathbf{\nabla }_{h}\phi -_{nh}=\vec{\mathbf{G}}_{\vec{v}_{h}} \label{hor-mtm} +_{nh}=\vec{\mathbf{G}}_{\vec{v}_{h}} \label{eq:mom-h} \end{equation} \begin{equation} -\frac{\partial \phi _{hyd}}{\partial r}=-b \label{hydro} +\frac{\partial \phi _{hyd}}{\partial r}=-b \label{eq:hydrostatic} \end{equation} \begin{equation} -\epsilon _{nh}\frac{\partial \dot{r}}{\partial t}+\frac{\partial \phi _{nh}}{% -\partial r}=G_{\dot{r}} \label{vertmtm} +\epsilon _{nh}\frac{\partial \dot{r}}{\partial t}+\frac{\partial \phi _{nh}}{ +\partial r}=G_{\dot{r}} \label{eq:mom-w} \end{equation} Here $\epsilon _{nh}$ is a non-hydrostatic parameter. The $\left( \vec{\mathbf{G}}_{\vec{v}},G_{\dot{r}}\right) $ in eq(\ref {hor-mtm}) and (\ref{vertmtm}) represent advective, metric and Coriolis -terms in the momentum equations. In spherical coordinates they take the form% -\footnote{% +terms in the momentum equations. In spherical coordinates they take the form +\footnote{ In the hydrostatic primitive equations (\textbf{HPE}) all underlined terms in (\ref{Gu}), (\ref{Gv}) and (\ref{Gw}) are omitted; the singly-underlined terms are included in the quasi-hydrostatic model (\textbf{QH}). The fully @@ -285,11 +293,11 @@ \left. \begin{tabular}{l} $G_{u}=-\vec{\mathbf{v}}.\nabla u$ \\ -$-\left\{ \underline{\frac{u\dot{r}}{{r}}}-\frac{uv\tan lat}{{r}}% +$-\left\{ \underline{\frac{u\dot{r}}{{r}}}-\frac{uv\tan lat}{{r}} \right\} $ \\ -$-\left\{ -2\Omega v\sin lat+\underline{\underline{2\Omega \dot{r}\cos lat}}% +$-\left\{ -2\Omega v\sin lat+\underline{\underline{2\Omega \dot{r}\cos lat}} \right\} $ \\ -$+\mathcal{F}_{u}\mathcal{+D}_{u}$% +$+\mathcal{F}_{u}\mathcal{+D}_{u}$ \end{tabular} \right\} \left\{ \begin{tabular}{l} @@ -298,17 +306,17 @@ \textit{Coriolis} \\ \textit{\ Forcing/Dissipation} \end{tabular} -\right. \qquad \label{Gu} +\right. \qquad \label{eq:gu-speherical} \end{equation} \begin{equation} \left. \begin{tabular}{l} $G_{v}=-\vec{\mathbf{v}}.\nabla v$ \\ -$-\left\{ \underline{\frac{v\dot{r}}{{r}}}-\frac{u^{2}\tan lat}{{r}% +$-\left\{ \underline{\frac{v\dot{r}}{{r}}}-\frac{u^{2}\tan lat}{{r} }\right\} $ \\ $-\left\{ -2\Omega u\sin lat\right\} $ \\ -$+\mathcal{F}_{v}\mathcal{+D}_{v}$% +$+\mathcal{F}_{v}\mathcal{+D}_{v}$ \end{tabular} \right\} \left\{ \begin{tabular}{l} @@ -317,7 +325,7 @@ \textit{Coriolis} \\ \textit{\ Forcing/Dissipation} \end{tabular} -\right. \qquad \label{Gv} +\right. \qquad \label{eq:gv-spherical} \end{equation} \qquad \qquad \qquad \qquad \qquad @@ -325,10 +333,10 @@ \left. \begin{tabular}{l} $G_{\dot{r}}=-\vec{\mathbf{v}}.\nabla \dot{r}$ \\ -$+\left\{ \frac{u^{_{^{2}}}+v^{2}}{{{r}}}% +$+\left\{ \frac{u^{_{^{2}}}+v^{2}}{{{r}}} \right\} $ \\ ${+2\Omega u\cos lat}$ \\ -$\mathcal{F}_{\dot{r}}\mathcal{+D}_{\dot{r}}$% +$\mathcal{F}_{\dot{r}}\mathcal{+D}_{\dot{r}}$ \end{tabular} \right\} \left\{ \begin{tabular}{l} @@ -337,15 +345,15 @@ \textit{Coriolis} \\ \textit{\ Forcing/Dissipation} \end{tabular} -\right. \label{Gw} +\right. \label{eq:gw-spherical} \end{equation} \qquad \qquad \qquad \qquad \qquad -In the above `${r}$' is the distance from the center of the earth and `$% +In the above `${r}$' is the distance from the center of the earth and `$ lat$' is latitude. Grad and div operators in spherical coordinates are defined in appendix -OPERATORS.% +OPERATORS. \marginpar{ Fig.6 Spherical polar coordinate system.} @@ -358,7 +366,7 @@ These are discussed at length in Marshall et al (1997a). In the `hydrostatic primitive equations' (\textbf{HPE)} all the underlined -terms in Eqs. (\ref{Gu} $\rightarrow $\ \ref{Gw}) are neglected and `${r% +terms in Eqs. (\ref{Gu} $\rightarrow $\ \ref{Gw}) are neglected and `${r }$' is replaced by `$a$', the mean radius of the earth. Once the pressure is found at one level - e.g. by inverting a 2-d Elliptic equation for $\phi _{s} $ at $r=R_{moving}$ - the pressure can be computed at all other levels @@ -367,7 +375,7 @@ In the `quasi-hydrostatic' equations (\textbf{QH)} strict balance between gravity and vertical pressure gradients is not imposed. The $2\Omega u\cos \phi $ Coriolis term are not neglected and are balanced by a non-hydrostatic -contribution to the pressure field: only the terms underlined twice in Eqs. (% +contribution to the pressure field: only the terms underlined twice in Eqs. ( \ref{Gu} $\rightarrow $\ \ref{Gw}) are set to zero and, simultaneously, the shallow atmosphere approximation is relaxed. In \textbf{QH}\ \textit{all} the metric terms are retained and the full variation of the radial position @@ -391,7 +399,7 @@ \paragraph{Non-hydrostatic Ocean} -In the non-hydrostatic ocean model all terms in equations (\ref{Gu} $% +In the non-hydrostatic ocean model all terms in equations (\ref{Gu} $ \rightarrow $\ \ref{Gw}) are retained. A three dimensional elliptic equation must be solved subject to Neumann boundary conditions (see below). It is important to note that use of the full \textbf{NH} does not admit any new @@ -404,12 +412,12 @@ \paragraph{Quasi-nonhydrostatic Atmosphere} -In the non-hydrostatic version of our atmospheric model we approximate $\dot{% +In the non-hydrostatic version of our atmospheric model we approximate $\dot{ r}$ in the vertical momentum eqs(\ref{vertmtm}) and (\ref{Gw}) (but only here) by: \begin{equation} -\dot{r}=\frac{Dp}{Dt}=\frac{Dp_{hyd}}{Dt} \label{quasinonhydro} +\dot{r}=\frac{Dp}{Dt}=\frac{Dp_{hyd}}{Dt} \label{eq:quasi-nh-w} \end{equation} where $p_{hy}$ is the hydrostatic pressure. @@ -440,20 +448,20 @@ \subparagraph{Non-hydrostatic } -Non-hydrostatic forms of the incompressible Boussinesq equations in $z-$% +Non-hydrostatic forms of the incompressible Boussinesq equations in $z-$ coordinates are supported. \subsection{Solution strategy} -The method of solution employed in the \textbf{HPE}, \textbf{QH} and \textbf{% -NH} models are summarized in Fig.7.% +The method of solution employed in the \textbf{HPE}, \textbf{QH} and \textbf{ +NH} models are summarized in Fig.7. \marginpar{ Fig.7 Solution strategy} Overview paragraph...... There is no penalty in implementing \textbf{QH} over \textbf{HPE} except, of -course, some complication that goes with the inclusion of $\cos \phi \ $% +course, some complication that goes with the inclusion of $\cos \phi \ $ Coriolis terms and the relaxation of the shallow atmosphere approximation. But this leads to negligible increase in computation. In \textbf{NH}, in contrast, one additional elliptic equation - a three-dimensional one - must @@ -469,13 +477,8 @@ dividing the total (pressure/geo) potential in to three parts, a surface part, $\phi _{s}(x,y)$, a hydrostatic part $\phi _{hyd}(x,y,r)$ and a non-hydrostatic part $\phi _{nh}(x,y,r)$, as in (\ref{pressuresplit}), and -writing the momentum equation in the form -\begin{equation} -\frac{\partial }{\partial t}\vec{\mathbf{v}_{h}}+\mathbf{\nabla }_{h}\phi -_{s}+\mathbf{\nabla }_{h}\phi _{hyd}+\epsilon _{nh}\mathbf{\nabla }\phi -_{nh}=\vec{\mathbf{G}}_{\vec{v}} \label{mtm-split} -\end{equation} -as in (\ref{hor-mtm}). +writing the momentum equation +as in (\ref{eq:mom-h}). \subsubsection{Hydrostatic pressure} @@ -488,17 +491,18 @@ \] and so -\[ +\begin{equation} \phi _{hyd}(x,y,r)=\int_{r}^{R_{o}}bdr -\] +\label{eq:hydro-phi} +\end{equation} \subsubsection{Surface pressure} -The surface pressure equation can be obtained by integrating continuity, (% +The surface pressure equation can be obtained by integrating continuity, ( \ref{incompressible})c, vertically from $r=R_{fixed}$ to $r=R_{moving}$ \[ -\int_{R_{fixed}}^{R_{moving}}\left( \mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}% +\int_{R_{fixed}}^{R_{moving}}\left( \mathbf{\nabla }_{h}\cdot \vec{\mathbf{v} }_{h}+\partial _{r}\dot{r}\right) dr=0 \] @@ -506,23 +510,24 @@ \[ \frac{\partial \eta }{\partial t}+\vec{\mathbf{v}}.\nabla \eta -+\int_{R_{fixed}}^{R_{moving}}\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}}% ++\int_{R_{fixed}}^{R_{moving}}\mathbf{\nabla }_{h}\cdot \vec{\mathbf{v}} _{h}dr=0 \] -where $\eta =R_{moving}-R_{o}$ is the free-surface $r$-anomaly in units of $% +where $\eta =R_{moving}-R_{o}$ is the free-surface $r$-anomaly in units of $ r $. The above can be rearranged to yield, using Leibnitz's theorem: \begin{equation} \frac{\partial \eta }{\partial t}+\mathbf{\nabla }_{h}\cdot \int_{R_{fixed}}^{R_{moving}}\vec{\mathbf{v}}_{h}dr=0 -\label{integralcontinuity} +\label{eq:free-surface} \end{equation} Whether $\phi $ is pressure (ocean model, $p/\rho _{c}$) or geopotential (atmospheric model), in (\ref{mtm-split}), the horizontal gradient term can be written \begin{equation} -\mathbf{\nabla }_{h}\phi _{s}=\mathbf{\nabla }_{h}b\eta \label{link} +\mathbf{\nabla }_{h}\phi _{s}=\mathbf{\nabla }_{h}b\eta +\label{eq:phi-surf} \end{equation} where $b$ is the buoyancy. @@ -533,14 +538,14 @@ \subsubsection{Non-hydrostatic pressure} -Taking the horizontal divergence of (\ref{hor-mtm}) and adding $\frac{% +Taking the horizontal divergence of (\ref{hor-mtm}) and adding $\frac{ \partial }{\partial r}$ of (\ref{vertmtm}), invoking the continuity equation (\ref{incompressible}), we deduce that: \begin{equation} -\nabla _{3}^{2}\phi _{nh}=\nabla .\vec{\mathbf{G}}_{\vec{v}}-\left( \mathbf{% -\nabla }_{h}^{2}\phi _{s}+\mathbf{\nabla }^{2}\phi _{hyd}\right) =\nabla .% -\vec{\mathbf{F}} \label{3dinvert} +\nabla _{3}^{2}\phi _{nh}=\nabla .\vec{\mathbf{G}}_{\vec{v}}-\left( \mathbf{ +\nabla }_{h}^{2}\phi _{s}+\mathbf{\nabla }^{2}\phi _{hyd}\right) =\nabla . +\vec{\mathbf{F}} \label{eq:3d-invert} \end{equation} For a given rhs this 3-d elliptic equation must be inverted for $\phi _{nh}$ @@ -559,7 +564,7 @@ \end{equation} where $\widehat{n}$ is a vector of unit length normal to the boundary. The kinematic condition (\ref{nonormalflow}) is also applied to the vertical -velocity at $r=R_{moving}$. No-slip $\left( v_{T}=0\right) \ $or slip $% +velocity at $r=R_{moving}$. No-slip $\left( v_{T}=0\right) \ $or slip $ \left( \partial v_{T}/\partial n=0\right) \ $conditions are employed on the tangential component of velocity, $v_{T}$, at all solid boundaries, depending on the form chosen for the dissipative terms in the momentum @@ -569,7 +574,7 @@ \begin{equation} \widehat{n}.\nabla \phi _{nh}=\widehat{n}.\vec{\mathbf{F}} -\label{inhomneumann} +\label{eq:inhom-neumann-nh} \end{equation} where @@ -579,7 +584,7 @@ \] presenting inhomogeneous Neumann boundary conditions to the Elliptic problem (\ref{3dinvert}). As shown, for example, by Williams (1969), one can exploit -classical 3D potential theory and, by introducing an appropriately chosen $% +classical 3D potential theory and, by introducing an appropriately chosen $ \delta $-function sheet of `source-charge', replace the inhomogenous boundary condition on pressure by a homogeneous one. The source term $rhs$ in (\ref{3dinvert}) is the divergence of the vector $\vec{\mathbf{F}}.$ By @@ -598,7 +603,7 @@ {inhomneumann}) the modified boundary condition becomes: \begin{equation} -\widehat{n}.\nabla \phi _{nh}=0 \label{homneuman} +\widehat{n}.\nabla \phi _{nh}=0 \label{eq:hom-neumann-nh} \end{equation} If the flow is `close' to hydrostatic balance then the 3-d inversion @@ -622,10 +627,11 @@ Many forms of momentum dissipation are available in the model. Laplacian and biharmonic frictions are commonly used: -\[ -D_{V}=A_{h}\nabla _{h}^{2}v+A_{v}\frac{\partial ^{2}v}{\partial z^{2}}% +\begin{equation} +D_{V}=A_{h}\nabla _{h}^{2}v+A_{v}\frac{\partial ^{2}v}{\partial z^{2}} +A_{4}\nabla _{h}^{4}v -\] +\label{eq:dissipation} +\end{equation} where $A_{h}$ and $A_{v}\ $are (constant) horizontal and vertical viscosity coefficients and $A_{4}\ $is the horizontal coefficient for biharmonic friction. These coefficients are the same for all velocity components. @@ -634,18 +640,19 @@ The mixing terms for the temperature and salinity equations have a similar form to that of momentum except that the diffusion tensor can be -non-diagonal and have varying coefficients. $\qquad $% -\[ +non-diagonal and have varying coefficients. $\qquad $ +\begin{equation} D_{T,S}=\nabla .[\underline{\underline{K}}\nabla (T,S)]+K_{4}\nabla _{h}^{4}(T,S) -\] -where $\underline{\underline{K}}\ $is the diffusion tensor and the $K_{4}\ $% +\label{eq:diffusion} +\end{equation} +where $\underline{\underline{K}}\ $is the diffusion tensor and the $K_{4}\ $ horizontal coefficient for biharmonic diffusion. In the simplest case where the subgrid-scale fluxes of heat and salt are parameterized with constant horizontal and vertical diffusion coefficients, $\underline{\underline{K}}$, reduces to a diagonal matrix with constant coefficients: -\[ +\begin{equation} \qquad \qquad \qquad \qquad K=\left( \begin{array}{ccc} K_{h} & 0 & 0 \\ @@ -653,7 +660,8 @@ 0 & 0 & K_{v} \end{array} \right) \qquad \qquad \qquad -\] +\label{eq:diagonal-diffusion-tensor} +\end{equation} where $K_{h}\ $and $K_{v}\ $are the horizontal and vertical diffusion coefficients. These coefficients are the same for all tracers (temperature, salinity ... ). @@ -664,10 +672,10 @@ {hor-mtm}) and (\ref{vertmtm}) in the (so-called) `vector invariant' form: \begin{equation} -\frac{D\vec{\mathbf{v}}}{Dt}=\frac{\partial \vec{\mathbf{v}}}{\partial t}% +\frac{D\vec{\mathbf{v}}}{Dt}=\frac{\partial \vec{\mathbf{v}}}{\partial t} +\left( \nabla \times \vec{\mathbf{v}}\right) \times \vec{\mathbf{v}}+\nabla \left[ \frac{1}{2}(\vec{\mathbf{v}}\cdot \vec{\mathbf{v}})\right] -\label{vecinvariant} +\label{eq:vi-identity} \end{equation} This permits alternative numerical treatments of the non-linear terms based on their representation as a vorticity flux. Because gradients of coordinate