_ALE.dox

/*! \page ALE Vertical Lagrangian method: conceptual

\section section_ALE Lagrangian and ALE

As discussed by Adcroft and Hallberg (2008) \cite adcroft2006 and

Griffies, Adcroft and Hallberg (2020) \cite Griffies_Adcroft_Hallberg2020,

we can conceive of two general classes

of algorithms that frame how hydrostatic ocean models are

formulated. The two classes differ in how they treat the vertical

direction. Quasi-Eulerian methods follow the approach traditionally

used in geopotential coordinate models, whereby vertical motion is

diagnosed via the continuity equation. Quasi-Lagrangian methods are

traditionally used by layered isopycnal models, with the vertical

Lagrangian approach specifying motion that crosses coordinate

surfaces. Indeed, such dia-surface flow can be set to zero using

Lagrangian methods for studies of adiabatic dynamics. MOM6 makes use

of the vertical Lagrangian remap method, as pioneered for ocean

modeling by Bleck (2002) \cite bleck2002 and further documented by

\cite Griffies_Adcroft_Hallberg2020, with this method a limit case of

the Arbitrary-Lagrangian-Eulerian method (\cite hirt1997). Dia-surface

transport is implemented via a remapping so that the method can be

summarized as the Lagrangian plus remap approach and so it is a

one-dimensional version of the incremental remapping of Dukowicz (2000)

\cite dukowicz2000.

\image html  ALE_general_schematic.png "Schematic of the 3d Lagrangian regrid/remap method"  width=70%

\image latex ALE_general_schematic.png "Schematic of the 3d Lagrangian regrid/remap method"  width=0.7\textwidth

Refer to the above figure taken from Griffies, Adcroft, and Hallberg

(2020) \cite Griffies_Adcroft_Hallberg2020.  It shows a schematic of

the Lagrangian-remap method as well as the Arbitrary

Lagrangian-Eulerian (ALE) method. The first panel shows a square fluid

region and square grid used to represent the fluid, along with

rectangular subregions partitioned by grid lines. The second panel

shows the result of evolving the fluid region and evolving the

grid. The grid can evolve according to the fluid flow, as per a

Lagrangian method, or it can evolve according to some specified grid

evolution, as per an ALE method. The right panel depicts the grid

reinitialization onto a target grid (the regrid step). A regrid step

necessitates a corresponding remap step to estimate the ocean state on

the target grid, with conservative remapping required to preserve

integrated scalar contents (e.g., potential enthalpy, salt mass, and

seawater mass). The regrid/remap steps are needed for Lagrangian

methods in order for the grid to retain an accurate representation of

the ocean state.  Ideally, the remap step does not affect any changes

to the fluid state; rather, it only modifies where in space the fluid

state is represented. However, any numerical realization incurs

interpolation inaccuracies that lead to unphysical (spurious) state

changes.

\section section_ALE_MOM Vertical Lagrangian regrid/remap method

We now get a bit more specific to the vertical Lagrangian method.

For this purpose, recall recall the basic dynamical equations (those

equations with a time derivative) of MOM6 discussed in

\ref General_Coordinate

\f{align}

\rho_0

\left[ \frac{\partial \mathbf{u}}{\partial t} + \frac{( f + \zeta )}{h} \,

\hat{\mathbf{z}} \times h \, \mathbf{u} + \underbrace{ \dot{r} \,

\frac{\partial \mathbf{u}}{\partial r} }

\right]

&= -\nabla_r \, (p + \rho_{0} \, K) -

\rho \nabla_r \, \Phi + \mathbf{\mathcal{F}}

&\mbox{horizontal momentum}

\label{eq:h-horz-momentum-vlm}

\\

\frac{\partial h}{\partial t} + \nabla_r \cdot \left( h \, \mathbf{u} \right) +

\underbrace{ \delta_r ( z_r \dot{r} ) }

 &= 0

&\mbox{thickness}

\label{eq:h-thickness-equation-vlm}

\\

\frac{\partial ( \theta \, h )}{\partial t} + \nabla_r \cdot \left( \theta h \,

\mathbf{u} \right) + \underbrace{ \delta_r ( \theta \, z_r \dot{r} ) }

&=

h \mathbf{\mathcal{N}}_\theta^\gamma - \delta_r J_\theta^{(z)}

&\mbox{potential/Conservative temp}

\label{eq:h-temperature-equation-vlm}

 \\

\frac{\partial ( S \, h )}{\partial t} + \nabla_r \cdot \left( S \, h \,

\mathbf{u} \right) + \underbrace{ \delta_r ( S \, z_r \dot{r} ) }

 &=

h \mathbf{\mathcal{N}}_S^\gamma - \delta_r J_S^{(z)}

&\mbox{salinity}

\label{eq:h-salinity-equation-vlm}

\f}

The MOM6 implementation of the vertical Lagrangian method makes

use of two general steps. The first evolves the ocean state forward in

time according to a vertical Lagrangian approach with with

\f$\dot{r}=0\f$. Hence, the horizontal momentum, thickness, and

tracers are time stepped with the underbraced terms removed in the

above equations. All advective transport occurs within a layer as

defined by constant \f$r\f$-surfaces so that the volume within each

layer is fixed. All other terms are retained in their full form,

including subgrid scale terms that contribute to the transfer of

tracer and momentum into distinct \f$r\f$ layers (e.g., dia-surface

diffusion of tracer and velocity). Maintaining constant volume within

a layer yet allowing for tracers to move between layers engenders no

inconsistency between tracer and thickness evolution. The reason is

that tracer diffusion, even dia-surface diffusion, does not transfer

volume.

The second step in the method comprises the generation of a new

vertical grid following a prescription, such as whether the grid

should align with isopcynals or constant \f$z^{*}\f$ or a combination.

This second step is known as the regrid step.  The ocean state is then

vertically remapped to the newly generated vertical grid. This

remapping step incorporates dia-surface transfer of properties, with

such transfer depending on the prescription given for the vertical

grid generation. To minimize discretization errors and the associated

spurious mixing, the remapping step makes use of the high order

accurate methods developed by \cite white2008 and \cite white2009.

\section section_ALE_MOM_numerics Outlining the numerical algorithm

The underlying algorithm for treatment of the vertical can be related

to operator-splitting of the underbraced terms in the above equations.

If we consider, for simplicity, an Euler-forward update for a

time-step \f$\Delta t\f$, the time-stepping for the thickness and

tracer equation (\f$C\f$ is an arbitrary tracer) can be summarized as

(from Table 1 in Griffies, Adcroft and Hallberg (2020)

\cite Griffies_Adcroft_Hallberg2020)

\f{align}

\label{html:ale-equations}\notag

\\

 \delta_{r} w^{\scriptstyle{\mathrm{grid}}}

 &= -\nabla_{r} \cdot [h \, \mathbf{u}]^{(n)}

 &\mbox{layer motion via horz conv}

\\

 h^{\dagger} &= h^{(n)} + \Delta t \, \delta_{r} w^{\scriptstyle{\mathrm{grid}}}

= h^{(n)} - \Delta t \, \nabla_{r} \cdot [h \, \mathbf{u}]^{(n)}

 &\mbox{update thickness via horz advect}

\\

 [h \, C]^{\dagger} &= [h \, C]^{(n)} -\Delta t \, \nabla_{r} \cdot [ h \, C \, \mathbf{u} ]^{(n)}

 &\mbox{update tracer via horz advect}

\\

 h^{(n+1)} &=  h^{\scriptstyle{\mathrm{target}}}

 &\mbox{regrid to the target grid}

\\

 \delta_{r} w^{(\dot{r})} &= -(h^{\scriptstyle{\mathrm{target}}}  -  h^{\dagger})/\Delta t

 &\mbox{diagnose dia-surface transport}

\\

 [h \, C]^{(n+1)} &= [h \, C]^{\dagger}  - \Delta t  \, \delta_{r} ( w^{(\dot{r})} \, C^{\dagger})

 &\mbox{remap tracer via dia-surface transport}

\f}

The first three equations constitute the Lagrangian portion of the

algorithm.  In particular, the second equation provides an

intermediate or predictor value for the updated thickness,

\f$h^{\dagger}\f$, resulting from the vertical Lagrangian update.

Similarly, the third equation performs a Lagrangian update of the

thickness-weighted tracer to intermediate values, again operationally

realized by dropping the \f$w^{(\dot{r})}\f$ contribution.

The fourth equation is the regrid step, which is the key step in the

algorithm with the new grid defined by the new thickness

\f$h^{(n+1)}\f$.  The new thickness is prescribed by the target values

for the vertical grid,

\f{align}

 h^{(n+1)} = h^{\scriptstyle{\mathrm{target}}}.

\f}

The prescribed target grid thicknesses are then used to diagnose the

dia-surface velocity according to

\f{align}

 \delta_{r} w^{(\dot{r})} = -(h^{\scriptstyle{\mathrm{target}}} - h^{\dagger})/\Delta t.

\f}

This step, and the remaining step for tracers, constitute the

remapping portion of the algorithm.  For example, if the prescribed

coordinate surfaces are geopotentials, then \f$w^{(\dot{r})}\f$ and

\f$h^{\scriptstyle{\mathrm{target}}} = h^{(n)}\f$, in which case the

remap step reduces to Cartesian vertical advection.

Within the above framework for evolving the ocean state, we make use

of a standard split-explicit time stepping method by decomposing the

horizontal momentum equation into its fast (depth integrated) and slow

(deviation from depth integrated) components. Furthermore, we follow

the methods of Hallberg and Adcroft (2009) \cite hallberg2009 to

ensure that the free surface resulting from time stepping the depth

integrated thickness equation (i.e., the free surface equation) is

consistent with the sum of the thicknesses that result from time

stepping the layer thickness equations for each of the discretized

layers; i.e., \f$\sum_{k} h = H + \eta\f$.

*/