expand continuous equations

Author: mark-petersen

components/omega/doc/design/OmegaV1GoverningEqns.md

@@ -11,7 +11,7 @@ Add later, if it seems necessary. There is a toc on the left bar.
 -->
 
 
-## Overview
+## 1. Overview
 
 This design document describes the governing equations for Omega, the Ocean Model for E3SM Global Applications. Overall, Omega is an unstructured-mesh ocean model based on TRiSK numerical methods ([Thuburn et al. 2009](https://www.sciencedirect.com/science/article/pii/S0021999109004434)) that is specifically designed for modern exascale computing architectures. The algorithms in Omega will be mostly identical to those in MPAS-Ocean, but it will be written in c++ rather than Fortran in order to take advantage of libraries to run on GPUs, such as Kokkos ([Trott et al. 2022](https://ieeexplore.ieee.org/document/9485033)). Significant differences between MPAS-Ocean and Omega are:
 
@@ -24,28 +24,28 @@ The planned versions of Omega are:
 - **Omega-1.0: Layered ocean, idealized, no surface fluxes.** This adds active temperature, salinity, and density as a function of pressure in the vertical. Vertical advection and diffusion terms are added to the momentum and tracer equations. An equation of state and simple vertical mixing, such as constant coefficient, are needed. Capability and testing are similar to [Petersen et al. 2015](http://www.sciencedirect.com/science/article/pii/S1463500314001796). Tests include overflow, internal gravity wave, baroclinic channel, seamount, and vertical merry-go-round.
 - **Omega-1.1: Layered ocean, idealized, with surface fluxes.** Addition of simple vertical mixing scheme such as Pacanowski & Philander; nonlinear equation of state (TEOS10); tracer surface restoring to a constant field; constant wind forcing; and flux-corrected transport for horizontal advection. Testing will be with the baroclinic gyre and single column tests of surface fluxes and vertical mixing. 
 - **Omega-2.0: Coupled within E3SM, ocean only.** Ability to run C cases (active ocean only) within E3SM. Requires addition of E3SM coupling infrastructure; simple analysis (time-averaged output of mean, min, max); split baroclinic-barotropic time; global bounds checking on state. Testing and analysis similar to [Ringler et al. 2013](https://www.sciencedirect.com/science/article/pii/S1463500313000760), except Omega uses a non-Boussinesq formulation.
-- **Omega-2.1: E3SM fully coupled** Ability to run G cases (active ocean and sea ice) and B cases (all components active) within E3SM. This will include: a full vertical mixing scheme, such as KPP; frazil ice formation in the ocean; and a submosescale parameterization. Parameterizations mean for lower resolutions, such as Gent-McWilliams and Redi mixing, will not be in this version. Simulations can be run at high resolution and compared to previous E3SM simulations, including [Petersen et al. 2019](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018MS001373) and [Caldwell et al. 2019](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001870).
+- **Omega-2.1: E3SM fully coupled** Ability to run G cases (active ocean and sea ice) and B cases (all components active) within E3SM. This will include: a full vertical mixing scheme, such as KPP; frazil ice formation in the ocean; and a submosescale parameterization. Omega 2.1 will mostly be run at eddy-resolving resolutions, and will not include parameterizations for lower resolutions such as Gent-McWilliams and Redi mixing.  Simulations will be compared to previous E3SM simulations, including [Caldwell et al. 2019](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001870) and [Petersen et al. 2019](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018MS001373).
 
 This document describes the governing equations for the layered ocean model, which are applicable for Omega-1.0 onwards. Specific terms, such as the pressure gradient, vertical mixing, and parameterizations, are left in their general form here but are described in more detail in other design documents.
 
-## Requirements
+## 2. Requirements
 
 The requirements in the [Omega-0 design document](OmegaV0ShallowWater) still apply. Additional design requirements for Omega-1 are:
 
 ### Omega will be an incompressible, hydrostatic, non-Boussinesq ocean model.
 See discussion in introduction. Non-Boussinesq changes the independant vertical variable to pressure, and mass conservation equation to be pressure-thickness, $h=\Delta p$.
 
 ### Omega will use TEOS10 for the equation of state.
-See additional EOS design document.
+See additional [EOS design document](EOS)
 
 ### Omega-1.0 will add new terms for the pressure gradient, vertical mixing, and vertical advection.
 See individual design documents for each term.
 
-## Governing Equations
+## 3. Governing Equations
 
 ### MPAS-Ocean Equations of Motion
 
-The MPAS-Ocean layered formulation solves the following equations for momentum, thickness, and tracers at layer $k$. These are continuous in the horizontal and discrete in the vertical.
+The MPAS-Ocean layered formulation are provided here for reference. MPAS-Ocean solves for momentum, thickness, and tracers at layer $k$. These are continuous in the horizontal and discrete in the vertical.
 
 $$
 \frac{\partial {\mathbf u}_k}{\partial t}
@@ -74,17 +74,15 @@ $$ (mpaso-continuous-tracer)
 The layer thickness $h$, vertical velocity $w$, pressure $p$, and tracer $\varphi$, are cell-centered quantities, while the horizontal velocity ${\mathbf u}$ and $e$ superscript are variables interpolated to the cell edges.
 
 ### Continuous Equations
-The continuous form of the conservation equations are as follows [Kundu et al. 2024](https://www.amazon.com/Fluid-Mechanics-Pijush-K-Kundu/dp/012405935X), chapter 4, eqns 4.7 and 4.22. This is before any assumptions are made, so this is a compressible, non-hydrostatic, non-Boussinesq fluid. Here all variables are a function of $(x,y,z)$, $\bf u$ denotes the horizontal velocity vector, $w$ is the vertical velocity, ${\bf u} \otimes {\bf u} = {\bf u}{\bf u}^T$ is the tensor product, $\nabla$ is the horizontal gradient, and other variables defined in the [section below](#variable-definitions).
+The continuous form of the conservation equations are as follows [Kundu et al. 2024](https://www.amazon.com/Fluid-Mechanics-Pijush-K-Kundu/dp/012405935X), chapter 4, eqns 4.7 and 4.22. This is before any assumptions are made, so this is a compressible, non-hydrostatic, non-Boussinesq fluid. Here all variables are a function of $(x,y,z)$, ${\bf u}_{3D}$ denotes the three-dimensional velocity vector, ${\bf u}_{3D} \otimes {\bf u}_{3D} = {\bf u}_{3D}{\bf u}_{3D}^T$ is the tensor product, $\nabla_{3D}$ is the three-dimensional gradient, $D/Dt$ is the material derivative, and other variables defined in the [Variable Defitionion Section below](#variable-definitions).
 
 momentum:
 
 $$
-\frac{D \rho {\bf u} }{D t} \equiv
-\frac{\partial \rho {\bf u}}{\partial t}
- + \nabla \cdot \left( \rho {\bf u} \otimes {\bf u}  \right)
-+ w\frac{\partial \rho {\bf u}}{\partial z}
-  = - \nabla p 
-- \rho f {\bf k} \times {\bf u}
+\frac{D \rho {\bf u}_{3D} }{D t} \equiv
+\frac{\partial \rho {\bf u}_{3D}}{\partial t}
+ + \nabla_{3D} \cdot \left( \rho {\bf u}_{3D} \otimes {\bf u}_{3D}  \right)
+  = - \nabla_{3D} p 
 + \rho {\bf D}^u + \rho {\bf F}^u 
 $$ (continuous-momentum-R13App)
 
@@ -93,8 +91,7 @@ mass:
 $$
 \frac{D \rho}{D t} \equiv
 \frac{\partial \rho }{\partial t}
- + \nabla \cdot \left( \rho  {\bf u}  \right)
- + \frac{\partial}{\partial z} \left(  \rho  w \right)
+ + \nabla_{3D} \cdot \left( \rho  {\bf u}_{3D}  \right)
 = 0
 $$ (continuous-mass-R13App)
 
@@ -103,24 +100,23 @@ tracer:
 $$
 \frac{D \rho \varphi }{D t} \equiv
 \frac{\partial \rho \varphi}{\partial t}
- + \nabla \cdot \left( \rho \varphi {\bf u}  \right)
- + \frac{\partial}{\partial z} \left(  \rho \varphi w \right)
+ + \nabla_{3D} \cdot \left( \rho \varphi {\bf u}_{3D}  \right)
 = D^\varphi + S^\varphi
 $$ (continuous-tracer-R13App)
 
 Here we have express the following terms as a general operators, with examples of specific forms provided below: the dissipation ${\bf D}^u$, momemtum forcing ${\bf F}^u$, tracer diffusion $D^\varphi$, and tracer sources and sinks $S^\varphi$. The momentum equation can be rewritten using the product rule on $\rho {\bf u}$, mass conservation, and dividing by $\rho$, as:
 
 $$
-\frac{D {\bf u} }{D t} \equiv
-\frac{\partial {\bf u}}{\partial t}
- + \nabla \cdot \left( {\bf u} \otimes {\bf u}  \right)
-+ w\frac{\partial {\bf u}}{\partial z}
-  = - \frac{1}{\rho} \nabla p 
-- f {\bf k} \times {\bf u}
-+ {\bf D}^u + {\bf F}^u
+\frac{D {\bf u}_{3D} }{D t} \equiv
+\frac{\partial {\bf u}_{3D}}{\partial t}
+ + \nabla_{3D} \cdot \left( {\bf u}_{3D} \otimes {\bf u}_{3D}  \right)
+  = - \frac{1}{\rho} \nabla_{3D} p 
++ {\bf D}^u + {\bf F}^u 
 $$ (continuous-momentum-rho)
 
-Note that the advective term may be replaced with the vector identity
+Geophysical fluids such as the ocean and atmosphere are rotating and stratified, and horizontal velocities are orders of magnitude larger than vertical velocities. It is therefore convenient to separate the horizontal and vertical as ${\bf u}_{3D} = \left( {\bf u}, w \right)$ and $\nabla_{3D} = \left( \nabla, d/dz \right)$ where $z$ is the vertical direction in a local coordinate frame on the earth's surface and $w$ is the vertical velocity. The gravitational force $g {\bf k}$ is added in the vertical. The rotating frame of reference results in the Coriolis force $f {\bf k} \times {\bf u} \equiv f {\bf u}^\perp$, where $f$ is the Coriolis parameter and ${\bf u}^\perp$ is the horizontal velocity rotated $90^o$ from $\bf u$. See any textbook in the [References](#references) for a full derivation.
+
+The horizontal parts of the advective term may be replaced with the vector identity
 
 $$
 \nabla \cdot \left( {\bf u} \otimes {\bf u}  \right) = 
@@ -129,7 +125,7 @@ $$
 &= \zeta \boldsymbol{u}^{\perp} + \nabla K.
 $$ (advection-identity)
 
-where $\zeta$ is relative vorticity and $K$ is kinetic energy. This is used to produce the final formulation below.
+where $\zeta$ is relative vorticity and $K$ is kinetic energy. This step separates the horizontal advection into non-divergent and non-rotational components, which is useful in the final TRiSK formulation.
 
 For a primitive equation ocean model, we assume the fluid is incompressible and hydrostatic. For Omega we are not making the Boussinesq assumption, so all density terms are the full density. In particular, the density coefficient of the pressure gradient in [](#continuous-momentum-R13App) is not a constant, as it is in primitive equation models like POP and MPAS-Ocean.
 
@@ -147,7 +143,7 @@ $$
   - \frac{1}{\rho} \frac{\partial p}{\partial z} - g + {\bf D}^u + {\bf F}^u
 $$ (continuous-vert-mom)
 
-assume that advection of vertical momentum $Dw/Dt$, dissipation, and forcing are small, and that the first order balance is between pressure gradient and gravity,
+assume that advection of vertical momentum $Dw/Dt$, dissipation, and forcing are small, and that the first order balance is between pressure gradient and buoyancy,
 
 $$
 \frac{\partial p}{\partial z}
@@ -164,35 +160,56 @@ The constitutive equation is the equation of state,
 $$
 \rho = f_{eos}(p,T,S).
 $$ (continuous-eos-R13App)
-where are temperature and salinity ($T$ and $S$) are examples of tracers $\varphi$.
+where temperature and salinity ($T$ and $S$) are examples of tracers $\varphi$.
 
-Momentum diffusion terms include
+The final form of the continuous conservations equations for a non-Boussinesq, incompressible, non-hydrostatic ocean are
+
+momentum:
 
 $$
-{\bf D}^u =  {\bf D}^u_h + {\bf D}^u_v
-$$ (continuous-mom_diff-R13App)
+\frac{\partial {\bf u}}{\partial t}
++ f {\bf u}^\perp + \zeta {\bf u}^\perp
++ w\frac{\partial {\bf u}}{\partial z}
+  = - \frac{1}{\rho} \nabla p 
+  - \nabla K
++ {\bf D}^u + {\bf F}^u 
+$$ (continuous-momentum-final)
+
+mass:
+
 $$
- {\bf D}^u_v=\frac{\partial }{\partial z}
-\left( \nu_v \frac{\partial {\bf u}}{\partial z} \right)
-$$ (continuous-v_mom_diff-R13App)
+\frac{\partial \rho }{\partial t}
+ + \nabla \cdot \left( \rho  {\bf u}  \right)
+ + \frac{\partial}{\partial z} \left(  \rho  w \right)
+= 0
+$$ (continuous-mass-final)
+
+tracer:
+
 $$
- {\bf D}^u_h=\nu_h \nabla^2 \bf u = \nu_h(\nabla \delta + {\bf k}\times \nabla \eta)
-$$ (continuous-h_mom_diff-R13App)
-and may also include a del-4 term, Rayleigh drag, and bottom drag. Momentum forcing is due to wind stress.
+\frac{\partial \rho \varphi}{\partial t}
+ + \nabla \cdot \left( \rho \varphi {\bf u}  \right)
+ + \frac{\partial}{\partial z} \left(  \rho \varphi w \right)
+= D^\varphi + S^\varphi
+$$ (continuous-tracer-final)
 
-Tracer diffusion terms include
+Momentum diffusion terms include Laplacian (del2), biharmonic (del4), and vertical viscosity, 
 
 $$
-D^\varphi =  D^\varphi_h + D^\varphi_v
-$$ (continuous-tracer_diff-R13App)
+{\bf D}^u 
+= \nu_2 \nabla^2 {\bf u} - \nu_4 \nabla^4 {\bf u}
+ + \frac{\partial }{\partial z} \left( \nu_v \frac{\partial {\bf u}}{\partial z} \right)
+$$ (continuous-h_mom_diff)
+and may also include a Rayleigh drag and eventually parameterizations. Momentum forcing is due to wind stress and bottom drag. Similarly, the tracer diffusion terms include Laplacian (del2), and vertical viscosity,
+
 $$
- D^\varphi_v = \rho \frac{\partial }{\partial z}
+D^\varphi =   
+ \nabla\cdot\left(\rho \kappa_2 \nabla\varphi \right)
++ \rho \frac{\partial }{\partial z}
   \left( \kappa_v \frac{\partial \varphi}{\partial z} \right),
 $$ (continuous-v_tr_diff-R13App)
-$$
- D^\varphi_h = \nabla\cdot\left(\rho \kappa_h \nabla\varphi \right).
-$$ (continuous-h_tr_diff-R13App)
-and may also include a del-4 term. Sources and sinks include surface fluxes from the atmosphere and land, and bio-geo-chemical reactions.
+and may also include a biharmonic (del4) term and parameterizations such as Redi mixing. Sources and sinks include surface fluxes from the atmosphere and land, and bio-geo-chemical reactions.
+All of the diffusion and forcing terms are written in more detail with the [Discrete Equations](#discrete-equations) below.
 
 ### Layered Equations
 
@@ -365,11 +382,10 @@ Table 1. Definition of variables
 | $C_W$               | wind stress coefficient     | 1/m      | constant |   |                         |
 | $D$             | divergence | 1/s      | cell | Divergence  | $D=\nabla\cdot\mathbf u$ |
 | $\mathcal{D} $      | drag                        | m/s$^2$    | edge     |   |                       |
-| $f$                 | Coriolis parameter          | 1/s      | vertex   | FVertex  |                  |
-| ${\bf F}^u $      | momentum forcing                     | m/s$^2$    | edge     |   |                       |
 |${\bf D}^u $, ${\bf D}^u_h$, ${\bf D}^u_v$  | momentum diffusion terms | | edge | $h$ horizonal, $v$ vertical|
 |$ D^\varphi$, $ D^\varphi_h$, $ D^\varphi_v$ | tracer diff. terms | | cell | |
-|$f$ | Coriolis parameter | | vertex |
+| $f$                 | Coriolis parameter          | 1/s      | vertex   | FVertex  |  $f = 2\Omega sin(\phi)$, $\Omega$ rotation rate, $\phi$ latitude|
+| ${\bf F}^u $      | momentum forcing                     | m/s$^2$    | edge     |   |                       |
 |$f_{eos}$ | equation of state | -  |
 |$F$ | thickness flux | | edge | $F=hu$ |
 |$g$ | grav. acceleration | | constant  |
@@ -390,7 +406,7 @@ Table 1. Definition of variables
 |$T^u,T^h,T^\varphi$ | tendencies | - |  |
 | ${\mathbf u}$   | velocity, vector form       | m/s      | edge     |   |                                                              |k
 | $u_{e,k}$   | velocity, normal to edge      | m/s      | edge     | NormalVelocity  |                                                              |k
-| $u^\perp_{e,k}$   | velocity, tangential to edge      | m/s      | edge     | TangentialVelocity  |                                                              |k
+| $u^\perp_{e,k}$   | velocity, tangential to edge      | m/s      | edge     | TangentialVelocity  | $u^\perp = {\bf k} \times {\bf u}|
 | ${\mathbf u}_W$ | wind velocity               | m/s      | edge     |   |                                                              |
 | $\eta$              | absolute vorticity          | 1/s      | vertex   |   | $\eta=\omega + f$ |
 | $\kappa_2$          | tracer diffusion            | m$^2$/s    | cell     |   |                                                              |
@@ -415,7 +431,7 @@ Table 1. Definition of variables
 ### Operator Formulation
 The horizontal operator stencils remain the same as those given in [Omega-0 design document](OmegaV0ShallowWater.md#operator-formulation).
 
-## Verification and Testing
+## 4. Verification and Testing
 
 Capability and testing are similar to [Petersen et al. 2015](http://www.sciencedirect.com/science/article/pii/S1463500314001796). The following tests are in idealized domains and do not require surface fluxes or surface restoring. For the following tests to show results comparable to those published with other models, the full dynamic sequence of density, pressure, momentum, and advection must work correctly. The successful completion of the following tests is a validation of the primitive equation functions in Omega 1.0. All of the following tests may exercise a linear equation of state or the nonlinear TEOS10. The first four tests quantify the anomalous mixing caused by the numerical schemes. The first five are on cartesian planes with regular hexagon meshes.
 

components/omega/doc/index.md

@@ -88,7 +88,6 @@ devGuide/TridiagonalSolvers
 
 ```{toctree}
 :caption: Design documents
-:numbered:
 :maxdepth: 1
 
 design/OmegaV0ShallowWater