In the system (1)–(3), $h(x,t)$ represents the total water height at each point $x\in \mathrm{\Omega }\subset ℝ$, and time $t\ge 0$, where $\mathrm{\Omega }$ is the considered (horizontal) domain. The water height is decomposed along the vertical axis into a prescribed number of layers $L\ge 1$ (see 1). For any layer $\alpha$, its thickness will be assumed to be $h_{\alpha }=l_{\alpha }h$, for some values $l_{\alpha }\in (0,1)$ such that ${\sum <!--\; FUNCTION\; APPLICATION\; -->}_{\alpha =1}^L{l}_{\alpha }=1$. Usually, $l_{\alpha }=1∕L$ is selected, although any other choice is possible.

The upper and lower interfaces of layer $\alpha$ will be represented by $z_{\alpha +1∕2}$ and $z_{\alpha -1∕2}$, respectively, that is,

The uppermost interface corresponds to the sea surface, denoted by $\eta (x,t)=h(x,t)+z_b(x,t)$; the lowermost corresponds to the seafloor basin represented by $z_b(x,t)$, which is supposed to be perturbed by the earthquake; and finally,

denotes the level of the middle point of the layers. The depth-averaged velocities in the horizontal and vertical directions are written as $u_{\alpha }(x,t),$ and $w_{\alpha }(x,t)$, respectively. Finally, $p_{\alpha +1∕2}$ denotes the non-hydrostatic pressure at the interface $z_{\alpha +1∕2},$ which is assumed to be $0$ at the free surface. $ū={\sum <!--\; FUNCTION\; APPLICATION\; -->}_{\alpha =1}^L{l}_{\alpha }{u}_{\alpha }$ is the mean of the depth-averaged horizontal velocities.

Moreover, for any variable $f\in \{u,w\}$ of the system, we denote

As usual $g=9.81m∕s^2$ denotes the gravity acceleration and ${\mathrm{\Gamma }}_{\alpha +1∕2}$ parametrizes the mass transfer across interfaces

where we assume no mass transfer through the seafloor or the free surface (${\mathrm{\Gamma }}_{1∕2}={\mathrm{\Gamma }}_{L+1∕2}=0$). Each layer is supplemented with the following divergence-free constraint

where

and the term ${\partial }_t{z}_b$ accounts for the movement of the bottom interface.

Note that system (1)–(3) is endorsed with some extra dissipation models for accounting friction with the bottom (${\tau }_{\alpha }^u$), viscous terms that model the shear stresses between the layers ($K_{\alpha \pm 1∕2}$), and breaking of the waves near the coast (${\tau }_{\alpha }^w$). Here, we propose the following dissipation models (see [13]):

• For the friction effects between the water and the seafloor, we use a standard Gauckler-Manning friction formula that is applied at the lower layer

  $${\tau }_{\alpha }^u=\{\begin{array}{cc}\mathit{gL}{n}^2|u_1|{\displaystyle \frac{h{u}_1}{h^{4∕3}}},\hfill & \text{if }\alpha =1,\hfill \\ 0,\hfill & \text{if }\alpha \in \{2,\dots ,L\}.\hfill \end{array}$$

  (5)

• We follow a simplified model of the one presented in [1] for the shear stress between the layers

  $$K_{\alpha +1∕2}=-\nu {\displaystyle \frac{u_{\alpha +1}-u_{\alpha }}{(h_{\alpha +1}+h_{\alpha })∕2}},$$

  (6)

  where $\nu =5g\cdot 10^2\cdot n^2$ is a constant kinematic viscosity, and $K_{1∕2}=K_{L+1∕2}=0.$

• For the breaking dissipation model, we consider here an extension of the simple, efficient, and robust model considered in [6] for a two-layer model:

  $${\tau }_{\alpha }^w=C{w}_{\alpha }|\mathit{\partial x}(h{u}_{\alpha })|,\alpha \in \{1,\dots ,L\},$$

  (7)

  where the coefficient $C\equiv C(x,t)$ defines breaking criteria to switch on/off the dissipation of the energy due to the presence of a breaking wave (see [16]). Here, we shall use

  $$C=\{\begin{array}{cc}35\left({\displaystyle \frac{|ū|}{0.4\sqrt{\mathit{gh}}}}-1\right),\hfill & \text{if }|ū|>\sqrt{\mathit{gh}},\hfill \\ 0,\hfill & \text{if }|ū|\le \sqrt{\mathit{gh}}.\hfill \end{array}$$

Concerning the mathematical properties of the previous model, one can show that the system (1)–(4) satisfies an energy balance equation (see [10]) and that it has good dispersive properties. Indeed, using a standard Stokes-type Fourier analysis for the linearized previous system around the water at rest steady-states, phase, group velocities, and linear shoaling gradient are determined and compared with the Airy or Stokes linear theory for different numbers of layers (see 2, where relative errors are shown for the phase and group velocities, as well as for the shoaling gradient in percents). One can prove uniform convergence for the analytical values when the number of layers increases (see [10]).