Omega V1 Governing Equations Design Document

[mark-petersen]

Derivation of non-Boussinesq primitive equations for Omega Version 1.

See compiled document at
https://portal.nersc.gov/project/e3sm/lvroekel/omega_docs/html/design/OmegaV1GoverningEqns.html

[xylar]

Some comments, mostly about generalizing the gravitational potential and adding it a few places I think it's missing.

Also, a request to remove the commented-out text to make review easier.

[mark-petersen]

I have added substantial updates based on feedback and conversations this past week. In particular, the document now uses layer mass-thickness $h$ in kg/m^2. Squashing commits, rebasing and pushing a single commit.

[vanroekel]

the buoyancy term is missing here rho g khat

[mark-petersen]

The first equation set is the most general, so gravity would be within ${\bf F}^u$ but is not yet specified. A few paragraphs down I begin "Geophysical fluids such as the ocean and atmosphere are rotating and stratified..." and start to make assumptions and call out specific terms. I have the sentence "The gravitational force $g {\bf k}$ is included in ${\bf F}^u$".

[vanroekel]

oh sorry, totally missed that. Thanks for the explanation

[xylar]

I realize that one could include gravity in ${\bf F}^u$ but I would much prefer that we not do so. In my view, it makes a lot of the subsequent derivation clearer if we include $\nabla_{3D} \Phi$ here from the start. I have been pushing for this the whole time and will continue to do so.

[vanroekel]

All, I've made a pretty substantive change in the subgrid scale parameterization section that addresses how to do vertical derivatives in a finite volume framework and recasting the temperature vertical diffusion discussion to be in line with the tri-diagonal solve.

The main big question in my mind is to double/triple check the pressure gradient derivation as it is the main difference between this doc and MPAS-Ocean.

[xylar]

I didn't get a chance to do a last read-through today as I had hoped. But I'm happy to approve things based on the hard work of @mark-petersen, @vanroekel, @sbrus89, @katsmith133, @cbegeman, @alicebarthel and others, and my glancing at things on Wednesday (the pressure-gradient terms in particular).

[hyungyukang]

Suggested change

	Rearranging, is it useful to note that
	Rearranging, it is useful to note that

[katsmith133]

Most of my comments are to fix typos, broken links, and formatting. The only major thing is that I was not able to follow the Vertical Tracer Diffusion section. I followed the vertical derivatives part and the Vertical Momentum Dissipation section just fine, but then was confused on how and why things were different in the Vertical Tracer Diffusion section. It was a bit confusing and lacked some explanations. Otherwise, I think the rest read and flowed well.

[katsmith133]

"In an equation for $\tilde{h}_k \varphi$, this term is different than in the velocity equation for vertical turbulent diffusion."

I think we should explain why its different and we don't just treat it the same way.

[katsmith133]

Define the Adv and Hdiff terms?

[katsmith133]

This section in general is a bit confusing and could benefit from some rewriting or more explanation.

[hyungyukang]

I only have a few minor comments. The document was generally easy to follow. I agree with @katsmith133 that the vertical tracer diffusion section could be clarified.

[hyungyukang]

I tried to improve the vertical tracer diffusion section. @vanroekel, please review this first and feel free to revise.

[vanroekel]

I've made a few updates to W_tr to try be a bit more clear. Since perlmutter is down the docs are here -

https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.vanroekel/omega_doc/html/design/OmegaV1GoverningEqns.html

but all my changes are pushed to this PR.

[mark-petersen]

I am concerned about the pressure gradient term, e.g. in (50)

The pressure gradient should be negative. That is, flow is induced from high to low pressure.

The pseudothickness inside the gradient does not look physical. If you have uniform horizontal pressure, uniform density, and level layers that are thinning horizontally, this term would induce velocity and it should not. Perhaps that is cancelled by the grad z terms that follow?

What if we apply the product rule so that

grad( h alpha p) = alpha p grad (h) + h grad(alpha p)

and gather the first term with the grad(z) terms? Then the three terms have similar form, just evaluated at (top, mid, bot).

[vanroekel]

I can definitely fix the negative, that's easy and correct to do.

I'm not sure I like the alternative proposed. I think this form is physical, just the layer integrated pressure gradient. I'd agree it looks different than most expect, but I think it's physically correct. I don't have a super strong feeling about the alternate form. I can see how that would be cleaner to understand for people.

[mark-petersen]

If you write h_k = z^{top} - z^{bot} then the grad( h_k ) combines with the grad( z ) terms and only differs by the evaluation of alpha p at top vs layer avg and bot vs layer avg. That's interesting.

[vanroekel]

that is very interesting. if we use (20) can we write the first term (alpha p \grad h) as

alpha p grad (ztop-zbot)

?
this looks even closer to the bottom line terms but still not quite there I think due to alpha and p being layer averages

[cbegeman]

Suggested change

	The second is the surface force ${\bf f}$, which acts on the boundary of the control volume and includes pressure and viscous stresses. These forces appear as surface integrals over the boundary and drive momentum exchange between adjacent fluid parcels or between the fluid and its environment.
	The second is the surface force ${\bf \tau}$, which acts on the boundary of the control volume and includes pressure and viscous stresses. These forces appear as surface integrals over the boundary and drive momentum exchange between adjacent fluid parcels or between the fluid and its environment.

[cbegeman]

It isn't clear why the notation $z'$ is used. Also see comment below.

[cbegeman]

Is there a point later on where this is helpful? If so, can we cite an equation?

[katsmith133]

Added a schematic for the pseudo-velocities: https://portal.nersc.gov/project/e3sm/katsmith/omega_docs/html/design/OmegaV1GoverningEqns.html

Feel free to change wording. I could not for the life of me figure out how to label and reference it. I tried:
```{figure} images/tilde_w.jpeg
:label: tilde-w-figure
:align: center
:width: 300 px

Schematic of pseudo-velocities.
```
and referencing with [](#tilde-w-figure), but that did not seem to work.

[cbegeman]

The work here is really impressive! It's great to have this level of detail in the derivation. My comments are pretty minor and mostly pertain to clarity. Thanks all, esp @vanroekel for the latest edits.

[cbegeman]

Suggested change

	\tilde{z}(z) = -\frac{1}{\rho_0 g} \left( p^{\text{surf}} + \int_{z}^{z^{\text{surf}}} \rho(z') g \, dz' \right)
	\tilde{z}(x, y, z, t) = -\frac{1}{\rho_0 g} \left( p^{\text{surf}(x, y, z, t)} + \int_{z}^{z^{\text{surf}}} \rho(x, y, z', t) g \, dz' \right)

I think this aids clarity even though it is perhaps overkill.

[cbegeman]

Suggested change

	- Vertical transport of mass, tracers, and momentum will use $\tilde{w}$ rather than volume flux $w$.
	- Vertical transport of mass, tracers, and momentum will use $\tilde{w}$ rather than the vertical velocity in height coordinates $w$.

It's weird to call w a volume flux since it has units of m/s.

[cbegeman]

Suggested change

	This approach ensures consistency with mass conservation and simplifies vertical discretization.
	This approach ensures consistency with mass conservation and simplifies derivation of the discrete equations.

Is this what you mean? If not, can we say why this simplifies vertical discretization? It isn't obvious to me.

[cbegeman]

Is $\rho(x, y, z, t)$ different from $\rho(x, y, z', t)$?

[cbegeman]

Suggested change

	We partition the control surface of a finite-volume cell into the side walls $\partial V^{\text{side}}$ (which are fixed in space) and the upper and lower pseudo-height surfaces $\partial V^{\text{top}}(t)$ and $\partial V^{\text{bot}}(t)$, which evolve in time.
	We partition the control surface of a finite-volume cell into the side walls $\partial V^{\text{side}}$ (which are fixed in time) and the upper and lower pseudo-height surfaces $\partial V^{\text{top}}(t)$ and $\partial V^{\text{bot}}(t)$, which evolve in time.

Is this a typo or am I misunderstanding what you mean?

[cbegeman]

Suggested change

	This form can be interfaced with the Omega [tridiagongal solver](TridiagonalSolver.md) routine.
	This form can be interfaced with the Omega [tridiagonal solver](TridiagonalSolver.md) routine.

[cbegeman]

Suggested change

-\frac{\rho_0}{\left[\tilde{h}_k\right]_e} \left\{ \left[\left<u^\prime \tilde{w}_{tr}^\prime\right> \right]_{e,k} - \left[\left<u^\prime \tilde{w}_{tr}^\prime\right> \right]_{e,k+1} \right\} = -\frac{\rho_0}{\left[\tilde{h}_k\right]_e} \left\{ \frac{\left(u_{e,k-1} - u_{e,k}\right)}{0.5 \left(\tilde{h}_k-1 + \tilde{h}_{k}\right)} - \frac{\left(u_{e,k} - u_{e,k+1}\right)}{0.5 \left(\tilde{h}_k + \tilde{h}_{k+1}\right)} \right\}.

-\frac{\rho_0}{\left[\tilde{h}_k\right]_e} \left\{ \left[\left<u^\prime \tilde{w}_{tr}^\prime\right> \right]_{e,k} - \left[\left<u^\prime \tilde{w}_{tr}^\prime\right> \right]_{e,k+1} \right\} = -\frac{\rho_0}{\left[\tilde{h}_k\right]_e} \left\{ \frac{\left(u_{e,k-1} - u_{e,k}\right)}{0.5 \left(\tilde{h}_{k-1} + \tilde{h}_{k}\right)} - \frac{\left(u_{e,k} - u_{e,k+1}\right)}{0.5 \left(\tilde{h}_k + \tilde{h}_{k+1}\right)} \right\}.

[cbegeman]

Suggested change

	Since, Omega will predict layer average quantities (e.g., $u_k$), it's not immediately clear that a traditional discretization is appropriate as there is a discontinuity in the predicted variables at the layer edge. To circumvent this problem, we turn to weak derivatives instead of traditional pointwise forms. A weak derivative of $\varphi$ is defined as
	Since, Omega will predict layer average quantities (e.g., $u_k$), it's not immediately clear that a traditional discretization is appropriate as there is a discontinuity in the predicted variables at layer interfaces. To circumvent this problem, we turn to weak derivatives instead of traditional pointwise forms. A weak derivative of $\varphi$ is defined as

[cbegeman]

Suggested change

	If [](#weak-derivative) is averaged between $\tilde{z}_{k-1/2}^{\text{top}}$ and $\tilde{z}_{k+1/2}^{\text{top}}$, we arrive at the expected form of the gradient
	If [](#weak-derivative) is averaged between $\tilde{z}_{k-1/2}^{\text{top}}$ and $\tilde{z}_{k+1/2}^{\text{top}}$, we arrive at the expected form of the gradient centered on layer interfaces

[cbegeman]

I don't know if this level of detail is important here, but in MPAS-O, the $\left|u_{e,k}\right|$ term is replaced with the kinetic energy, computed first on cells and then averaged back to edges. I imagine this was motivated by a desire to have a large stencil (velocity on all edges spanning two cells). Perhaps coincidentally, this is also the stencil that is used for sea ice drag.

[hyungyukang]

I’m approving this document based on my review, edits, and the feedback from others.
It’s great to have such a detailed document for the governing equations. This is an excellent work! Many thanks to everyone who contributed to this document.