## 2 Numerical Model
### Numerical Flux Computation  
Among the different numeric fluxes developed in the past, [39] proposed the HLL scheme (Harten-Lax-Leer) [40]) as the most stable method and the one best able to compute fluxes F and G for this case. They suggested the following HLL fluxes at the cell interfaces.  
* The computation of the F flux: \[\text{F}(\text{U}_{\text{L}},\text{U}_{\text{R}})=\left\{\begin{matrix}\text{F} (\text{U}_{\text{L}})&\text{if }0\leq\text{C}_{1}^{\text{x}}\\ \frac{\text{C}_{2}^{\text{x}}\text{F}(\text{U}_{\text{L}})-\text{C}_{1}^{ \text{x}}\text{F}(\text{U}_{\text{R}})+\text{C}_{1}^{\text{x}}\text{C}_{2}^ {\text{x}}(\text{U}_{\text{R}}-\text{U}_{\text{L}})}{\text{C}_{2}^{\text{x} }-\text{C}_{1}^{\text{x}}}&\text{if }\text{C}_{1}^{\text{x}}<0<\text{C}_{2}^{\text{x}}\\ \text{F}(\text{U}_{\text{R}})&\text{if }\text{C}_{2}^{\text{x}}\leq 0\end{matrix}\right.\] (26) where \(\text{U}_{\text{L}}\) and \(\text{U}_{\text{R}}\) are, respectively, the left and right approximations of the solution in the cell. The coefficients \(\text{C}_{1}^{\text{x}}\) at \(\text{C}_{2}^{\text{x}}\) are calculated as: \[\text{C}_{1}^{\text{x}}=\min_{\text{U}=\text{U}_{\text{L}}\cup\text{U}_{ \text{R}}}(\min_{\text{j}=1,2,3}\beta_{1}^{\text{x}}(\text{U}));\,\text{C}_{2 }^{\text{x}}=\max_{\text{U}=\text{U}_{\text{L}},\text{U}_{\text{R}}}(\max_{ \text{j}=1,2,3}\beta_{1}^{\text{x}}(\text{U}))\] (27) where the term \(\beta_{1}^{\text{x}}\), with j = 1, 2, 3, represents the eigenvalues of the Jacobian of F and can be expressed as: \[\beta_{1}^{\text{x}}=\text{u}-\sqrt{\text{gh}},\,\,\beta_{2}^{\text{x}}=\text {u},\,\,\beta_{3}^{\text{x}}=\text{u}+\sqrt{\text{gh}}\,.\] (28)
* The computation of the G flux: \[\text{G}(\text{U}_{\text{D}},\text{U}_{\text{U}})=\left\{\begin{matrix}\text{G} (\text{U}_{\text{D}})&\text{if }0\leq\text{C}_{1}^{\text{y}}\\ \frac{\text{C}_{1}^{\text{y}}\text{C}_{2}^{\text{y}}\text{G}(\text{U}_{\text {D}})-\text{C}_{1}^{\text{y}}\text{G}(\text{U}_{\text{U}})+\text{C}_{1}^{ \text{y}}\text{C}_{2}^{\text{y}}(\text{U}_{\text{U}}-\text{U}_{\text{D}})}{ \text{C}_{2}^{\text{x}}-\text{C}_{1}^{\text{y}}}&\text{if }\text{C}_{1}^{\text{y}}<0<\text{C}_{2}^{ \text{y}}\\ \text{G}(\text{U}_{\text{U}})&\text{if }\text{C}_{2}^{\text{y}}\leq 0\end{matrix}\right.\] (29) where \(\text{U}_{\text{D}}\) and \(\text{U}_{\text{U}}\) are, respectively, the top and bottom approximations of the solution in the cell. The coefficients \(\text{C}_{1}^{\text{y}}\) at \(\text{C}_{2}^{\text{y}}\) are calculated as: \[\text{C}_{1}^{\text{y}}=\min_{\text{U}=\text{U}_{\text{D}},\text{U}_{\text{U}} }(\min_{\text{j}=1,2,3}\beta_{1}^{\text{y}}(\text{U}));\,\text{C}_{2}^{\text{y} }=\max_{\text{U}=\text{U}_{\text{D}},\text{U}_{\text{U}}}(\max_{\text{j}=1,2,3 }\beta_{1}^{\text{y}}(\text{U}))\] (30) where the term \(\beta_{1}^{\text{y}}\) with j = 1, 2, 3 is the eigenvalues of the Jacobian of G and can be expressed as: \[\beta_{1}^{\text{y}}=\text{v}-\sqrt{\text{gh}},\,\,\beta_{2}^{\text{y}}=\text {v}\,,\,\,\beta_{3}^{\text{y}}=\text{v}+\sqrt{\text{gh}}\,.\] (31)Here, h is the water height; g is the gravitational acceleration; u and v are the depth-averaged velocity components in the \(x\) and \(y\) directions.