update for Xylars comments

Author: mark-petersen

components/omega/doc/design/OmegaV1GoverningEqns.md

@@ -200,9 +200,27 @@ Omega-1 will be non-Boussinesq, which means that density is not assumed to be co
 
 For non-Boussinesq layered equations we begin by defining the pressure-thickness of layer $k$ as
 
-$$ (sigma-int-to-surf)
+$$
+h_k(x,y,t)
+ \equiv \int_{s_k^{bot}}^{s_k^{top}} \rho g dz = \int_{\sigma_k^{top}}^{\sigma_k^{bot}} dp.
+$$ (def-h)
+
+The letter $h$ is used because this is the familiar variable for thickness in Boussinesq primitive equations. In the Boussinesq case $h$ is volume divided by area, so that the thickness equation describes conservation of volume. Here for the non-Boussinesq case, $h_k$ is mass per unit area, multiplied by $g$. Since horizontal cell area and $g$ remain constant, the thickness equation describes conservation of mass. 
+
+Throughout this derivation we can write all equations equivalently in $z$-coordinates (depth), or in $p$-coordinates (pressure). From the hydrostatic equation, any quantity $\varphi$ may be integrated in $z$ or $p$ as
 
-where $\eta$ is the sea surface height and $p^{surf}$ is the surface pressure. Note that the negative sign means that pressure increases with depth, so positive $p$ points downward, so that the $top$ and $bot$ extents of the integration limits are flipped in [](#depth-pressure-integral-conversion).
+$$
+\int_{s_k^{bot}}^{s_k^{top}} \varphi \rho g dz
+= \int_{\sigma_k^{bot}}^{\sigma_k^{top}} \varphi dp
+$$(depth-pressure-integral-conversion)
+
+where we can convert from an interfacial depth surface to a pressure surface with
+
+$$
+\sigma(x,y) = - \int_{s(x,y)}^{\eta(x,y)} \rho g dz + p^{surf}(x,y)
+$$ (def-sigma)
+
+where $\eta$ is the sea surface height and $p^{surf}$ is the surface pressure. Note that the negative sign means that pressure increases with depth.  This means that positive $p$ points downward, so that the $top$ and $bot$ extents of the integration limits are flipped in [](#depth-pressure-integral-conversion).
 
 For any three-dimensional quantity $\varphi(x,y,z,t)$, the pressure-thickness-averaged quantity in layer $k$ is
 
@@ -218,7 +236,8 @@ Our derivation has not made any assumptions about density, and may be used for b
 $$
 h_k^{Bouss}(x,y,t)
  \equiv \int_{s_k^{bot}}^{s_k^{top}} dz.
-$$
+$$ (def-h-bouss)
+
 The definition of the thickness-weighted average remains the same as in [](#def-thickness-pressure-average), but in a Boussinesq fluid the variations in density are assumed to be small, so that $\rho(x,y,z,t)\sim\rho_0$ where $\rho_0$ is a constant. 
 
 
@@ -290,12 +309,12 @@ $$
 + q_k\left(h\mathbf{u}_k^{\perp}\right) 
 + \nabla K_k
 = 
-- \frac{1}{\rho} \nabla_r p - \nabla_r \phi
+- \frac{1}{\rho} \nabla_r p - \nabla_r \Phi
 + {\bf D}_k^u + {\bf F}_k^u
 $$ (layered-momentum)
 
 $$
-\frac{\partial h_k}{\partial t} + \nabla_r \cdot \left(h_k \mathbf{u}_k\right) + \omega_{bot} - \omega_{top}= 0,
+\frac{\partial h_k}{\partial t} + \nabla_r \cdot \left(h_k \mathbf{u}_k\right) + \omega_{bot} - \omega_{top}= S^h_k
 $$ (layered-mass)
 
 $$
@@ -310,15 +329,17 @@ The discretized versions of the governing equations are
 
 $$
 \frac{\partial u_e}{\partial t} + \left[ \frac{{\mathbf k} \cdot \nabla \times u_e +f_v}{[h_i]_v}\right]_e\left([h_i]_e u_e^{\perp}\right)
-= -g\nabla(h_i-b_i) - \nabla K_i + \nu_2 \nabla^2 u_e - \nu_4 \nabla^4 u_e + \mathcal{D}_e + \mathcal{F}_e
+= 
+- \frac{1}{\rho_i} \nabla_r p_i - \nabla_r \Phi_i
+- \nabla K_i + \nu_2 \nabla^2 u_e - \nu_4 \nabla^4 u_e + \mathcal{D}_e + \mathcal{F}_e
 $$ (discrete-momentum)
 
 $$
 \frac{\partial h_i}{\partial t} + \nabla \cdot \left([h_i]_e u_e\right) = 0,
 $$ (discrete-thickness)
 
 $$
-\frac{\partial h_i \phi_i}{\partial t} + \nabla \cdot \left(u_e [h_i \phi_i]_e \right) = \kappa_2 h_i \nabla^2 \phi_i - \kappa_4 h_i \nabla^4 \phi_i,
+\frac{\partial h_i \varphi_i}{\partial t} + \nabla \cdot \left(u_e [h_i \varphi_i]_e \right) = \kappa_2 h_i \nabla^2 \varphi_i - \kappa_4 h_i \nabla^4 \varphi_i,
 $$ (discrete-tracer)
 
 where subscripts $i$, $e$, and $v$ indicate cell, edge, and vertex locations ($i$ was chosen for cell because $c$ and $e$ look similar).  Here square brackets $[\cdot]_e$ and $[\cdot]_v$ are quantities that are interpolated to edge and vertex locations. The interpolation is typically centered, but may vary by method, particularly for advection schemes. For vector quantities, $u_e$ denotes the normal component at the center of the edge, while $u_e^\perp$ denotes the tangential component. In the discrete system, the normal component $u_e$ points positively from the lower cell index to the higher cell index, while the tangential component $u_e^\perp$ points positively $90^o$ to the left of $u_e$ (for unit vectors, ${\mathbf n}_e^\perp = {\mathbf k}\times {\mathbf n}_e$).
@@ -372,7 +393,7 @@ Table 1. Definition of variables
 | $\kappa_4$          | biharmonic tracer diffusion | m$^4$/s    | cell     |   |                                                              |
 | $\nu_2$             | viscosity                   | m$^2$/s    | edge     |   |                                                              |
 | $\nu_4$             | biharmonic viscosity        | m$^4$/s    | edge     |   |                                                              |
-| $\phi$              | tracer                      | varies | cell     |   | units may be kg/m$^3$ or similar | 
+| $\varphi$              | tracer                      | varies | cell     |   | units may be kg/m$^3$ or similar | 
 | $\omega$             | relative vorticity          | 1/s      | vertex   |  RelativeVorticity | $\omega={\mathbf k} \cdot \left( \nabla \times {\mathbf u}\right)$ |
 |$w$ | vertical transport | | cell  | determined by coord. type |
 |$z$ | vertical coordinate | - | positive upward |
@@ -385,9 +406,10 @@ Table 1. Definition of variables
 |$\rho$ | density | | cell  |
 |$\rho_0$ | reference density | | constant | |
 |$\varphi$ | tracer | | cell | e.g. $\Theta$, $S$ |
+|$\Phi$ | geopotential| | cell | mrp todo define here |
 
-## Operator Formulation
-See
+### Operator Formulation
+The horizontal operator stencils remain the same as those given in [Omega-0 design document](OmegaV0ShallowWater.md#operator-formulation).
 
 ## Verification and Testing
 
@@ -414,7 +436,7 @@ See [Petersen et al. 2015](http://www.sciencedirect.com/science/article/pii/S146
 This is a 3D domain with a seamount in the center, where temperature and salinity are stratified in the vertical and constant in the horizontal. The test is simply that an initial velocity field of zero remains zero. For z-level layers the velocity trivially remains zero because the horizontal pressure gradient is zero. For tilted layers, this is a test of the pressure gradient error and the velocity is never exactly zero. This is a common test for sigma-coordinate models like ROMS because the bottom layers are extremely tilted along the seamount, but it is a good test for any model with tilted layers. Omega will use slightly tilted layers in p-star mode (pressure layers oscillating with SSH) and severely tilted layers below ice shelves, just like MPAS-Ocean. See [Ezer et al. 2002](https://www.sciencedirect.com/science/article/pii/S1463500302000033), [Haidvogel et al. 1993](https://journals.ametsoc.org/view/journals/phoc/23/11/1520-0485_1993_023_2373_nsofaa_2_0_co_2.xml), [Shchepetkin and McWilliams 2003](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2001JC001047), and previous MPAS-Ocean [confluence page](https://acme-climate.atlassian.net/wiki/spaces/OCNICE/blog/2015/11/19/40501447/MPAS-O+Sigma+coordinate+test+sea+mount).
 
 ### Cosine Bell on the Sphere
-This test uses a fixed horizontal velocity field to test horizontal tracer advection. It is repeated from Omega-0 and is important to conduct again as we convert Omega to a layered primitive-equation model. See `cosine_bell` case in both compass and polaris.
+This test uses a fixed horizontal velocity field to test horizontal tracer advection. It is repeated from [Omega-0 design document](OmegaV0ShallowWater) and is important to conduct again as we convert Omega to a layered primitive-equation model. See `cosine_bell` case in both compass and polaris.
 
 ### Merry-Go-Round
 This is an exact test for horizontal and vertical tracer advection. A fixed velocity field is provided, and a tracer distribution is advected around a vertical plane. See the `merry_go_round` test in compass, and the results on the [merry-go-round pull request](https://github.com/MPAS-Dev/compass/pull/108) and [compass port pull request](https://github.com/MPAS-Dev/compass/pull/452).