Re-evaluating efficiency of first-order numerical schemes for two-layer shallow water systems by considering different eigenvalue solutions
RRe-evaluating efficiency of first-order numerical schemesfor two-layer shallow water systems by consideringdifferent eigenvalue solutions
Nino Krvavica a, ∗ a University of Rijeka, Faculty of Civil Engineering, Radmile Matejcic 3, 51000 Rijeka,Croatia
Abstract
The efficiency of several first-order numerical schemes for two-layer shallowwater equations (SWE) are evaluated here by considering different eigenvaluesolutions. This study is a continuation of our previous work [25] in which wehave proposed an efficient implementation of a Roe solver for two-layer SWEbased on analytical expressions for eigenvalues and eigenvectors. In this work,the accuracy and computational cost of numerical, analytical, and approximatedeigenvalue solvers are compared when implemented in Roe, Intermediate FieldCaPturing (IFCP) and Polynomial Viscosity Matrix (PVM) schemes. Severalnumerical tests are performed to examine the overall efficiency of numericalschemes with different eigenvalue solvers when computing two-layer shallow-water flows. The results confirm that analytical eigenvalue solutions are muchfaster than numerical solvers, with a computational cost closer to approximateexpressions. Consequently, the Roe scheme with analytical solutions to theeigenstructure is equally efficient as the IFCP scheme. On the other hand,IFCP and PVM schemes with analytical solutions to eigenvalues are found to beequally efficient as those with approximated expressions. Analytical eigenvaluesshow slightly better results when dealing with larger density differences betweenthe layers. ∗ Corresponding author
Email address: [email protected] (Nino Krvavica)
Preprint submitted to Elsevier January 15, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J a n eywords: shallow-water equations, two-layer flow, eigenvalues, Roe scheme, IFCPscheme, PVM scheme, finite volume method
1. Introduction
Two-layer shallow-water equations (SWE) are widely used to simulate geo-physical flows in stratified conditions. Some examples of a two-layer configu-ration include exchange flows in sea straits [15, 14], highly stratified estuaries[23, 24], as well as various types of gravity currents [27, 3], such as mudflows[7], debris flows [35, 30], submarine avalanches [20, 29, 36], and pyroclastic flows[19]. Although such processes can be described more accurately by 3D Navier-Stokes equations, two-layer shallow water models make a popular alternativebecause of their simplicity and a significantly lower computational cost.Two-layer SWEs are mathematically defined as hyperbolic systems of cou-pled conservation laws with source terms, or so-called balance laws [15]. Theseequations are challenging to solve numerically because of the layer coupling andnon-conservative source terms accounting for the variable geometry, friction, orentrainment. Over the last two decades, a numerical resolution of two-layerSWE has been an object of intense research [15, 14, 26, 16, 31, 6, 21, 28, 8, 18].A popular choice for numerical resolution of two-layer SWEs are finite vol-ume methods (FVM), and among them a family of path-conserving schemes[39, 5, 15, 34, 33].Implementation of the path-conservative schemes involves a numerical vis-cosity matrix, which is usually derived from some or all eigenvalues of a cor-responding Jacobian of the flux matrix. The choice of the numerical viscositymatrix determines the numerical diffusion and accuracy of the scheme. Since an-alytical expressions for eigenvalues of two-layer SWE systems were consideredunavailable until recently [14, 35, 1, 21, 37], either approximate expressions[38, 1] or numerical algorithms were used instead. Unfortunately, numerical al-gorithms, such as root-finding and eigensolver methods, make schemes computa-2ionally (too) demanding. For example, several studies evaluating the efficiencyof first-order schemes in solving two-layer flows [16, 10, 21], find Roe schemes,which are based on all eigenvalues, to be the most accurate, but far less efficientthen some other first-order schemes, mainly because of excessive computationalcosts when performing a full spectral decomposition by a numerical eigensolver.Recently, a new solution [25] to the efficiency problem for the Roe schemewas proposed. This new approach introduced a semi-analytical implementationof the Roe scheme based on simple closed-form solutions for the eigenvalues andeigenvectors [25]. New scheme, named A-Roe, was found to be much faster thanthe numerical implementation of the Roe scheme, while producing equally ac-curate results [25]. An additional advantage is that closed-form solutions enablea direct and accurate prediction of complex eigenvalues and implementation ofa corrective algorithm for the loss of hyperbolicity [25].Our previous paper [25] evaluated the performance of the A-Roe schemebased on analytical expressions for eigenvalues and eigenvectors, and comparedits efficiency to Lax-Friedrichs (LF) and GFORCE scheme, as well as the Poly-nomial Viscosity Matrix redefinition of the Roe scheme (PVM-Roe) [10] andIntermediate Field CaPturing scheme (IFCP) [21]. The former two do notuse any eigenvalue information, whereas the latter two use all four eigenvalues[16, 10, 21]. In the case of PVM-Roe and IFCP schemes, the classical imple-mentation based on approximated eigenvalues was considered. Furthermore,our previous paper [25] focused on overall advantages of using analytical eigen-values; primarily, increased computational speed and hyperbolicity correctionalgorithm.The present study is a natural continuation of the previous paper [25]. Oneevident question emerged from the main conclusions of the previous paper - ifanalytical expressions for the eigenstructure significantly improved the efficiencyof the Roe scheme, what effect will they have on the performance of other first-order schemes that are also based on some or all eigenvalues? Clearly, theefficiency of those numerical schemes should be re-evaluated.The main aim of this study is to investigated the sensitivity of several other3umerical schemes, namely PVM-2U [10] and IFCP [21], on the choice of eigen-values and to evaluate potential benefits of using analytical eigenvalues insteadof recommended approximated expressions. For this purpose, the accuracy andcomputational speed of recently proposed closed-form eigenvalue solutions arecarefully compared against two available alternatives - the numerical eigen-solvers and approximated expressions for eigenvalues. Next, the sensitivity ofnumerical schemes to the choice of an eigenvalue solver is assessed. And finally,the overall computational efficiency of Roe, IFCP, and PVM-2U schemes withdifferent eigenvalues is evaluated by performing several numerical tests whichconsider different density differences between the layers and different channel ge-ometries. Some significant remarks on the implementation of numerical schemes,which increase their efficiency, are also presented.
2. Two layer shallow-water flow: Theory, eigenvalues, and numericalschemes
A one-dimensional (1D) two-layer shallow-water flow in prismatic channelswith rectangular cross-sections of constant width is considered for all tests. Thegoverning system of equations written in a general vector form is repeated herein a more compact form for context and reproducibility [15]: ∂ w ∂t + ∂ f ( w ) ∂x = B ( w ) ∂ w ∂x + g ( w ) , (1)where x refers to the axis of the channel and t is time. The vector of conservedquantities w is defined as: w = (cid:110) h q h q (cid:111) T , (2)where h j is the layer thickness (or depth), q j = h j u j is the layer flow rate perunit width, and index j = 1 , f(w) is: f ( w ) = (cid:110) q q h + g h q q h + g h (cid:111) T , (3)4here g is acceleration of gravity. Matrix B(w) is a result of coupling thetwo-layer system, defined as [15]: B ( w ) = − c
00 0 0 0 − rc , (4)where r = ρ /ρ < ρ and thelower layer density ρ , and c j = gh j is propagation celerity of internal andexternal perturbations (waves), for j = 1 ,
2. Finally, the bathymetry sourceterm g(w) is defined as follows [15]: g ( w ) = (cid:110) − gh b d x − gh b d x (cid:111) T , (5)where b is the bed elevation. The system given by Eq. (1) can be rewritten in the following quasi-linearform [15]: ∂ w ∂t + A ( w ) ∂ w ∂x = g ( w ) , (6)where A (w ) = ∂ f ( w ) ∂ w − B ( w ) = J ( w ) − B ( w ) (7)is the pseudo-Jacobian matrix that contains the flux gradient terms as well asthe coupling terms: A (w) = c − u u c
00 0 0 1 rc c − u u . (8)The four eigenvalues of A (w) define the propagation speeds of barotropic(external) and baroclinic (internal) perturbations. In most geophysical flows,one of two external eigenvalues is negative λ − ext <
0, while the other is positive λ + ext > .2.1. Numerical eigenvalues Eigenvalues and eigenvectors of A (w) can be numerically obtained by solv-ing the following equation: A K = KΛ . (9)where Λ is a 4 × λ k , k =1 , ..,
4, and K is matrix whose columns are the corresponding right eigenvectors.Usually a QR algorithm is used for this purpose [4]. The following approximation derived under the assumption of r ≈ u ≈ u are usually used for computing the internal and external eigenvalues[38]: λ ± ext = U ± (cid:112) g ( h + h ) (10) λ ± int = U ± (cid:115) g (1 − r ) h h h + h (cid:20) − ( u − u ) g (1 − r )( h + h ) (cid:21) , (11)with U = h u + h u h + h and U = h u + h u h + h . (12)Note that Eqs. (10) and (11) are valid only when dealing with two layersof similar densities ( r = ρ /ρ ≈
1) and when velocities in both layers arecomparable ( u ≈ u ). Those conditions are found in some stratified flows innature, such as exchange flows through sea straits [14, 17] or some cases ofhighly stratified estuaries [24]. However, for geophysical flows characterized bya larger relative density difference, such as granular, debris or mud flows, theapproximated values may significantly deviate from exact values and cannotaccurately predict a possible hyperbolicity loss [1, 37, 25]. Recently, a simple closed-form approach for computing real roots of thecharacteristic quartic Eq. (13) was proposed [25]. The solutions are based onFerrari’s formulas [2], they consist of eight simple evaluations, and are repeatedhere for consistency and clarity. A detailed derivation of these equations is6vailable in [25]. The proposed closed-form solutions to eigenvalues are given interms of coefficients a, b, c and d of a characteristic polynomial of matrix A (w) [25]: p ( λ ) = λ + aλ + bλ + cλ + d (13)with: a = − u + u ) , (14) b = u − c + 4 u u + u − c , (15) c = − u (cid:0) u − c (cid:1) − u (cid:0) u − c (cid:1) , (16) d = (cid:0) u − c (cid:1) (cid:0) u − c (cid:1) − rc c . (17)Real eigenvalues are then computed by the following expressions [25]: λ ± ext = λ , = − a ± √ Z + (cid:113) − A − Z ∓ B √ Z , (18) λ ± int = λ , = − a ± √ Z − (cid:113) − A − Z ∓ B √ Z . (19)where Z = 13 (cid:18) (cid:112) ∆ cos φ − A (cid:19) , (20) φ = arccos (cid:18) ∆ √ ∆ (cid:19) , (21)with A = 2 b − a , (22) B = 2 c − ab + a . (23)and ∆ = b + 12 d − ac, (24)∆ = 27 a d − abc + 2 b − bd + 27 c . (25)Note, that it is possible to combine these equations into a single explicit so-lution in terms of conserved variables, but the resulting formula would certainlybe too extensive to be presented in a journal format, and probably not opti-mized to be implemented in a computational algorithm. However, an exampleof such expanded formulation is available [40].7 .2.4. Nature of eigenvalues Since the relative density difference r has a major influence on internal eigen-values, they are usually smaller than the external ones. Therefore, the followingindexing and order of eigenvalues will be used herein: λ = λ − ext , λ = λ − int , λ = λ + int , λ = λ + ext (26)It is worth mentioning that the external eigenvalues λ ± ext are always real[32, 1]; however, at sufficiently large relative velocities ∆ u = | u − u | , as wellas for very low or very high relative densities r , the internal eigenvalues λ ± int may become complex and the governing system loses its hyperbolic character[32, 12].Although analytical expressions given by Eqs. (18) and (19) are valid onlyfor real roots, a simple hyperbolicity correction algorithm can be incorporated,which verifies whether the term in Eq. (21) satisfies the condition | ∆ √ ∆ | < A class of path-conservative schemes are considered here to approximate thegoverning equations for two-layer shallow water flow [33]. A first order accuratepath-conservative scheme for Eq. (6) may be written as follows [33]: w n +1 i = w ni − ∆ t ∆ x (cid:16) D + i − / + D − i +1 / (cid:17) (27)where ∆ x and ∆ t are the respective spatial and time increment (consideredconstant here for simplicity), w ni denotes the approximate cell-averaged values ofthe exact solution obtained by the numerical scheme at cell I i = [ x i − / , x i +1 / ]in time t n = n ∆ t , and matrices D ± i +1 / are continuous functions of conservedvariables D ± (cid:0) w ni , w ni +1 (cid:1) .For the governing system of equations, a generalized numerical scheme based8n Roe linearisation [39] may be written by Eq. (27), with: D ± i +1 / = 12 (cid:104) f ( w ni +1 ) − f ( w ni ) + B i +1 / ( w ni +1 − w ni ) − g i +1 / ± Q i +1 / (cid:16) w ni +1 − w ni − A − i +1 / g i +1 / (cid:17) (cid:105) (28)where Q i +1 / is a numerical viscosity matrix that determines the numericaldiffusivity of the results, and whose choice depends on a particular scheme [21].Matrices and vectors B i +1 / , A i +1 / , and g i +1 / correspond to B ( w ni +1 , w ni ), A ( w ni +1 , w ni ), and g ( w ni +1 , w ni ), respectively, evaluated at the cell interface aftera suitable Roe linearization is performed (see [15] and [25] for details).The viscosity matrix in Roe schemes coincides with the absolute pseudo-Jacobian matrix [15]: Q i +1 / = | A i +1 / | (29)and the absolute value of A i +1 / can be directly obtained from: | A i +1 / | = K i +1 / | Λ i +1 / | K − i +1 / . (30)where | Λ i +1 / | is a 4 × | λ k | , k = 1 , .., A i +1 / . The latter iscomputational expensive, while the former is less accurate. In both cases, anadditional re-composition of the viscosity matrix from eigenstructure is required(see Eq. (30)), which imposes an additional computational cost. This drawbackhas motivated a development of numerical schemes in which the viscosity matrix Q i +1 / is directly approximated from pseudo-Jacobian matrix.One possible way to construct such a scheme is to approximate the viscositymatrix by a polynomial function. Those class of methods are called PolynomialViscosity Matrix (PVM) [10]. PVM methods define the viscosity matrix througha general polynomial evaluation of the pseudo-Jacobian matrix, given by: Q i +1 / = p l (cid:0) A i +1 / (cid:1) (31)9here p l ( x ) is a polynomial of degree lp l ( x ) = l (cid:88) j =0 α j x j (32)The main advantage of PVM schemes is that a full spectral decomposition isnot required.In [10] several different PVM schemes were presented and tested, some ofthem were derived from other popular Riemann solvers redefined under thePVM formulation, and some of them were newly proposed. In our previouspaper [25] we chose a PVM redefinition of the Roe scheme (PVM-Roe) thatis based on all four eigenvalues, and here we choose a scheme called PVM-2Uthat is based only on two external eigenvalues. In [10], PVM-2U proved to bethe most efficient throughout their numerical tests and is, therefore, selectedhere for further analysis and comparison regarding the effects of the eigenvaluesolvers.The viscosity matrix of PVM-2U may be computed as [10]: Q i +1 / = α Id + α A i +1 / + α A i +1 / , (33)where coefficients α k , k = 0 , , α k , k = 0 , , t ∆ x max | λ j,i +1 / | = γ ≤ There are some modifications of numerical schemes that can be made tooptimize the algorithms and improve their computational performance.The first optimization deals with the Roe scheme and re-composition of theviscosity matrix given by Eq. (30). Although analytical closed-form solutionsto K − are available in the Appendix of [25], it is computationally faster torewrite Eq. (30) as K Ti +1 / | A i +1 / | T = ( K i +1 / | Λ i +1 / | ) T , (35)which corresponds to a general matrix equation Ax = B , solve it numerically for | A i +1 / | T (for example, by a LAPACK routine gesv [4]), and then transpose it.This is about 2-3 times faster than finding the inverse of K i +1 / and performingmatrix multiplication to obtain the viscosity matrix as written in Eq. (30).A second optimization is available for the family of PVM schemes. UsingEq. (7) and the usual Roe linearization of the pseudo-Jacobian matrix: J ( w i +1 , w i ) · ( w i +1 − w i ) = f ( w i +1 ) − f ( w i ) (36)term Q i +1 / (cid:0) w ni +1 − w ni (cid:1) given in Eq. (28) may be replaced by: C i +1 / (cid:2) f ( w ni +1 ) − f ( w ni ) − B i +1 / ( w ni +1 − w ni ) (cid:3) (37)where C i +1 / = Q i +1 / A − i +1 / lowers the order of a viscosity matrix by one,and is already needed for the source term discretization. In other words, fullviscosity matrix Q i +1 / is not required for the family of PVM schemes; instead,only C i +1 / is computed. For both PVM-2U and IFCP, C i +1 / is defined as: C i +1 / = α A − i +1 / + α Id + α A i +1 / , (38)11n this way, computation of the square of A i +1 / (see Eq. (33)) is avoided.
3. Results
Several numerical tests are presented to evaluate the efficiency of Roe, IFCP,and PVM-2U schemes with different eigenvalue solvers. First, the accuracyand computational speed of numerical, analytical and approximated eigenvaluesolvers is analysed. Next, the performance of numerical schemes in computingthe numerical viscosity matrix is examined, as well as their sensitivity to thechoice of eigenvalues. Finally, five numerical tests are given to analyse the overallefficiency of different numerical schemes in computing two-layer shallow-waterflows.In particular, three eigensolver algorithms are examined: • N-Eig uses a numerical eigenvalue solver which decomposes a generalsquare matrix into a diagonal matrix Λ whose elements are eigenvalues,and matrix K whose columns are right eigenvectors. This algorithm isimplemented in Python using Numpy function numpy.linalg.eig whichis based on the geev LAPACK routines written in FORTRAN [4]. • A-Eig is an analytical eigenvalue solver based on a closed-form solution tothe roots of the characteristic quartic given by Eqs. (18) and (19) • E-Eig only estimates eigenvalues based on the approximations given byEqs. (10) and (11)By combining different eigensolvers with Roe, IFCP, and PVM-2U schemes,the following numerical algorithms for computing the viscosity matrix and solv-ing two-layer SWEs are chosen for the efficiency analysis: N-Roe, A-Roe, E-Roe,A-IFCP, E-IFCP, A-PVM2, and E-PVM2 schemes.The accuracy of algorithms and schemes are evaluated by using either ab-solute error AE, relative bias error RBE, or root relative square error RRSE,defined respectively as: AE Φ = | Φ − Φ ref | (39)12 BE Φ = Φ − Φ ref Φ ref (40) RRSE Φ = (cid:113)(cid:80) Mn =1 [Φ( x n , t end ) − Φ ref ( x n , t end )] (cid:113)(cid:80) Mn =1 Φ ref ( x n , t end ) , (41)where Φ is the evaluated parameter (e.g. eigenvalue λ , numerical flux f , depth h or velocity u ), Φ ref are the corresponding reference values, and M is a numberof spatial points.All numerical algorithms have been implemented in Python 3.6 and vec-torized using the Numpy package. The tests have been performed on 64-bitWindows 10 machine with Intel Core i7-3770 3.4 GHz processor. All algorithmsfor computing the eigenstructure and viscosity matrix are freely available onGithub [22]. This test examines the accuracy and computational speed (CPU time) of cor-responding three different algorithms for computing eigenvalues N-Eig, A-Eig,and E-Eig. Although eigenvalues computed by the numerical eigensolver are not exact , their errors are of the order of round-off errors. Therefore, the numericalresults are used as a reference when evaluating the accuracy of analytical andapproximated eigenvalues.Since the main idea is to evaluate algorithms for solving eigenvalues as aintegral part of Roe, IFCP and PVM-2U schemes for two-layer SWEs, physicallyrealistic flow parameters (which always produce real eigenvalues) are chosen forthis test. Therefore, a large set of parameters ( N = 10 ) is randomly generatedfrom a given range: 1 . < h , < . − . < u , < . − . Differentdensity ratios are selected, namely r = 0 . , . , . , . , and 0 .
3, with g set to9 . − .Fig. 1 illustrates the statistical representation of the absolute errors of A-Eigand E-Eig algorithms computed by Eq. (39) for one million sets of independentflow parameters. The average error of the analytic solver is around 10 − andthe maximum errors are below 10 − which is close to a round-off error of13he numerical eigensolver. In other words, both these solvers produce almostidentical results. The average errors of the approximated expressions are morenoticeable, and they increase with density difference. For r = 0 .
98, averageerrors are around 10 − , whereas for r = 0 .
3, average errors are larger than10 − . This is expected, since eigenvalue approximations are derived under theassumption of r ≈ u ≈ u , and become less accurate as parametersdeviate from these assumptions.To examine the lack of accuracy of the eigenvalues approximations, Fig. 2illustrates the relative bias errors of individual eigenvalues for small r . It seemsthat the eigenvalue approximations always overestimate external eigenvalues,and underestimate the internal ones. Contrary to approximated expressionsgiven by Eq. 10 that suggest how r only affects the internal eigenvalues, it isclear from this analysis that r equally affects external and internal eigenvalues. r a b s o l u t e e rr o r A) analytical eigenvalues (A-Eig) r a b s o l u t e e rr o r B) eigenvalue approximations (E-Eig)
Figure 1: Boxplot of the absolute errors when computing eigenvalues by: A) the analyticalclosed-form solutions (A-Eig) and B) eigenvalue approximations (E-Eig). Boxes denote theinterquartile range and median value, while whiskers denote min and max values.
Table 1 shows the computational speed of three eigenvalue algorithms. Thefastest algorithm is the E-Eig (based on approximated eigenvalues) that needed0.11 s, followed by the analytical solver A-Eig with 0.29 s, and finally the numer-ical solver N-Eig with 4.14 s. Notice that both approximated and closed-formeigenvalue solvers are one order of magnitude faster than the numerical solver.These results are in agreement with [25] and suggest that the prevailing opinion14 xt int + int + ext eigenvalues0.200.150.100.050.000.050.100.150.20 r e l a t i v e b i a s e rr o r A) approximated eigenvalue errors r=0.5 ext int + int + ext eigenvalues0.200.150.100.050.000.050.100.150.20 r e l a t i v e b i a s e rr o r B) approximated eigenvalue errors r=0.3
Figure 2: Boxplot of the relative bias errors when using approximated eigenvalues E-Eig for:A) r = 0 . r = 0 .
3. Boxes denote the interquartile range and median value, whilewhiskers denote min and max values. in the scientific community about the ”computational complexity” of the ana-lytical solver is not justified. Although A-Eig needs double the time of E-Eig,it is considerably more accurate, as shown in Fig. 1.
Table 1: CPU times when computing one million sets of eigenvalues by numerical N-Eig,analytical A-Eig, and approximated E-Eig eigenvalue solvers (best of 5 runs).
Eigensolver N-Eig A-Eig E-EigCPU time (ms) 4136.5 286.8 109.7
The next test examines the accuracy and the computational speed of differ-ent implementations of Roe, IFCP, and PVM-2U schemes. The implementationdiffers in the choice of a particular eigenvalue algorithm, presented in the pre-vious subsection. The main goal is to investigate how do eigenvalue algorithmsaffect the efficiency of Roe, IFCP, and PVM-2U schemes when computing thecorresponding numerical viscosity matrix, and to analyse their sensitivity to aparticular choice of the eigenvalue solver.Same as in the previous example, physically realistic flow parameters arechosen. One million sets of parameters are randomly generated from a given15ange: 1 . < h , < . − . < u , < . − . Different density ratiosare selected, namely r = 0 . , . , . , . , and 0 .
3, with g set to 9 . − .Note that Roe schemes compute Q , whereas IFCP and PVM-2U schemescompute only C , which is additionally multiplied by A for comparison. Theerrors are estimated by a relative bias error given by Eq. (40), with N-Roebeing the reference value.Fig. 3 illustrates the statistical representation of relative bias errors forone million sets of independent flow parameters. Results obtained by A-Roe(Fig. 3a) and E-Roe (Fig. 3b) schemes correspond to the accuracy of the eigen-value solvers they use. The average error of the A-Roe method is of the order10 − and the maximum errors are always below 10 − . On the other hand, theaverage errors of the E-Roe scheme are much higher and increase with densitydifference. They amount to 10 − for r = 0 .
98 and just under 1 . r = 0 . − , whereas for E-IFCP the relative error grows with decreasing r . Itis interesting to note that E-IFCP, in contrast to E-Roe, underestimates theviscosity matrix. In PVM-2U scheme, the differences between the analyticaland approximated eigenvalues are negligible, although the errors generated byapproximated eigenvalues are somewhat higher for smaller values of r .It is interesting to note that the errors for r < . .98 0.9 0.7 0.5 0.3density ratio r r e l a t i v e b i a s e rr o r
1e 15 a) A-Roe method 0.98 0.9 0.7 0.5 0.3density ratio r r e l a t i v e b i a s e rr o r b) E-Roe method r r e l a t i v e b i a s e rr o r c) A-IFCP method 0.98 0.9 0.7 0.5 0.3density ratio r r e l a t i v e b i a s e rr o r d) E-IFCP method r r e l a t i v e b i a s e rr o r e) A-PVM2 method 0.98 0.9 0.7 0.5 0.3density ratio r r e l a t i v e b i a s e rr o r f) E-PVM2 method Figure 3: Boxplot of relative bias errors in computing one million numerical viscosity matricesby: a) A-Roe, b) E-Roe, c) A-IFCP, d) E-IFCP, e) A-PVM2, f) E-PVM2. Boxes denote theinterquartile range and median value, while whiskers denote min and max values.
Since all considered schemes (N-Roe, A-Roe, E-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2) are theoretically stable only for exact eigenvalues [15, 21],the well-balance properties of these schemes with differently computed eigen-values should be verified in practice. The first numerical test is designed to testthe well-balanced properties of the schemes when simulating water at rest in a18 able 2: Total CPU time (including time required only for eigenvalues Λ and eigenvectors K )when computing one million numerical viscosity matrices by N-Roe, A-Roe, E-Roe, A-IFCP,E-IFCP, A-PVM2, and E-PVM2. (best of 5 runs) Scheme N-Roe A-Roe E-Roe A-IFCP E-IFCP A-PVM2 E-PVM2 Λ (s) 4.13 0.29 0.11 0.29 0.11 0.29 0.11 K (s) 3.38 0.28 0.27 - - - -Total (s) 9.50 2.11 1.93 1.58 1.41 1.57 1.38 channel with smooth non-flat bed. To do so, the spatial domain is set to [0, 1],and a bed function is defined by a single bump: b ( x ) = cos ( π ( x − . / .
1) + 14 m , if 0 . < x < . . , otherwise (42)The initial condition is given by: h ( x,
0) = 0 . , h ( x,
0) = 0 . − b ( x ) (43) u ( x,
0) = u ( x,
0) = 0 m s − (44)Non-reflective conditions are imposed at the boundaries, and the relativedensity ratio is set to r = 0 .
98. A single grid density of ∆ x = 1/100 m and afixed time step ∆ t = 0 .
002 s is chosen.Figure 4 shows the interface profile and the lower layer velocities obtainedby N-Roe, A-Roe, E-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 numericalschemes at t = 0 . x = 1 /
100 m. All schemes satisfy the well-balanceproperties for water at rest (the velocities are close to the order of computationalprecision), except E-Roe, where noticeable spurious oscillations appear at 0 . 1e 15 N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 Figure 4: Test I: Computed interface and lower layer velocity obtained by N-Roe, A-Roe,E-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 scheme at t = 0 . x = 1 / 100 m r =0.98 This test investigates the performance of numerical schemes when simulatinga two-layer flow through a rectangular channel with a flat bottom topography.The flow structure is defined by two internal dam-breaks which eventually collideand produce a superimposed wave. A small density difference between the layersis considered, namely r = 0 . 98. The spatial domain is set to [0, 100], and theinitial condition is given by: h ( x, 0) = . , if 40 < x < 60 m0 . , otherwise h ( x, 0) = 1 . m − h ( x, 0) (45) u ( x, 0) = u ( x, 0) = 0 m s − (46)Non-reflective conditions are imposed at the boundaries. Several mesh sizesare considered, namely ∆ x = 1, 1/2, 1/4, and 1/8 m. A variable time step∆ t is evaluated at each step to satisfy CF L = 0 . 9. The reference solution iscomputed using the A-Roe scheme and a dense grid with ∆ x = 1/16 m.Figure 5 shows the temporal evolution of the interface and free-surface pro-files for the reference solution. Solutions obtained by other methods are not20llustrated, because they are almost undistinguishable at this scale. Therefore,a detail of solutions for the interface depth and lower layer velocity are shownin Fig. 6 where N-Roe, A-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 nu-merical schemes are compared against the reference solution at t = 25 s with∆ x = 1 / r = 0 . d e p t h ( m ) Figure 5: Test II: Temporal evolution of the interface and surface profile (the reference solu-tion) A CPU time vs. relative root square error E Φ is presented in Fig. 7. Theresults suggest that the E-IFCP is the most efficient scheme, closely followedby A-IFCP and then A-Roe. In this test, the A-PVM2 and E-PVM2 are theleast efficient schemes, even less than N-Roe. For all schemes, the differencesin accuracy between the analytical and approximated eigenvalue solvers areinsignificant. Approximated eigenvalues produce almost equal results as analyt-ical ones, but need slightly less computational time. This can be explained by21 h ( m ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 35 40 45 50 55 60 65x (m)0.40.20.00.20.4 u ( m / s ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 Figure 6: Test II: A detail of the interface depth and lower layer velocity obtained by N-Roe,A-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 scheme, compared to the reference solutionat t = 25 s with ∆ x = 0 . satisfactory accuracy of approximated eigenvalues for r values close to one. CPU time (s)10 E h N-RoeA-RoeA-IFCPE_IFCPA-PVM2E-PVM2 CPU time (s)10 E u Figure 7: Test II: CPU time vs error for N-Roe, A-Roe, A-IFCP, E-IFCP, A-PVM2, andE-PVM2 scheme, compared to the reference solution r =0.98 This test considers an internal dam-break over uneven bottom, which resultsin a transition from subcritical to supercritical flow. A small density differencebetween the layers is considered, namely r = 0 . 98. The spatial domain is set to[0, 10], and the bottom topography is defined by a 0.5 m high sill located in thecenter of the channel, given by the function: b ( x ) = 0 . (cid:0) − ( x − (cid:1) . (47)22he initial condition is given by: h ( x, 0) = . , if x < . , otherwise h ( x, 0) = 1 . m − b ( x ) − h ( x, 0) (48) u ( x, 0) = u ( x, 0) = 0 m s − (49)Non-reflective conditions are imposed at the boundaries. Several mesh sizesare considered, namely ∆ x = 1/10, 1/20, 1/40, and 1/80 m. A variable timestep ∆ t is evaluated at each step to satisfy CF L = 0 . 9. The reference solutionis computed using the A-Roe scheme and a dense grid with ∆ x = 1/160 m.Figure 8 shows the temporal evolution of the interface and free-surface pro-files for the reference solution. A detail of solutions for the interface depth andlower layer velocity are shown in Fig. 9 where N-Roe, A-Roe, A-IFCP, E-IFCP,A-PVM2, and E-PVM2 numerical schemes are compared against the referencesolution at t = 15 s with ∆ x = 1 / 20 m. The results indicate that Roe andIFCP schemes are equally accurate, whereas PVM-2U schemes produce notice-ably more diffused results. Same as in the previous case, the implementation ofthe eigenvalue solver has no apparent influence on the accuracy of the results,since approximated eigenvalues exhibit negligible errors for r = 0 . E Φ is presented in Fig. 10. Theresults are very similar to the previous test case, shown in Fig. 7. The E-IFCP isthe most efficient scheme, closely followed by A-IFCP and then A-Roe. The A-PVM2 and E-PVM2 are the least efficient schemes. The differences in efficiencybetween the analytical and approximated eigenvalue solvers are insignificant forIFCP and PVM-2U schemes. r =0.4 A two-layer flow through a rectangular channel with a flat bottom topog-raphy is considered again. In contrast to a similar internal dam-break scenariopresented in [25], here we consider a much larger density difference between thelayers, namely r = 0 . 4. The spatial domain is set to [0, 50], and the initial23 .00.51.01.5 t=0.0 s t=3.0 s t=6.0 s0 2 4 6 8 100.00.51.01.5 t=9.0 s 0 2 4 6 8 10t=12.0 s 0 2 4 6 8 10t=15.0 slength (m) d e p t h ( m ) Figure 8: Test III: Temporal evolution of the interface and surface profile (the referencesolution) h ( m ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 u ( m / s ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 Figure 9: Test III: A detail of the interface depth and lower layer velocity obtained by N-Roe,A-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 scheme, compared to the reference solutionat t = 15 s with ∆ x = 0 . 05 m CPU time (s)10 E h N-RoeA-RoeA-IFCPE_IFCPA-PVM2E-PVM2 CPU time (s)10 E u Figure 10: Test III: CPU time vs error for N-Roe, A-Roe, A-IFCP, E-IFCP, A-PVM2, andE-PVM2 scheme, compared to the reference solution condition is given by: h ( x, 0) = . , if x < 25 m0 . , otherwise h ( x, 0) = 1 . m − h ( x, 0) (50) u ( x, 0) = u ( x, 0) = 0 m s − (51)Non-reflective conditions are imposed at the boundaries. Several mesh sizesare considered, namely ∆ x = 1/2, 1/4, 1/8, and 1/16 m. A variable time step∆ t is evaluated at each step to satisfy CF L = 0 . 9. The reference solution iscomputed using the A-Roe scheme and a dense grid with ∆ x = 1/32 m.Figure 11 shows the temporal evolution of the interface and free-surface pro-files for the reference solution. A detail of solutions for the interface depth andlower layer velocity are shown in Fig. 12 where N-Roe, A-Roe, A-IFCP, E-IFCP,A-PVM2, and E-PVM2 numerical schemes are compared against the referencesolution at t = 5 s with ∆ x = 1 / r ≈ 1, which is not the case here.A CPU time vs. relative root square error E Φ is presented in Fig. 13. Ingeneral, the results show that the A-IFCP, E-IFCP and A-Roe are the most effi-25 .000.250.500.751.00 t=0.0 s t=1.0 s t=2.0 s0 10 20 30 40 500.000.250.500.751.00 t=3.0 s 0 10 20 30 40 50t=4.0 s 0 10 20 30 40 50t=5.0 slength (m) d e p t h ( m ) Figure 11: Test IV: Temporal evolution of the interface and surface profile (the referencesolution) 28 30 32 34 36x (m)0.200.250.300.350.400.450.50 h ( m ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 28 30 32 34 36x (m)0.20.40.60.81.0 u ( m / s ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 Figure 12: Test IV: A detail of the interface depth and lower layer velocity obtained by N-Roe,A-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 scheme, compared to the reference solutionat t = 5 s with ∆ x = 0 . 25 m CPU time (s)10 E h N-RoeA-RoeA-IFCPE_IFCPA-PVM2E-PVM2 CPU time (s)10 E u Figure 13: Test IV: CPU time vs error for N-Roe, A-Roe, A-IFCP, E-IFCP, A-PVM2, andE-PVM2 scheme, compared to the reference solution r =0.4 This test presents a lower-layer column collapse on a sloped bottom. Againa larger density difference between the layers is considered, namely r = 0 . 4. Thespatial domain is set to [0, 40], and the initial condition is given by: h ( x, 0) = . , if x < 25 m0 . , otherwise h ( x, 0) = 1 . m − h ( x, 0) (52) u ( x, 0) = u ( x, 0) = 0 m s − (53)Non-reflective conditions are imposed at the boundaries. Several mesh sizesare considered, namely ∆ x = 2/5, 1/5, 1/10, and 1/20 m. A variable time step∆ t is evaluated at each step to satisfy CF L = 0 . 9. The reference solution iscomputed using the A-Roe scheme and a dense grid with ∆ x = 1/40 m.Figure 14 shows the temporal evolution of the interface and free-surfaceprofiles for the reference solution. A detail of solutions for the interface depth27nd lower layer velocity are shown in Fig. 15 where N-Roe, A-Roe, A-IFCP,E-IFCP, A-PVM2, and E-PVM2 numerical schemes are compared against thereference solution at t = 5 s with ∆ x = 1 / d e p t h ( m ) Figure 14: Test V: Temporal evolution of the interface and surface profile (the referencesolution) A CPU time vs. relative root square error E Φ is presented in Fig. 13. Similarto the previous test, E-IFCP and A-IFCP are the most efficient schemes, closelyfollowed by A-Roe scheme. Both PVM-2U schemes and N-Roe are noticeablyless efficient. The differences in efficiency between the analytical and approx-imated eigenvalue solvers are not significant. For IFCP scheme they seems tobe almost identical, but for PVM-2U scheme, the analytical implementation ismore efficient. 28 h ( m ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 24 26 28 30 32 34x (m)0.20.40.60.81.01.2 u ( m / s ) Ref. sol.N-RoeA-RoeA-IFCP E_IFCPA-PVM2E-PVM2 Figure 15: Test V: A detail of the interface depth and lower layer velocity obtained by N-Roe,A-Roe, A-IFCP, E-IFCP, A-PVM2, and E-PVM2 scheme, compared to the reference solutionat t = 5 s with ∆ x = 0 . CPU time (s)10 E h N-RoeA-RoeA-IFCPE_IFCPA-PVM2E-PVM2 CPU time (s)10 E u Figure 16: Test V: CPU time vs error for N-Roe, A-Roe, A-IFCP, E-IFCP, A-PVM2, andE-PVM2 scheme, compared to the reference solution . Discussion and Conclusion This study re-evaluated the efficiency of Roe, IFCP, and PVM-2U schemesfor two-layer shallow water systems with different solutions for eigenvalues. Forthis purpose, numerical, approximated, and recently-proposed analytical solu-tions for the eigenstructure were considered, combined with Roe, PVM-2U andIFCP numerical schemes. The choice of eigenvalues in numerical schemes fortwo-layer SWE models were evaluated in three stages.First analysis focused only on the accuracy and computational time of dif-ferent eigenvalue solvers. Numerical and analytical eigenvalue solvers producealmost identical results. Approximated expressions, on the other hand, areless accurate, and the errors grow with density difference between the layers.This increase in errors appears because approximate solutions deviate from theinitial assumption that the density ratio is close to one. Regarding the computa-tional time, analytical and approximated expressions are one order of magnitudefaster than the numerical solver, with approximate expressions being two timesfaster than analytical ones. This additionally confirms our previous study [25]where analytical expressions were found to be equally accurate but one orderof magnitude faster than numerical solvers (approximated expressions were notconsidered in the previous study).The second analysis shifted the focus from eigensolvers to the numerical vis-cosity matrix, which was computed using different numerical schemes in com-bination with different eigenvalue solvers. The results revealed that the Roemethod is highly sensitive to the choice of eigenvalues, IFCP method is some-what sensitive, and PVM-2U method shows very little sensitivity to the choiceof eigenvalue. In general, Roe method is the most accurate, followed by IFCP,and then PVM-2U method. This is expected since IFCP approximates the vis-cosity matrix using four eigenvalues, and PVM-2U using only two eigenvalues.Also, all schemes are more accurate when analytical eigenvalues are used incomparison to approximated values. It is interesting that the Roe method withapproximated eigenvalues is the least accurate method overall for large density30ifferences. The computational time needed to obtain the numerical viscositymatrix is the longest for Roe scheme with numerical eigensolver. Other com-binations of numerical schemes and eigensolvers are several times faster, andvery close to each other. The rankings from the fastest to the slowest combina-tion are: PVM-2U scheme with approximated eigenvalues, IFCP with approxi-mated eigenvalues, PVM-2U with analytical eigenvalues, IFCP with analyticaleigenvalues, Roe with approximated eigenvalues, and finally Roe with analyticaleigenvalues. It is important to note that only 10-20% of total time needed toconstruct the viscosity matrix is spent on computing the eigenvalues.In the third analysis several numerical tests were performed to investigate theoverall performance of different numerical schemes in combination with differenteigenvalue solvers. A total of five tests were designed to evaluate more realisticscenarios, including different density ratios between the layers and various chan-nel geometries. The results revealed that the Roe scheme with approximatedeigenvalues is not well-balanced, and should not be considered in two-layer mod-elling. Remaining four test showed that IFCP (with approximated and analyt-ical eigenvalues) and Roe scheme (with analytical eigenvalues) are very closein performance, with IFCP being slightly better. PVM-2U was noticeably lessefficient, regardless of the choice for the eigensolver.These findings instil more confidence into findings from out previous study[25] that showed how Roe method with analytical eigenvalues is very close toIFCP scheme with approximated eigenvalues. Note that the previous studydid not considered the modifications to the numerical schemes presented inSection 2.4, which may improve their computational speed. Also note that theconclusions from the previous study were based only on two numerical exampleswith the same (small) density difference between the layers and used a fixedtime step. Furthermore, the previous study did not evaluate the impact ofimplementing analytical eigenvalues into PVM and IFCP schemes. Also, thePVM-2U is considered in the present study, whereas [25] evaluated PVM-Roewhich is a redefinition of the Roe scheme under the PVM paradigm. PVM-2U,on the other hand, is a new PVM method based on two external eigenvalues.31verall, when modelling layers of smaller density difference, it seems thatapproximated eigenvalues are a more efficient choice for IFCP and PVM-2Uschemes. For larger density differences, however, the analytical eigenvalues areas efficient. It should be emphasized that analytical eigenvalues are more preciseand can help with an accurate prediction of hyperbolicity losses for all densityratios, without producing any overhead in computational time.The extension of all three considered schemes (Roe, IFCP, PVM-2U) to ahigher order is straightforward, following a general approach presented in [13].An extension to a two-dimensional case is also possible following the procedurefrom [11]. A high order extension of the family of PVM methods is similar tothe Roe method [10]. For example, Castro et al. [9] presented the extensionof IFCP method to solve a two-layer Savage-Hutter type model and simulatetsunamis generated by landslides. Since the eigenvalues affect only the viscositymatrix, the extension of each considered scheme to a higher order is the sameregardless of the eigenvalue solver. On the other hand, some impact of theeigenvalues on the efficiency of higher order methods is expected and should befurther analysed.Although the efficiency of analytical solutions to the eigenstructure havebeen assessed here for two-layer shallow-water flows, these closed-form solutionsto eigenvalues expressed in coefficient of a characteristic quartic, can directlybe applied to some other non-conservative hyperbolic systems defined by fourcoupled partial differential equations, such as two-phase granular flows. Acknowledgements This work has been fully supported by the University of Rijeka under theproject number 17.06.2.1.02 (River-Sea Interaction in the Context of ClimateChange). 32 eferences [1] Abgrall, R., Karni, S., 2009. Two-layer shallow water system: a relaxation ap-proach. SIAM Journal on Scientific Computing 31, 1603–1627. URL: https://doi.org/10.1137/06067167X .[2] Abramowitz, M., Stegun, I.A., 1965. Handbook of mathematical functions withformulas, graphs, and mathematical tables. Dover Publications, New York.[3] Adduce, C., Sciortino, G., Proietti, S., 2011. Gravity currents produced by lockexchanges: Experiments and simulations with a two-layer shallow-water modelwith entrainment. Journal of Hydraulic Engineering 138, 111–121. URL: https://doi.org/10.1061/(ASCE)HY.1943-7900.0000484 .[4] Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J.,Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, D.,1999. LAPACK Users’ Guide. Third ed., Society for Industrial and AppliedMathematics, Philadelphia, PA.[5] Bermudez, A., Vazquez, M.E., 1994. Upwind methods for hyperbolic conservationlaws with source terms. Computers & Fluids 23, 1049–1071. URL: https://doi.org/10.1016/0045-7930(94)90004-3 .[6] Bouchut, F., Zeitlin, V., 2010. A robust well-balanced scheme for multi-layershallow water equations. Discrete and Continuous Dynamical Systems-Series B13, 739–758. URL: http://dx.doi.org/10.3934/dcdsb.2010.13.739 .[7] Canestrelli, A., Fagherazzi, S., Lanzoni, S., 2012. A mass-conservative centeredfinite volume model for solving two-dimensional two-layer shallow water equationsfor fluid mud propagation over varying topography and dry areas. Advancesin Water Resources 40, 54–70. URL: https://doi.org/10.1016/j.advwatres.2012.01.009 .[8] Canestrelli, A., Toro, E.F., 2012. Restoration of the contact surface in FORCE-type centred schemes II: Non-conservative one-and two-layer two-dimensionalshallow water equations. Advances in Water Resources 47, 76–87. URL: https://doi.org/10.1016/j.advwatres.2012.03.018 . 9] Castro, M., de la Asunci´on, M., Macias, J., Par´es, C., Fern´andez-Nieto, E.,Gonz´alez-Vida, J., de Luna, T.M., 2012. Ifcp riemann solver: Application totsunami modelling using gpus. Numerical Methods for Hyperbolic Equations:Theory and Applications , 237–244.[10] Castro, M.J., Fern´andez-Nieto, E.D., 2012. A class of computationally fast firstorder finite volume solvers: PVM methods. SIAM Journal on Scientific Comput-ing 34, A2173–A2196. URL: https://dx.doi.org/10.1137/100795280 .[11] Castro, M.J., Fern´andez-Nieto, E.D., Ferreiro, A., Garc´ıa-Rodr´ıguez, J.A., Par´es,C., 2009. High order extensions of roe schemes for two-dimensional nonconserva-tive hyperbolic systems. Journal of Scientific Computing 39, 67–114.[12] Castro, M.J., Fern´andez-Nieto, E.D., Gonz´alez-Vida, J.M., Par´es-Madronal, C.,2011. Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system. Journal of Scientific Computing 48, 16–40. URL: https://doi.org/10.1007/s10915-010-9427-5 .[13] Castro, M.J., Gallardo, J., Par´es, C., 2006. High order finite volume schemesbased on reconstruction of states for solving hyperbolic systems with nonconser-vative products. applications to shallow-water systems. Mathematics of compu-tation 75, 1103–1134.[14] Castro, M.J., Garcıa-Rodrıguez, J.A., Gonz´alez-Vida, J.M., Macıas, J., Par´es, C.,V´azquez-Cend´on, M.E., 2004. Numerical simulation of two-layer shallow waterflows through channels with irregular geometry. Journal of Computational Physics195, 202–235. URL: https://doi.org/10.1016/j.jcp.2003.08.035 .[15] Castro, M.J., Mac´ıas, J., Par´es, C., 2001. A Q-scheme for a class of systemsof coupled conservation laws with source term. Application to a two-layer 1-Dshallow water system. ESAIM: Mathematical Modelling and Numerical Analysis35, 107–127. URL: https://doi.org/10.1051/m2an:2001108 .[16] Castro, M.J., Pardo, A., Par´es, C., Toro, E., 2010. On some fast well-balancedfirst order solvers for nonconservative systems. Mathematics of Computation 79,1427–1472. URL: http://doi.org/10.1090/S0025-5718-09-02317-5 . 17] Chakir, M., Ouazar, D., Taik, A., 2009. Roe scheme for two-layer shallow waterequations: Application to the Strait of Gibraltar. Mathematical Modelling of Nat-ural Phenomena 4, 114–127. URL: https://doi.org/10.1051/mmnp/20094508 .[18] Chiapolino, A., Saurel, R., 2018. Models and methods for two-layer shallowwater flows. Journal of Computational Physics 371, 1043–1066. URL: https://doi.org/10.1016/j.jcp.2018.05.034 .[19] Doyle, E.E., Hogg, A.J., Mader, H.M., 2011. A two-layer approach to mod-elling the transformation of dilute pyroclastic currents into dense pyroclasticflows. Proceedings of the Royal Society of London A: Mathematical, Physicaland Engineering Sciences 467, 1348–1371. URL: http://dx.doi.org/10.1098/rspa.2010.0402 .[20] Fern´andez-Nieto, E.D., Bouchut, F., Bresch, D., Castro, M.J., Mangeney, A.,2008. A new Savage–Hutter type model for submarine avalanches and generatedtsunami. Journal of Computational Physics 227, 7720–7754. URL: https://doi.org/10.1016/j.jcp.2008.04.039 .[21] Fern´andez-Nieto, E.D., Castro, M.J., Par´es, C., 2011. On an intermediate fieldcapturing Riemann solver based on a parabolic viscosity matrix for the two-layer shallow water system. Journal of Scientific Computing 48, 117–140. URL: https://doi.org/10.1007/s10915-011-9465-7 .[22] Krvavica, N., 2019. Q-2L-SWE, Viscosity matrix for two-layer shallow waterequations. URL: https://github.com/NKrvavica/Q-2L-SWE . [Online; accessed25-February-2019].[23] Krvavica, N., Koˇzar, I., Travaˇs, V., Oˇzani´c, N., 2017a. Numerical modelling oftwo-layer shallow water flow in microtidal salt-wedge estuaries: Finite volumesolver and field validation. Journal of Hydrology and Hydromechanics 65, 49–59.URL: https://doi.org/10.1515/johh-2016-0039 .[24] Krvavica, N., Travaˇs, V., Oˇzani´c, N., 2017b. Salt-Wedge Response to VariableRiver Flow and Sea-Level Rise in the Microtidal Rjeˇcina River Estuary, Croa-tia. Journal of Coastal Research 33, 802–814. URL: https://doi.org/10.2112/JCOASTRES-D-16-00053.1 . 25] Krvavica, N., Tuhtan, M., Jeleni´c, G., 2018. Analytical implementation of roesolver for two-layer shallow water equations with accurate treatment for loss ofhyperbolicity. Advances in Water Resources 122, 187–205. URL: https://doi.org/10.1016/j.advwatres.2018.10.017 .[26] Kurganov, A., Petrova, G., 2009. Central-upwind schemes for two-layer shallowwater equations. SIAM Journal on Scientific Computing 31, 1742–1773. URL: https://doi.org/10.1137/080719091 .[27] La Rocca, M., Adduce, C., Sciortino, G., Pinzon, A.B., Boniforti, M.A., 2012.A two-layer, shallow-water model for 3D gravity currents. Journal of HydraulicResearch 50, 208–217. URL: https://doi.org/10.1080/00221686.2012.667680 .[28] Lee, W.K., Borthwick, A.G., Taylor, P.H., 2011. A fast adaptive quadtree schemefor a two-layer shallow water model. Journal of Computational Physics 230, 4848–4870. URL: https://doi.org/10.1016/j.jcp.2011.03.007 .[29] Luca, I., Hutter, K., Kuo, C., Tai, Y., 2009. Two-layer models for shallowavalanche flows over arbitrary variable topography. International Journal ofAdvances in Engineering Sciences and Applied Mathematics 1, 99–121. URL: https://doi.org/10.1007/s12572-010-0006-7 .[30] Majd, M.S., Sanders, B.F., 2014. The LHLLC scheme for two-layer and two-phase transcritical flows over a mobile bed with avalanching, wetting and dry-ing. Advances in water resources 67, 16–31. URL: https://doi.org/10.1016/j.advwatres.2014.02.002 .[31] Murillo, J., Garc´ıa-Navarro, P., 2010. An Exner-based coupled model for two-dimensional transient flow over erodible bed. Journal of Computational Physics229, 8704–8732. URL: https://doi.org/10.1016/j.jcp.2010.08.006 .[32] Ovsyannikov, L., 1979. Two-layer shallow water model. Journal of Applied Me-chanics and Technical Physics 20, 127–135.[33] Par´es, C., 2006. Numerical methods for nonconservative hyperbolic systems: atheoretical framework. SIAM Journal on Numerical Analysis 44, 300–321. URL: https://doi.org/10.1137/050628052 . 34] Par´es, C., Castro, M.J., 2004. On the well-balance property of roe’s methodfor nonconservative hyperbolic systems. applications to shallow-water systems.ESAIM: Mathematical Modelling and Numerical Analysis 38, 821–852. URL: https://doi.org/10.1051/m2an:2004041 .[35] Pelanti, M., Bouchut, F., Mangeney, A., 2008. A Roe-type scheme for two-phaseshallow granular flows over variable topography. ESAIM: Mathematical Modellingand Numerical Analysis 42, 851–885. URL: https://doi.org/10.1051/m2an:2008029 .[36] Pitman, E.B., Le, L., 2005. A two-fluid model for avalanche and debris flows.Philosophical Transactions of the Royal Society of London A: Mathematical,Physical and Engineering Sciences 363, 1573–1601. URL: https://doi.org/10.1098/rsta.2005.1596 .[37] Sarno, L., Carravetta, A., Martino, R., Papa, M., Tai, Y.C., 2017. Some consider-ations on numerical schemes for treating hyperbolicity issues in two-layer models.Advances in Water Resources 100, 183–198. URL: https://doi.org/10.1016/j.advwatres.2016.12.014 .[38] Schijf, J., Sch¨onfled, J., 1953. Theoretical considerations on the motion of saltand fresh water, in: Proceedings Minnesota International Hydraulic Convention,IAHR.[39] Toumi, I., 1992. A weak formulation of roe’s approximate riemann solver. Journalof computational physics 102, 360–373.[40] Wikipedia, 2018. Quartic function — wikipedia, the free encyclope-dia. URL: https://upload.wikimedia.org/wikipedia/commons/9/99/Quartic_Formula.svg . [Online; accessed 19-February-2019].. [Online; accessed 19-February-2019].