We briefly describe the numerical discretization of system (1)–(4) assuming a constant split of the layers $l_{\alpha }=1∕L$. Here, we follow the procedure described in [5]. The numerical scheme is based on a two-step projection-correction method, similar to the standard Chorin’s projection method for Navier-Stokes equations ([4]). That is a standard procedure when dealing with dispersive systems (see [9, 5, 12, 15] and references therein).

First, we shall solve the non-conservative hydrostatic underlying system (1)–(3) given by the compact equation

$$\mathit{\partial t}U+\mathit{\partial x}F(U)+B(U)\mathit{\partial x}U+G(U)\mathit{\partial x\eta }=0,$$

(8)

where the following compact notation has been used:

$$U=\left(\begin{array}{c}\hfill h\hfill \\ \hfill h{u}_1\hfill \\ \hfill ⋮\hfill \\ \hfill h{u}_L\hfill \\ \hfill h{w}_1\hfill \\ \hfill ⋮\hfill \\ \hfill h{w}_L\hfill \end{array}\right),F(U)=\left(\begin{array}{c}\hfill \mathit{hu}\hfill \\ \hfill u_{1}h{u}_1\hfill \\ \hfill ⋮\hfill \\ \hfill u_{L}h{u}_L\hfill \\ \hfill u_{1}h{w}_1\hfill \\ \hfill ⋮\hfill \\ \hfill u_{L}h{w}_L\hfill \end{array}\right),G(U)=\left(\begin{array}{c}\hfill 0\hfill \\ \hfill \mathit{gh}\hfill \\ \hfill ⋮\hfill \\ \hfill \mathit{gh}\hfill \\ \hfill 0\hfill \\ \hfill ⋮\hfill \\ \hfill 0\hfill \end{array}\right),$$

and $B$ is a matrix such $B\mathit{\partial x}U$ contains the non-conservative products related to the mass transfer across interfaces appearing at the momentum equations.

Then, in a second step, the non-hydrostatic terms from the right-hand side of equations (2)-(3) collected by the pressure vector

$$T(h,\mathit{\partial xh},z_b,\mathit{\partial x}{z}_b,P,\mathit{\partial x}P)=-L\left(\begin{array}{c}\hfill 0\hfill \\ \hfill \frac{1}{2}{h}_1{\partial }_x\left(p_{3∕2}+p_{1∕2}\right)-\left(p_{3∕2}-p_{1∕2}\right){\partial }_x{z}_1\hfill \\ \hfill ⋮\hfill \\ \hfill \frac{1}{2}{h}_L{\partial }_x\left(0+p_{L-1∕2}\right)-\left(0-p_{L-1∕2}\right){\partial }_x{z}_L\hfill \\ \hfill h_1\left(p_{3∕2}-p_{1∕2}\right)\hfill \\ \hfill ⋮\hfill \\ \hfill h_L\left(0-p_{L-1∕2}\right)\hfill \end{array}\right),P=\left(\begin{array}{c}\hfill p_{1∕2}\hfill \\ \hfill p_{3∕2}\hfill \\ \hfill ⋮\hfill \\ \hfill p_{L-1∕2}\hfill \end{array}\right),$$

as well as the divergence constraints at each layer given in (4) will be taken into account. Concerning the constraints, we will equivalently impose

$$B_{\alpha }=0,\text{ where }{B}_1=I_1,B_2=I_2-I_1,\dots ,B_L=I_L-I_{L-1},$$

so that the divergence impositions read as follows

$$\begin{array}{cc}B\left(U^{(k)},\partial x{U}^{(k)},z_b,\mathit{\partial x}{z}_b\right)=& \\ & \\ \left(\begin{array}{c}\hfill -h_1{\partial }_x(h_1{u}_1)+2{\partial }_x{z}_1{h}_1{u}_1-2h_1{w}_1+2h_1{\partial }_t{z}_b\hfill \\ \hfill -h_2{\partial }_x(h_2{u}_2)-h_1{\partial }_x(h_1{u}_1)+2h_2{u}_2{\partial }_x{z}_2-2h_1{u}_1{\partial }_x{z}_1-2h_2{w}_2+2h_1{w}_1\hfill \\ \hfill ⋮\hfill \\ \hfill -h_L{\partial }_x(h_L{u}_L)-h_{L-1}{\partial }_x(h_{L-1}{u}_{L-1})+2h_L{u}_L{\partial }_x{z}_L-2h_{L-1}{u}_{L-1}{\partial }_x{z}_{L-1}-2h_L{w}_L+2h_{L-1}{w}_{L-1}\hfill \end{array}\right).& \end{array}$$

System (8) is discretized using a second-order finite volume PVM positive-preserving well-balanced path-conservative method [2]. As usual, we consider a set of $N$ finite volume cells $I_i=[x_{i-1∕2},x_{i+1∕2}]$ with constant lengths $\mathrm{\Delta }x$ and define

$$U_i(t)={\displaystyle \frac{1}{\mathrm{\Delta }x}}{\int \; <!--\; nolimits\; -->}_{I_i}U(x,t)\mathit{dx},$$

the cell average of the function $U(x,t)$ on cell $I_i$ at time $t.$ Concerning non-hydrostatic terms, we consider mid-points $x_i$ of each cell $I_i$ and denote the point values of the function $P$ at time $t$ by

$$P_i(t)=\left(\begin{array}{c}\hfill p_{1∕2}(x_i,t)\hfill \\ \hfill p_{3∕2}(x_i,t)\hfill \\ \hfill ⋮\hfill \\ \hfill p_{L-1∕2}(x_i,t)\hfill \end{array}\right).$$

Non-hydrostatic terms will be approximated by second-order compact finite differences.

Let us detail the time stepping procedure followed. Assume given time steps $\mathrm{\Delta }{t}^n,$ and denote $t^n={\sum <!--\; FUNCTION\; APPLICATION\; -->}_{k\le n}\mathrm{\Delta }{t}^k.$ To obtain second-order accuracy in time, the two-stage second-order TVD Runge-Kutta scheme is adopted. At the $k$th stage, $k\in \{1,2\},$ the two-step projection-correction method is given by

$$\begin{array}{rcll}\frac{U^{(\tilde{k})}-U^{(k-1)}}{\mathrm{\Delta }t}+\mathit{\partial x}F(U^{(k-1)})+B(U^{(k-1)})\partial x{U}^{(k-1)}+G(U^{(k-1)})\partial x{z}_b& =& 0,& \text{(9)}\\ \frac{U^{(k)}-U^{(\tilde{k})}}{\mathrm{\Delta }t}-T(h^{(k)},\partial x{h}^{(k)},z_b,\mathit{\partial x}{z}_b,P^{(k)},\partial x{P}^{(k)})& =& 0& \text{(10)}\\ B(U^{(k)},\partial x{U}^{(k)},z_b,\mathit{\partial x}{z}_b)& =& 0& \text{(11)}\end{array}$$

where $U^{(0)}$ is $U$ at the time level $t^n$, $U^{(\tilde{k})}$ is an intermediate value in the two-step projection-correction method that contains the numerical solution of the hyperbolic system (8) at the corresponding $k$th stage of the Runge-Kutta, and $U^{(k)}$ is the $k$th stage estimate. After that, a final value of the solution at the $t^{n+1}$ time level is obtained:

$$U^{n+1}=\frac{1}{2}{U}^n+\frac{1}{2}{U}^{(2)}.$$

Observe that equations (10-11) require, at each stage of the calculation respectively, to solve a Poisson-like equation for each one of the variables contained in $P^{(k)}.$ The resulting linear system is solved using an iterative Jacobi method combined with a scheduled relaxation (see [5, 9]). Note that the usual CFL restriction must be imposed for the computation of the time step $\mathrm{\Delta }t.$

When friction with the bottom (5), viscous shear stress (6), and the breaking model (7) are considered, they will be computed semi-implicitly at the end of the second step of the projection-correction method at each $k$th stage of the TVD Runge-Kutta method as it is done in [6]. Note that the resolution of a straightforward tridiagonal system on the vertical for each volume $I_i$ is exclusively required for the viscosity model. In contrast, the friction and breaking models can be considered by solving conventional algebraic problems, as is commonly the practice for friction models.

The resulting numerical scheme is well-balanced for the water at rest solution and is linearly $L^{\infty }$-stable under the usual CFL condition related to the hydrostatic system. It is also worth mentioning that the numerical scheme is positive preserving and can deal with emerging topographies.