Multilayer VAM

MLVAM is a non-hydrostatic model for dispersive free-surface flow. The depth is split into n_layers layers that move with the free surface, and each layer is a vertically averaged and moment (VAM) column: horizontal velocity, vertical velocity and non-hydrostatic pressure as polynomials in the layer. Papers: the VAM column of Escalante, Morales de Luna, Cantero-Chinchilla and Castro-Orgaz 2024, the layer decomposition with interface mass transfer of Escalante, Fernández-Nieto, Morales de Luna and Castro 2019.

Starting point

The incompressible Navier-Stokes system in \((x, z)\) for velocity \((u, w)\), total pressure \(p_T\) and deviatoric stress \(\sigma\):

\[\begin{split} \begin{cases} \partial_x u + \partial_z w = 0, \\ \partial_t u + u \partial_x u + w \partial_z u + \partial_x p_T = \partial_x \sigma_{xx} + \partial_z \sigma_{xz}, \\ \partial_t w + u \partial_x w + w \partial_z w + \partial_z p_T = -g + \partial_x \sigma_{zx} + \partial_z \sigma_{zz} \end{cases} \end{split}\]

with the stress at the free surface \(\eta = h + b\) and at the bed \(b\), the hydrostatic and non-hydrostatic pressure, and the kinematic conditions

\[ \sigma_{xz}|_{z=\eta} = 0, \quad \sigma_{xz}|_{z=b} = \tau, \qquad p_T = p_H + p, \quad \partial_z p_H = -g, \qquad u|_{z=b} \partial_x b = w|_{z=b}, \quad \partial_t \eta + u|_{z=\eta} \partial_x \eta = w|_{z=\eta} \]

The depth is split into layers of fixed fractions, which exchange mass across the interface (two layers in the paper):

\[\begin{split} h_\alpha = l_\alpha h, \quad l_1 + l_2 = 1, \qquad \begin{cases} \partial_t h_1 + \partial_x (h_1 u_1) = -\Gamma_I, \\ \partial_t h_2 + \partial_x (h_2 u_2) = \Gamma_I, \end{cases} \qquad \Gamma_I = l_1 \partial_x (h_2 u_2) - l_2 \partial_x (h_1 u_1) \end{split}\]

In each layer, on the Legendre basis \(\phi_j\) of the layer coordinate \(\xi \in [0, 1]\) with \(\phi_j(0) = 1\):

\[ u = u_0\phi_0 + u_1\phi_1, \quad w = w_0\phi_0 + w_1\phi_1 + w_2\phi_2, \quad p = p_0\phi_0 + p_1\phi_1 + p_2\phi_2 \]

The equations are the mass and momentum balances tested with \(\phi_j\) and integrated over the layer.

Assumptions

• The in-plane viscous stresses are dropped, only \(\sigma_{xz}\) is kept.

• Per layer \(u \in \mathbb{P}_1\) and \(w, p \in \mathbb{P}_2\). Zoomy generalizes this to \(u \in \mathbb{P}_N\), \(w, p \in \mathbb{P}_{N+1}\) with \(N\) = level.

• The non-hydrostatic pressure vanishes at the free surface.

• The layers are fixed fractions of the depth. Zoomy takes any number of layers and keeps the pressure continuous across an interface, the paper allows a jump (its (16)).

• The interface velocity of the exchange term is the mean of the two layers (the paper’s choice, Zoomy’s default MeanInterface).

Generated equations

MLVAM(n_layers=2, level=1, dimension=2, closures=[Newtonian(), NavierSlip(), StressFree()])

State \(Q\): \(b\), \(h\), \(q_{1,0}\), \(q_{1,1}\), \(q_{2,0}\), \(q_{2,1}\), \(r_{1,0}\), \(r_{1,1}\), \(r_{2,0}\), \(r_{2,1}\), \(P_{1,0}\), \(P_{1,1}\), \(P_{2,0}\), \(P_{2,1}\)

Parameters: g = 9.81, rho = 1.0, nu = 0.0, lambda_s = 0.0, e_x = 0.0, l_1 = 0.5

Bed and total depth

bottom

\[ \partial_{t} b = 0 \]

continuity

\[ \partial_{t} h + \partial_{x} q_{1,0} + \partial_{x} q_{2,0} = 0 \]

Layer 1

momentum_x_1_0

\[\begin{split} \begin{aligned} &\partial_{t} q_{1,0} + \partial_{x} \left(\frac{g l_{1}^{2} h^{2}}{2}\right) + \partial_{x} \left(\frac{q_{1,0}^{2}}{l_{1} h}\right) + \partial_{x} \left(\frac{q_{1,1}^{2}}{3 l_{1} h}\right) \\ &\quad + \partial_{x} \left(\frac{l_{1} P_{1,0} h}{\rho}\right) + \frac{g l_{1} \rho h - 2 P_{1,1}}{\rho} \, \partial_{x} b \\ &\quad - \frac{l_{1} \left(g l_{1} \rho h - g \rho h - 2 P_{2,1}\right)}{\rho} \, \partial_{x} h \\ &\quad + \frac{l_{1} q_{1,0} + l_{1} q_{1,1} - l_{1} q_{2,0} + l_{1} q_{2,1} - q_{1,0} - q_{1,1}}{2 l_{1} h} \, \partial_{x} q_{1,0} \\ &\quad + \frac{l_{1} q_{1,0} + l_{1} q_{1,1} - l_{1} q_{2,0} + l_{1} q_{2,1} - q_{1,0} - q_{1,1}}{2 \left(l_{1} - 1\right) h} \, \partial_{x} q_{2,0} - \frac{2 \nu q_{1,1}}{l_{1}^{2} h^{2}} \\ &\quad + \frac{\lambda_{s} q_{1,0}}{l_{1} \rho h} - \frac{\lambda_{s} q_{1,1}}{l_{1} \rho h} - e_{x} g l_{1} h = 0 \end{aligned} \end{split}\]

momentum_x_1_1

\[\begin{split} \begin{aligned} &\partial_{t} q_{1,1} + \partial_{x} \left(\frac{l_{1} P_{1,1} h}{\rho} + \frac{2 q_{1,0} q_{1,1}}{l_{1} h}\right) + \frac{6 \left(P_{1,0} + P_{1,1} + 2 P_{2,1}\right)}{\rho} \, \partial_{x} b \\ &\quad + \frac{l_{1} \left(3 P_{1,0} + P_{1,1} + 6 P_{2,1}\right)}{\rho} \, \partial_{x} h \\ &\quad - \frac{3 l_{1} q_{1,0} - l_{1} q_{1,1} + 3 l_{1} q_{2,0} - 3 l_{1} q_{2,1} - 3 q_{1,0} + q_{1,1}}{2 l_{1} h} \, \partial_{x} q_{1,0} - \frac{q_{1,0}}{l_{1} h} \, \partial_{x} q_{1,1} \\ &\quad - \frac{3 l_{1} q_{1,0} - l_{1} q_{1,1} + 3 l_{1} q_{2,0} - 3 l_{1} q_{2,1} - 3 q_{1,0} + q_{1,1}}{2 \left(l_{1} - 1\right) h} \, \partial_{x} q_{2,0} + \frac{6 \nu q_{1,1}}{l_{1}^{2} h^{2}} \\ &\quad - \frac{3 \lambda_{s} q_{1,0}}{l_{1} \rho h} + \frac{3 \lambda_{s} q_{1,1}}{l_{1} \rho h} = 0 \end{aligned} \end{split}\]

momentum_z_1_0

\[\begin{split} \begin{aligned} &\partial_{t} r_{1,0} + \partial_{x} \left(\frac{q_{1,0} r_{1,0}}{l_{1} h}\right) + \partial_{x} \left(\frac{q_{1,1} r_{1,1}}{3 l_{1} h}\right) \\ &\quad + \frac{\left(l_{1} - 1\right) \left(q_{1,0} + q_{1,1}\right)}{l_{1} h} \, \partial_{x} b \partial_{x} q_{1,0} + \frac{q_{1,0} + q_{1,1}}{h} \, \partial_{x} b \partial_{x} q_{2,0} \\ &\quad + \frac{\left(l_{1} - 1\right) \left(q_{1,0} + q_{1,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{1,0} + \frac{l_{1} \left(q_{1,0} + q_{1,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{2,0} \\ &\quad - l_{1} \, \partial_{x} q_{1,0} \partial_{x} q_{2,0} - \left(l_{1} - 1\right) \left(\partial_{x} q_{1,0}\right)^{2} + \frac{2 P_{1,1}}{\rho} = 0 \end{aligned} \end{split}\]

momentum_z_1_1

\[\begin{split} \begin{aligned} &\partial_{t} r_{1,1} + \partial_{x} \left(\frac{5 q_{1,0} r_{1,1} + 3 q_{1,1} r_{1,0} + 2 q_{1,1} r_{1,1}}{5 l_{1} h}\right) \\ &\quad + \partial_{x} \left(\frac{2 \left(q_{1,0} - q_{1,1}\right) q_{1,1} \partial_{x} b}{5 l_{1} h}\right) - \frac{\left(l_{1} - 1\right) \left(3 r_{1,0} + r_{1,1}\right)}{l_{1} h} \, \partial_{x} q_{1,0} \\ &\quad - \frac{6 r_{1,0} - r_{1,1}}{5 l_{1} h} \, \partial_{x} q_{1,1} - \frac{3 r_{1,0} + r_{1,1}}{h} \, \partial_{x} q_{2,0} \\ &\quad + \frac{3 \left(l_{1} - 1\right) \left(q_{1,0} + q_{1,1}\right)}{l_{1} h} \, \partial_{x} b \partial_{x} q_{1,0} + \frac{q_{1,0} - q_{1,1}}{5 l_{1} h} \, \partial_{x} b \partial_{x} q_{1,1} \\ &\quad + \frac{3 \left(q_{1,0} + q_{1,1}\right)}{h} \, \partial_{x} b \partial_{x} q_{2,0} + \frac{3 \left(l_{1} - 1\right) \left(q_{1,0} + q_{1,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{1,0} \\ &\quad + \frac{3 l_{1} \left(q_{1,0} + q_{1,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{2,0} - 3 l_{1} \, \partial_{x} q_{1,0} \partial_{x} q_{2,0} \\ &\quad - 3 \left(l_{1} - 1\right) \, \left(\partial_{x} q_{1,0}\right)^{2} - \frac{12 P_{2,1}}{\rho} - \frac{6 P_{1,0}}{\rho} - \frac{6 P_{1,1}}{\rho} = 0 \end{aligned} \end{split}\]

constraint_1_0

\[\begin{split} \begin{aligned} &\partial_{x} \left(- q_{1,0}\right) + \frac{l_{1} h}{3} \, \partial_{x} \left(\frac{q_{1,1}}{l_{1} h}\right) + \frac{2 q_{1,0}}{l_{1} h} \, \partial_{x} b \\ &\quad + \frac{3 q_{1,0} + 2 q_{1,1}}{3 h} \, \partial_{x} h - \frac{2 r_{1,0}}{l_{1} h} = 0 \end{aligned} \end{split}\]

constraint_1_1

\[ \partial_{x} \left(- q_{1,0}\right) + \frac{2 q_{1,1}}{l_{1} h} \, \partial_{x} b + \frac{q_{1,0} + q_{1,1}}{h} \, \partial_{x} h - \frac{2 r_{1,1}}{l_{1} h} = 0 \]

Layer 2

momentum_x_2_0

\[\begin{split} \begin{aligned} &\partial_{t} q_{2,0} + \partial_{x} \left(- \frac{\left(l_{1} - 1\right) P_{2,0} h}{\rho}\right) + \partial_{x} \left(\frac{g l_{1}^{2} h^{2}}{2} - g l_{1} h^{2} + \frac{g h^{2}}{2}\right) \\ &\quad + \partial_{x} \left(\frac{q_{2,1}^{2}}{- 3 l_{1} h + 3 h}\right) + \partial_{x} \left(\frac{q_{2,0}^{2}}{- l_{1} h + h}\right) - \frac{g l_{1} \rho h - g \rho h + 2 P_{2,1}}{\rho} \, \partial_{x} b \\ &\quad - \frac{l_{1} \left(g l_{1} \rho h - g \rho h + 2 P_{2,1}\right)}{\rho} \, \partial_{x} h \\ &\quad - \frac{l_{1} q_{1,0} + l_{1} q_{1,1} - l_{1} q_{2,0} + l_{1} q_{2,1} - q_{1,0} - q_{1,1}}{2 l_{1} h} \, \partial_{x} q_{1,0} \\ &\quad - \frac{l_{1} q_{1,0} + l_{1} q_{1,1} - l_{1} q_{2,0} + l_{1} q_{2,1} - q_{1,0} - q_{1,1}}{2 \left(l_{1} - 1\right) h} \, \partial_{x} q_{2,0} \\ &\quad + \frac{2 \nu q_{2,1}}{l_{1}^{2} h^{2} - 2 l_{1} h^{2} + h^{2}} - \frac{e_{x} g h}{l_{1}^{2} - 2 l_{1} + 1} + \frac{e_{x} g l_{1}^{3} h}{l_{1}^{2} - 2 l_{1} + 1} \\ &\quad - \frac{3 e_{x} g l_{1}^{2} h}{l_{1}^{2} - 2 l_{1} + 1} + \frac{3 e_{x} g l_{1} h}{l_{1}^{2} - 2 l_{1} + 1} = 0 \end{aligned} \end{split}\]

momentum_x_2_1

\[\begin{split} \begin{aligned} &\partial_{t} q_{2,1} + \partial_{x} \left(- \frac{l_{1}^{2} P_{2,1} h^{2} - 2 l_{1} P_{2,1} h^{2} + 2 \rho q_{2,0} q_{2,1} + P_{2,1} h^{2}}{\rho \left(l_{1} - 1\right) h}\right) \\ &\quad + \frac{6 \left(P_{2,0} + P_{2,1}\right)}{\rho} \, \partial_{x} b + \frac{3 l_{1} P_{2,0} + 5 l_{1} P_{2,1} + 3 P_{2,0} + P_{2,1}}{\rho} \, \partial_{x} h \\ &\quad + \frac{3 l_{1} q_{1,0} + 3 l_{1} q_{1,1} + 3 l_{1} q_{2,0} + l_{1} q_{2,1} - 3 q_{1,0} - 3 q_{1,1}}{2 l_{1} h} \, \partial_{x} q_{1,0} \\ &\quad + \frac{3 l_{1} q_{1,0} + 3 l_{1} q_{1,1} + 3 l_{1} q_{2,0} + l_{1} q_{2,1} - 3 q_{1,0} - 3 q_{1,1}}{2 \left(l_{1} - 1\right) h} \, \partial_{x} q_{2,0} \\ &\quad + \frac{q_{2,0}}{\left(l_{1} - 1\right) h} \, \partial_{x} q_{2,1} + \frac{6 \nu q_{2,1}}{l_{1}^{2} h^{2} - 2 l_{1} h^{2} + h^{2}} = 0 \end{aligned} \end{split}\]

momentum_z_2_0

\[\begin{split} \begin{aligned} &\partial_{t} r_{2,0} + \partial_{x} \left(\frac{q_{2,1} r_{2,1}}{- 3 l_{1} h + 3 h}\right) + \partial_{x} \left(\frac{q_{2,0} r_{2,0}}{- l_{1} h + h}\right) \\ &\quad + \frac{q_{2,0} - q_{2,1}}{h} \, \partial_{x} b \partial_{x} q_{1,0} + \frac{l_{1} \left(q_{2,0} - q_{2,1}\right)}{\left(l_{1} - 1\right) h} \, \partial_{x} b \partial_{x} q_{2,0} \\ &\quad + \frac{l_{1} \left(q_{2,0} - q_{2,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{1,0} \\ &\quad + \frac{l_{1}^{2} \left(q_{2,0} - q_{2,1}\right)}{\left(l_{1} - 1\right) h} \, \partial_{x} h \partial_{x} q_{2,0} + l_{1} \, \partial_{x} q_{1,0} \partial_{x} q_{2,0} \\ &\quad + \left(l_{1} - 1\right) \left(\partial_{x} q_{1,0}\right)^{2} + \frac{2 P_{2,1}}{\rho} = 0 \end{aligned} \end{split}\]

momentum_z_2_1

\[\begin{split} \begin{aligned} &\partial_{t} r_{2,1} \\ &\quad + \partial_{x} \left(- \frac{2 l_{1} h q_{2,1} \partial_{x} q_{1,0} + 2 l_{1} q_{2,0} q_{2,1} \partial_{x} h - 2 l_{1} q_{2,1}^{2} \partial_{x} h - 2 h q_{2,1} \partial_{x} q_{1,0} + 5 q_{2,0} r_{2,1} + 3 q_{2,1} r_{2,0} + 2 q_{2,1} r_{2,1}}{5 \left(l_{1} - 1\right) h}\right) \\ &\quad + \partial_{x} \left(- \frac{2 \left(q_{2,0} - q_{2,1}\right) q_{2,1} \partial_{x} b}{5 \left(l_{1} - 1\right) h}\right) + \frac{3 r_{2,0} - r_{2,1}}{h} \, \partial_{x} q_{1,0} \\ &\quad + \frac{l_{1} \left(3 r_{2,0} - r_{2,1}\right)}{\left(l_{1} - 1\right) h} \, \partial_{x} q_{2,0} + \frac{6 r_{2,0} - r_{2,1}}{5 \left(l_{1} - 1\right) h} \, \partial_{x} q_{2,1} \\ &\quad - \frac{3 \left(q_{2,0} - q_{2,1}\right)}{h} \, \partial_{x} b \partial_{x} q_{1,0} - \frac{3 l_{1} \left(q_{2,0} - q_{2,1}\right)}{\left(l_{1} - 1\right) h} \, \partial_{x} b \partial_{x} q_{2,0} \\ &\quad - \frac{q_{2,0} - q_{2,1}}{5 \left(l_{1} - 1\right) h} \, \partial_{x} b \partial_{x} q_{2,1} - \frac{3 l_{1} \left(q_{2,0} - q_{2,1}\right)}{h} \, \partial_{x} h \partial_{x} q_{1,0} \\ &\quad - \frac{3 l_{1}^{2} \left(q_{2,0} - q_{2,1}\right)}{\left(l_{1} - 1\right) h} \, \partial_{x} h \partial_{x} q_{2,0} \\ &\quad - \frac{l_{1} \left(q_{2,0} - q_{2,1}\right)}{5 \left(l_{1} - 1\right) h} \, \partial_{x} h \partial_{x} q_{2,1} - 3 l_{1} \, \partial_{x} q_{1,0} \partial_{x} q_{2,0} \\ &\quad - \frac{1}{5} \, \partial_{x} q_{1,0} \partial_{x} q_{2,1} - 3 \left(l_{1} - 1\right) \, \left(\partial_{x} q_{1,0}\right)^{2} - \frac{6 P_{2,0}}{\rho} - \frac{6 P_{2,1}}{\rho} = 0 \end{aligned} \end{split}\]

constraint_2_0

\[\begin{split} \begin{aligned} &\partial_{x} \left(- q_{2,0}\right) - \frac{\left(l_{1} - 1\right) h}{3} \, \partial_{x} \left(\frac{q_{2,1}}{- l_{1} h + h}\right) - \frac{2 q_{2,0}}{\left(l_{1} - 1\right) h} \, \partial_{x} b \\ &\quad - \frac{3 l_{1} q_{2,0} - 2 l_{1} q_{2,1} + 3 q_{2,0} + 2 q_{2,1}}{3 \left(l_{1} - 1\right) h} \, \partial_{x} h - 2 \, \partial_{x} q_{1,0} + \frac{2 r_{2,0}}{l_{1} h - h} = 0 \end{aligned} \end{split}\]

constraint_2_1

\[\begin{split} \begin{aligned} &\partial_{x} \left(- q_{2,0}\right) - \frac{2 q_{2,1}}{\left(l_{1} - 1\right) h} \, \partial_{x} b + \frac{l_{1} q_{2,0} - l_{1} q_{2,1} - q_{2,0} - q_{2,1}}{\left(l_{1} - 1\right) h} \, \partial_{x} h \\ &\quad + \frac{2 r_{2,1}}{l_{1} h - h} = 0 \end{aligned} \end{split}\]

Notation: layer \(\ell\) has depth \(h_\ell = l_\ell h\) (parameter l_1, \(l_2 = 1 - l_1\)), \(q_{\ell,k} = (-1)^k h_\ell u_{\ell,k}\), \(r_{\ell,k} = (-1)^k h_\ell w_{\ell,k}\), \(P_{\ell,k} = (-1)^k \rho p_{\ell,k}\) (Zoomy’s basis has \(\phi_k(1) = 1\)). The top modes \(w_{\ell,N+1}\) are closed upward from the bed kinematic condition, \(p_{\ell,N+1}\) downward from \(p = 0\) at the surface, the interface mass flux by the layer fractions. MLVAM(n_layers=1, level=1) is the paper’s VAM system (4) with the constraints (5) row by row, inviscid and with \(e_x = 0\): the row of \(q_{1,1}\), \(r_{1,1}\) is \(-1\) times the paper’s, the constraint rows are \(-1\) times (5).

Options

setting
n_layers	number of layers
level	\(N\): \(u \in \mathbb{P}_N\), \(w, p \in \mathbb{P}_{N+1}\) per layer
dimension	2: (x, z). 3: (x, y, z), state q_x_l_k, q_y_l_k
closures	stress and interface pieces, for example [Newtonian(), NavierSlip(), StressFree()], plus UpwindInterface() for an upwinded interface velocity
stress	Stress.shallow the vertical shear, Stress.full the whole Newtonian tensor
vertical_traction	with Stress.full: the interface trace of the shear in the vertical rows (VerticalTraction.shared, one_sided, function)
parameters	g, rho, nu, lambda_s, e_x, layer fractions l_1, …, l_{n_layers-1}

The non-hydrostatic pressure makes the system a DAE: chorin_split gives the predictor, pressure and corrector systems of the Chorin solver.

from zoomy_core.model.models import MLVAM, Newtonian, NavierSlip, StressFree
from zoomy_core.systemmodel import SystemModel

mlvam = MLVAM(n_layers=2, level=1, dimension=2,
              closures=[Newtonian(), NavierSlip(), StressFree()])
model = SystemModel.from_model(mlvam)
split = mlvam.chorin_split(system_model=model)

In the GUI

There is no MLVAM card. Compose, Model card VAM, then in the ## Model section of case.py import MLVAM instead of VAM and replace VAM( by MLVAM(n_layers=2, . Run with the Split / Chorin (NumPy) solver, which takes the card’s split.