Re-implementing the HydrostaticFreeSurfaceModel with Generalized Vertical Coordinates (GVC)

[glwagner]

The equations of motion for most or all major ocean modeling software are implemented using a "Generalized Vertical Coordinate" (GVC). Generalized vertical coordinates contain a "fixed z" coordinate as a limiting case, but generalize to vertical coordinates that

• Are fixed in time but conform to topographic variations (σ coordinates);
• Vary in time in a diagnostic manner (for example, z-star coordinates, which vary with the sea surface displacement);
• Vary in time according to some prognostic equation (fully Lagrangian isopycnal coordinates, semi-Lagrangian hybrid coordinates / z-tilde coordinates that vary only with high frequency motion, and adaptive methods that prescribe GVC dynamics to obtain favorable properties like concentration in regions of strong stratification).

A fully general GVC typically also requires a "Lagrangian remapping" step to avoid extreme grid distortions in regions of persistent vertical velocities.

The implementation of GVC is likely a major refactoring of HydrostaticFreeSurfaceModel because it will change the equations of motion and could even potentially change the nature of its prognostic variables. For time-dependent GVC, tracer equations must be reformulated in terms of the "thickness-weighted" tracer concentration, which in our case means the tracer concentration normalized by the local grid spacing. This could mean either using the thickness-weighted tracer conservation as a state variable, or it could mean rewriting the time-stepping algorithm so that unweighted tracer conservation can be updated according to the conservation of thickness-weighted tracer.

A preliminary roadmap towards GVC in HydrostaticFreeSurfaceModel is

1. Introduce AbstractVerticalCoordinate and refactor the HydrostaticFreeSurfaceModel to integrate thickness-weighted equations. When using ZCoordinate, the resulting model produces identical results to the current implementation of HydrostaticFreeSurfaceModel. I think we will also want a rectilinear grid to use for testing to avoid modifying RegularRectilinearGrid, such as a refactored or reimplemented three-dimensionally StretchedRectilinearGrid from #1532 which accepts the use of a generalized vertical coordinates.

2. Introduce ZStarCoordinate (or maybe FreeSurfaceWeightedVerticalCoordinate). "z-star" vertical coordinate is a relatively simple GVC, with a small, diagnostic time-dependence. A successful implementation will require integrating thickness-weighted equations, correct vertical derivatives, and ensuring correct horizontal pressure gradients.

3. Implement remapping as an alternative to vertical advection and test using ZCoordinate and ZStarCoordinate.

Of course, steps 3 might evolve depending on our experience during step 2.

An important choice that's part of step 1. is whether to use thickness-weighted variables as state variables rather than unweighted variables (eg hu rather than u, and hc rather than c. I suspect it makes sense to use thickness-weighted state variables because this simplifies the underlying implementation (for example, we only need hⁿ rather than both hⁿ⁺¹ and hⁿ to evolve a tracer field, and we can evaluate the prognostic equations in the same order that we do now), and because it is easy to provide abstractions for users to compute, inspect, and output unweighted variables u = hu / h. The pros and cons of this important choice should be discussed.

The discussion on #1549 is related (I decided to start a new issue because GVC are more general than terrain-following coordinates).

Some references

• MOM6 description (Adcroft et al 2019)
• Primer on “remap” methods (Griffies et al, 2019)
• “z-tilde” GVC used by the NEMO ocean model (Leclair and Madec 2011)
• Evaluation of various GVC in MPAS-Ocean (Petersen et al 2015)
• GVC in a nonhydrostatic compressible atmospheric model ENDgame (Kavcic and Thuburn, 2018)
• Correct calculation of pressure gradients in a finite volume model (Adcroft et al 2008) An “adaptive”, semi-Lagrangian GVC (Burchard and Beckers 2004) See also a 3D description (Hofmeister et al 2010) and realistic test case (Hofmeister et al 2011) “z star” GVC in MITgcm (Adcroft and Campin 2004)

[francispoulin]

Sounds very exciting! As I mentioned, I have a student who would be very interested in working on the sigma-coordinate model, which would hopefully help with the end goal.

[glwagner]

A thorny detail is accurate calculate of pressure gradients for strongly warped coordinates.

Looking into the literature a little bit, it appears that the current pressure gradient force formulation may be sufficient for "weak" GVC, which include the "z-star" coordinate introduced by Adcroft and Campin (2004) (rescaling of the vertical coordinate by sea surface displacement) and the "z-tilde" coordinate introduced by Leclair and Madec (2011) (semi-Lagrangian vertical coordinate that follows high-pass filtered fast motions only and is restored to a target grid on a relatively short time scale).

Other formulations, including classical sigma coordinates and quasi-Lagrangian methods that involve grid warping severe enough to require remapping, may fail unless we improve our method for calculating the horizontal pressure gradient force. In particular, the method we use now is essentially finite difference and requires a vertical coordinate that exactly or "almost" coincides with a geopotential surface.

Finite volume treatment of the pressure gradient force is discussed by

• Lin (1997) in the context of an atmospheric model with an effectively linear equation of state
• Shchepetkin and McWilliams (2003) which implements a high-order method for evaluating a contour integral around momentum elements for computing the pressure gradient force
• Adcroft et al (2008) that utilizes a crucial "analytical integration" step valid for for a nonlinear equation of state that can be written in a "simplified" form, which avoids the relatively more expensive numerical integration of part of the contour (and is more accurate)
• Engwirda et al (2017) using high-order numerical integration techniques for nonlinear equations of state

[jm-c]

Just a detail:sigma coordinate does also moves in time with SSH (when the bottom is flat, sigma and z* are the same).

[glwagner]

Ah, that makes sense since the rescaling sets the surface (z=eta) as sigma=0, rather than z=0. Thanks for that clarification @jm-c !

[glwagner]

@simone-silvestri here's some description of the GVC problem. There's additional discussion on #1549.

Today, we discussed whether the prognostic variable in a model formulated with a generalized vertical coordinate is thickness weighted velocity hu, or velocity u (which has also come up on #2522).

Based on some discussion and consultation with papers it seems that we want both, because we evolve hu when we use "flux form" or "conservative form" momentum advection, but we evolve u when we use "vector invariant" advection. This choice has consequences beyond momentum advection though (I think).

Also, whichever momentum variable we prognose, we always prognose the thickness weighted tracer concentration hc.

Thoughts on this are very welcome @jm-c @jmbeckers !

In the "z-tilde" paper by Leclair and Madec 2011, we find

I think e3 above is our h.