Discussion on different numerical treatments on the loss of hyperbolicity for the two-layer shallow water system

Advances in Water Resources 182 (2023) 104587 Available online 23 November 2023 0309-1708/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/bync/4.0/). Contents lists available at ScienceDirect Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres Discussion on different numerical treatments on the loss of hyperbolicity for the two-layer shallow water system M.J. Castro Díaza, E.D. Fernández-Nietob, J. Garres-Díazc,∗, T. Morales de Lunaa aDpto. Análisis Matemático, Estadística e Investigación Operativa, y Matemática Aplicada, Universidad de Málaga, 29071 Málaga, Spain bDpto. Matemática Aplicada I, Universidad de Sevilla, 41012 Sevilla, Spain cDpto. Matemática Aplicada II, Universidad de Sevilla, 41092 Sevilla, Spain ARTICLE INFO Keywords: Two-layer shallow water system Hyperbolicity region Finite volume PVM methods Submarine landslides ABSTRACT This paper focus on the numerical approximation of two-layer shallow water system. First, a new approximation of the eigenvalues of the system is proposed, which satisfies some interesting properties. From this approximation, we give an accurate estimation of the non-hyperbolic region, which improves significantly the one computed with the classic eigenvalues approximation. In particular, we estimate both the lower and upper boundaries of the non-hyperbolic region. We also give a simple algorithm that allows us to compute bounds for the external eigenvalues, even when complex eigenvalues arise. We design an efficient FV solver depending on the hyperbolic nature of the system, that combines the PVM-2U-FL scheme and a new solver introduced in this work, named IFCP-FL, which degenerates to the Lax–Wendroff method in smooth areas. Different strategies for the numerical treatment for the loss of hyperbolicity are considered and discussed. Two of them are based on adding a friction term, depending on the classic or the new eigenvalues approximations, and the third one considers a simplified hyperbolic model in areas with complex eigenvalues. Some numerical tests are performed, including the case of two-layer fluids with different ratio of densities. The application to a two-layer model for submarine landslides is also considered. In the latter case, we show how the treatments based on adding friction are not appropriate for this kind of applications, whereas the treatment based on changing to a hyperbolic model produces much better results. 1. Introduction In this work we study the ideal case where a flow is composed of two immiscible layers of different densities in a shallow domain, that is, where the typical horizontal length is much larger than the vertical one. Such situations arise naturally in estuary flows, ocean currents, and atmospheric flows. This kind of flows are of interest from both a theoretical and practical point of view and they have frequently been used to study exchange flows through channels connecting two basins with different hydrological characteristics (see Armi (1986), Armi and Farmer (1986), Helfrich (1995) and Farmer and Armi (1986) among many others). In particular, one practical situation of this type of fluids can be found at the Strait of Gibraltar (see Izquierdo et al. (2001) and Tejedor et al. (1998)), where a two-layer exchange flows occurs: the less saline Atlantic water flowing at surface and penetrating into the Mediterranean, and the deeper, denser Mediterranean water flowing into the Atlantic. Another example occurs in the case of submarine avalanches, or aerial avalanches that reach a fluid (see Fernández-Nieto et al. (2008)). ∗Corresponding author. E-mail addresses: [email protected] (M.J. Castro Díaz), [email protected] (E.D. Fernández-Nieto), [email protected] (J. Garres-Díaz), [email protected] (T. Morales de Luna). The formulation corresponding to the two-layer system falls within the framework of hyperbolic conservation laws with non-conservative products and source terms. From the theoretical point of view, the model presents several difficulties. The first one is related precisely to the presence of non-conservative products. Classically, such products are not well-defined for discontinuous solutions, although they could be defined, for instance, by means of the theory developed in Dal Maso et al. (1995) (see also Hou and LeFloch (1994) and Crasta and LeFloch (2002)). The second one is related to the character of the system and its eigenvalues in particular. Indeed, although the model is hyperbolic in general, the eigenvalues may eventually become complex. Such situation is related to Kelvin–Helmholtz instability and the physical behaviour corresponds to mixing and turbulence (see for instance von Helmholtz,1868;Thomson,1871;Harris et al.,2020). Although one could say that in such situations the two-layer system is no longer valid, any numerical scheme should treat them so that computations do not https://doi.org/10.1016/j.advwatres.2023.104587 Received 3 July 2023; Received in revised form 13 November 2023; Accepted 16 November 2023 Advances in Water Resources 182 (2023) 104587 2 M.J. Castro Díaz et al. Fig. 1. Two-layer sketch: unknowns ℎ1,ℎ2and bottom 𝑏. The internal interface is 𝜂2 and the free surface 𝜂1. fail in the case of an eventual loss of hyperbolicity for a small region or time. Two-layer shallow water systems are defined by the depth-averaged mass and momentum conservation equations for two immiscible fluids, whose constant densities are 𝜌1, 𝜌2respectively, over a variable topography 𝑏(𝑥). We shall consider upper (index 𝑖= 1) and lower (index 𝑖= 2) layers with heights ℎ𝑖and discharges 𝑞𝑖, for 𝑖= 1,2(𝑢𝑖=𝑞𝑖∕ℎ𝑖 will denote the layer velocities), see Fig. 1. Let us remark that one of the layers could disappear in some parts of the domain in practical applications, as the case of aerial avalanches reaching a fluid area. The two-layer system reads (see e.g. Castro-Díaz et al. (2010)) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜕𝑡ℎ1+𝜕𝑥𝑞1= 0, 𝜕𝑡𝑞1+𝜕𝑥(𝑞2 1∕ℎ1)+𝑔ℎ1𝜕𝑥(ℎ1+ℎ2+𝑏)= 0, 𝜕𝑡ℎ2+𝜕𝑥𝑞2= 0, 𝜕𝑡𝑞2+𝜕𝑥(𝑞2 2∕ℎ2)+𝑔ℎ2𝜕𝑥(𝑟ℎ1+ℎ2+𝑏)= 0, (1) being 𝑔the gravitational acceleration and 𝑟=𝜌1∕𝜌2∈(0,1)the density ratio. For the sake of simplicity we have not considered friction terms between the two layers or the bottom and we shall focus on the 1D case, although the eigenvalues approximation introduced in this work can be also used in the 2D case, as it is explained in Section 3.3. In the case of submarine avalanches the right-hand side of the momentum equation would include a Coulomb friction term. The details and adaptation of the numerical scheme for submarine avalanches model is detailed in Section 3.4. For convenience, we shall denote by 𝜂1, 𝜂2the free surface and interface between the layers, described by 𝜂1=𝑏+ℎ1+ℎ2and 𝜂2=𝑏+ℎ2. Notice that the lake-at-rest steady solutions of system (1) are 𝑞1= 0, ℎ1=constant, 𝑞2= 0, ℎ2+𝑏=constant. These conditions are equivalent to having a flow at rest (𝑞𝑖= 0) and 𝜂1, 𝜂2constant. As said previously, system (1) includes non-conservative products, which correspond to the terms ℎ1𝜕𝑥ℎ2and ℎ2𝜕𝑥ℎ1, coming from the coupling between the layers. In this work we shall use the theory of Dal Maso et al. (1995) combined with the framework of pathconservative methods introduced in Parés (2006). The system may be written in a quasi-linear compact form as 𝜕𝑡𝒘+𝐴(𝒘)𝜕𝑥𝒘+𝑆(𝒘)𝜕𝑥𝑏= 0 (2) where 𝒘=(ℎ1, 𝑞1, ℎ2, 𝑞2)′is the conservative variable vector, and 𝐴(𝒘)∈4×4 and 𝑆(𝒘)∈R4are 𝐴(𝒘)=⎛⎜⎜⎜⎜⎝ 0100 −𝑢2 1+𝑔ℎ12𝑢1𝑔ℎ10 0001 𝑔𝑟ℎ20 −𝑢2 2+𝑔ℎ22𝑢2 ⎞⎟⎟⎟⎟⎠ , 𝑆 (𝒘)=⎛⎜⎜⎜⎜⎝ 0 𝑔ℎ1 0 𝑔ℎ2 ⎞⎟⎟⎟⎟⎠ . Concerning the eigenstructure of system (1), we will study the eigenvalues of the matrix 𝐴(𝒘), whose characteristic polynomial is 𝑃(𝜆) = (𝜆2− 2𝑢1𝜆+𝑢2 1−𝑔ℎ1)(𝜆2− 2𝑢2𝜆+𝑢2 2−𝑔ℎ2)−𝑟𝑔2ℎ1ℎ2.(3) A simple computation reveals that 𝜆= 0 is an eigenvalue in the particular case that the so-called composite Froude number (𝐺) satisfies 𝐺2=𝐹𝑟2 1+𝐹𝑟2 2− (1 − 𝑟)𝐹𝑟2 1𝐹𝑟2 2= 1,(4) where 𝐹𝑟𝑖,𝑖= 1,2are the internal Froude numbers given by 𝐹2 𝑖= 𝑢2 𝑖∕(𝑔(1 − 𝑟)ℎ𝑖). Otherwise, the eigenvalues will be non-zero when (4) is not verified. When solving numerically the two-layer model, it is important to have some estimation on the eigenvalues of the system. Although the eigenvalues correspond to the roots of a fourth order polynomial, which may be computed using the Cardano–Vieta formulas (Weisstein,2002), there is no easy explicit formulation for them. Therefore, from the numerical point of view, it is desirable to have some accurate, easy and explicit formulation of the eigenvalues. The eigenvalues of system (1) may be split into two internal (𝜆± int) and two external (𝜆± ext). Approximations of these eigenvalues were introduced in the classical work of Schijf and Schöfled (1953). They are a first order approximation of the true eigenvalues in the case 𝑟≈ 1 and 𝑢1≈𝑢2. Moreover, they allow to estimate when hyperbolicity of the system is lost. More explicitly, denoting by  𝜆𝑆𝑆 ± ext , 𝜆𝑆𝑆 ± int the approximation in Schijf and Schöfled (1953), they are given by  𝜆𝑆𝑆 ± ext =𝑢1ℎ1+𝑢2ℎ2 ℎ1+ℎ2 ±√𝑔(ℎ1+ℎ2),  𝜆𝑆𝑆 ± int =𝑢1ℎ2+𝑢2ℎ1 ℎ1+ℎ2 ±√𝑔(1 − 𝑟)ℎ1ℎ2 (ℎ1+ℎ2)(1 − (𝑢1−𝑢2)2 𝑔(1 − 𝑟)(ℎ1+ℎ2)), (5) where it holds that ⎛⎜⎜⎝  𝜆𝑆𝑆 ± ext  𝜆𝑆𝑆 ± int ⎞⎟⎟⎠ =⎛⎜⎜⎝ 𝜆± ext 𝜆± int⎞⎟⎟⎠ +(𝑢1−𝑢2). Remark that the approximated eigenvalues  𝜆𝑆𝑆 ± ext are always real values, but that is not the case for  𝜆𝑆𝑆 ± int , which become complex whenever √𝑔(1 − 𝑟)(ℎ1+ℎ2)<|𝑢1−𝑢2|.(6) Therefore, the hyperbolicity of the system is expected to be lost approximately when (6) is satisfied. Let us remark that the approximation of the eigenvalues (5) may be also seen as a particular case of the asymptotic analysis introduced in Sarno et al. (2021) in the case 𝑟≈ 1 and 𝑢1≈𝑢2. Notice that in Sarno et al. (2021), the asymptotic analysis of the eigenvalues and hyperbolicity region is analysed for several cases, namely for very close and very different layer heights. As said previously and explained in Castro-Díaz et al. (2010), this loss of hyperbolicity is related to the appearance of shear instabilities that may lead to intense mixing of the two layers (see Harris et al. (2020)). In Castro-Díaz et al. (2010) a numerical treatment to deal with this difficulty is proposed, which is based on adding an artificial friction between the layers. This friction would represent the energy that would be lost in a mixing situation, making ||𝑢1−𝑢2||smaller, so that we return back to the hyperbolic region of the system. Several other numerical treatments of the loss of hyperbolicity have been proposed in the bibliography, either using the estimations of the eigenvalues given in Schijf and Schöfled (1953), by approximating them numerically or using the Cardano-Vieta formula. Among many others we may cite (Castro et al.,2001,2004;Bouchut and Morales de Luna, 2008;Kurganov and Petrova,2009;Abgrall and Karni,2009;Castro- Díaz et al.,2010;Berthon et al.,2015;Sarno et al.,2017;Krvavica et al.,2018;Murillo et al.,2020). One of the goals of this work is to give a good approximation of the hyperbolicity region, as well as Advances in Water Resources 182 (2023) 104587 3 M.J. Castro Díaz et al. a novel numerical treatment when complex eigenvalues are found. It is expected that system (1) loses the hyperbolicity in many practical scenarios. Actually, our main motivation is the application of the two layer system to submarine avalanche simulations, where 𝑟≈ 0.33 and such scenario is easily found in the first times of the formation of the avalanche. The paper is organized as follows: In Section 2a new approximation of the eigenvalues of the system is proposed. It allows us to improve the estimation of the non-hyperbolic region. We also propose an algorithm that allows us to get bounds of the external eigenvalues. Section 3deals with the numerical discretization of the two-layer shallow water system by a finite volume method. We will propose a combination between two finite volume solvers, the PVM-2U-FL and IFCP-FL, being the last one an original contribution of this work. This new solver is introduced in Section 3.1.1. Three numerical treatments for the case of complex eigenvalues are presented in Section 3.2, two of them are based on adding a friction correction term, whereas the last one is based on changing locally the model. The extension to 2D horizontal domains is briefly discussed in Section 3.3. In Section 3.4 we explain how the two-layer system and the proposed new numerical scheme is adapted for the case of a submarine avalanches model. Then, Section 4will present several numerical tests. They are divided into four blocks of test cases: first a steady two-layer test to evaluate the new IFCP-FL method (Section 4.1), then one for high density ratio (Section 4.2) and one for low one (Section 4.3), and finally the application to submarine landslides (Section 4.4). In particular, in Sections 4.1 and 4.2.1 we show the better accuracy of the proposed IFCP-FL solver in comparison with other schemes. The difference between the three numerical treatments applied in the presence of complex eigenvalues are then analysed in the other test cases. Conclusions are presented in Section 5. 2. A new approximation for the eigenvalues and hyperbolicity region for the two-layer system This section is devoted to introducing a new approximation of the eigenvalues of system (1), whose properties are studied. In particular, we give an algorithm ensuring the approximation of bounds of the external eigenvalues of the system. Based on this new approximations, we are able to improve the estimation of the hyperbolicity region given by the classical approximation (6). Concretely, after an heuristic study, we propose to define the hyperbolicity region as a combination of two different curves. 2.1. Eigenvalues approximation Let us start with the approximation of the eigenvalues of system (1). Looking at the characteristic polynomial (3), we see that it is possible to easily compute the exact eigenvalues in some particular cases. Namely, if 𝑢1=𝑢2, then the exact eigenvalues are 𝜆± ext =𝑢1±√ √ √ √𝑔ℎ1+ℎ2 2+√(𝑔ℎ1+ℎ2 2)2 −𝑔2(1 − 𝑟)ℎ1ℎ2, 𝜆± int =𝑢1±√ √ √ √𝑔ℎ1+ℎ2 2−√(𝑔ℎ1+ℎ2 2)2 −𝑔2(1 − 𝑟)ℎ1ℎ2, (7) and in the case ℎ1=ℎ2, the exact eigenvalues are 𝜆± ext =𝑢1+𝑢2 2±√(𝑢1−𝑢2)2 4+𝑔ℎ1+√𝑟𝑔2ℎ2 1+ℎ1𝑔(𝑢1−𝑢2)2, 𝜆± int =𝑢1+𝑢2 2±√(𝑢1−𝑢2)2 4+𝑔ℎ1−√𝑟𝑔2ℎ2 1+ℎ1𝑔(𝑢1−𝑢2)2. (8) We look for an approximation of the eigenvalues in such a way that it matches with the exact ones in these particular cases (𝑢1=𝑢2or ℎ1=ℎ2). The definition for these new approximations are given in the following theorem: Theorem 1. Let  𝜆± ext, 𝜆± int be the approximations defined by  𝜆± ext =𝑢1ℎ1+𝑢2ℎ2 ℎ1+ℎ2 ±𝑐ext, 𝜆± int =𝑢1ℎ2+𝑢2ℎ1 ℎ1+ℎ2 ±𝑐int,(9a) where, 𝑐ext =√𝑔ℎ1+ℎ2 2+ℎ1ℎ2 (ℎ1+ℎ2)2(𝑢1−𝑢2)2+𝛾𝑐,(9b) 𝑐int =√𝑔ℎ1+ℎ2 2+ℎ1ℎ2 (ℎ1+ℎ2)2(𝑢1−𝑢2)2−𝛾𝑐,(9c) and 𝛾𝑐=√(𝑔ℎ1+ℎ2 2)2 + (𝑟− 1)𝑔2ℎ1ℎ2+2𝑔ℎ1ℎ2 ℎ1+ℎ2 (𝑢1−𝑢2)2.(9d) Then, the following properties are satisfied: 1. If 𝑢1=𝑢2or ℎ1=ℎ2, then  𝜆± ext and  𝜆± int match with the exact values of the external and internal eigenvalues, respectively. 2. Assuming (1 − 𝑟) ∼ (𝛿2)and (𝑢1−𝑢2) ∼ (𝜖), then: (a)  𝜆± ext and  𝜆± int are approximations of order (𝜖+𝛿)and (𝜖2+ 𝛿2+𝜖𝛿)of the classical approximations  𝜆𝑆𝑆 ± ext and  𝜆𝑆𝑆 ± int (5), respectively. (b)  𝜆± ext, and  𝜆± int are approximations of order (𝜖+𝛿)and (𝜖+𝛿2)of the exact external and internal eigenvalues, respectively. 3. Let us denote by 𝑃(𝑎)(𝜆)the characteristic polynomial associated to the proposed approximations  𝜆± ext and  𝜆± int, i.e., 𝑃(𝑎)(𝜆) = (𝜆− 𝜆− ext)(𝜆− 𝜆+ ext)(𝜆− 𝜆− int)(𝜆− 𝜆+ int).(10) Then: (a) If ℎ1=ℎ2or 𝑢1=𝑢2, then 𝑃(𝑎)(𝜆) = 𝑃(𝜆)for all 𝜆∈R, being 𝑃(𝜆)the characteristic polynomial associated to the exact eigenvalues of the system (3). (b) If ℎ1≠ℎ2and 𝑢1≠𝑢2, then there exist a unique value 𝜆0 such that 𝑃(𝑎)(𝜆0) = 𝑃(𝜆0). Moreover, we have that 𝜆0=𝑢1+𝑢2 2−ℎ1ℎ2(ℎ1−ℎ2) (ℎ1+ℎ2)3 (𝑢1−𝑢2)3 (2𝛾𝑐−𝑔(ℎ1+ℎ2)) . Proof. In what follows, the main steps of the proof for each statement are detailed. 1. It is easy to check that one gets that  𝜆± ext, 𝜆± int correspond to (7) or (8) in the case 𝑢1=𝑢2or ℎ1=ℎ2respectively. 2. First, note that the property (2b) is a straight consequence of (2a), since  𝜆𝑆𝑆 ± ext and  𝜆𝑆𝑆 ± int are first order approximations ((𝜖)) of the exact external and internal eigenvalues, respectively, provided that 𝑟≈ 1. Notice that here we have quantified (1 − 𝑟) ∼ (𝛿2). Then, it is enough to prove (2a). Using the Taylor’s theorem, we get √𝑎2+𝑥=𝑎+𝑥 2𝑎+(𝑥2)for any 𝑎 > 0. Then, we write 𝛾𝑐as 𝛾𝑐=𝑔ℎ1+ℎ2 2+(𝑟− 1)𝑔ℎ1ℎ2 ℎ1+ℎ2 +2ℎ1ℎ2(𝑢1−𝑢2)2 (ℎ1+ℎ2)2+(𝜖2+𝛿4+𝛿2𝜖2). Consequently, we get 𝑐ext =√𝑔(ℎ1+ℎ2) + (𝜀+𝛿), and 𝑐int =√(1 − 𝑟)𝑔ℎ1ℎ2 ℎ1+ℎ2(1 − (𝑢1−𝑢2)2 𝑔(1 − 𝑟)(ℎ1+ℎ2))+(𝜖2+𝛿2+𝛿𝜖), which proves (2a). Notice that  𝜆± int are indeed second order approximations of  𝜆𝑆𝑆 ± int , whereas  𝜆± ext are first order approximations of 𝜆± ext. Advances in Water Resources 182 (2023) 104587 4 M.J. Castro Díaz et al. 3. In order to prove that 𝑃(𝜆)intersects 𝑃(𝑎)in a unique point (𝜆0), it is enough to prove that 𝑃(𝜆)−𝑃(𝑎)(𝜆) = 𝑑1𝜆+𝑑0, with 𝑑0, 𝑑1∈R. First, the exact characteristic polynomial (3) is written as 𝑃(𝜆) = 𝜆4+(−2𝑢2− 2𝑢1)𝜆3+((𝑢1+𝑢2)2+ 2𝑢1𝑢2−𝑔(ℎ1+ℎ2))𝜆2 +(−2𝑢1𝑢2(𝑢1+𝑢2) + 2𝑢2𝑔ℎ1+ 2𝑢1𝑔ℎ2)𝜆+𝑢2 1𝑢2 2−𝑢2 1𝑔ℎ2 −𝑢2 2𝑔ℎ1+ (1 − 𝑟)𝑔2ℎ1ℎ2. We also develop the expression of the approximated polynomial 𝑃(𝑎)(𝜆) (10), and looking at the factors multiplying 𝜆3and 𝜆2, it is easy to check that −( 𝜆− int + 𝜆+ int + 𝜆− ext + 𝜆+ ext) = −2𝑢1− 2𝑢2, and ( 𝜆+ int)( 𝜆− int)+( 𝜆+ ext)( 𝜆− ext)+( 𝜆+ ext + 𝜆− ext)( 𝜆+ int + 𝜆− int) = (𝑢1+𝑢2)2+ 2𝑢1𝑢2−𝑔(ℎ1+ℎ2). Therefore, we get that 𝑃(𝜆) − 𝑃(𝑎)(𝜆) = 𝑑1𝜆+𝑑0, where, after some computations, we obtain that 𝑑1= (𝑢1−𝑢2)(ℎ1−ℎ2)(2𝛾𝑐 ℎ1+ℎ2 −𝑔), and 𝑑0=ℎ1ℎ2(ℎ1−ℎ2)2 (ℎ1+ℎ2)4(𝑢1−𝑢2)4−1 2(ℎ1−ℎ2)(𝑢2 1−𝑢2 2)(2𝛾𝑐 ℎ1+ℎ2 −𝑔) =ℎ1ℎ2(ℎ1−ℎ2)2 (ℎ1+ℎ2)4(𝑢1−𝑢2)4−𝑢1+𝑢2 2𝑑1. Notice that if 𝑢1=𝑢2or ℎ1=ℎ2, then  𝜆± ext, 𝜆± int equals the exact eigenvalues 𝜆± ext, 𝜆± int. Thus, we trivially have that 𝑃(𝜆)and 𝑃(𝑎)(𝜆)are the same polynomial. It can be also seen by noticing that 𝑑0=𝑑1= 0 in that case. In addition, in the case 𝑑0, 𝑑1≠0, the intersection between the two graphs is located at the point 𝜆0= −𝑑0∕𝑑1, which gives 𝜆0=𝑢1+𝑢2 2−ℎ1ℎ2(ℎ1−ℎ2) (ℎ1+ℎ2)3 (𝑢1−𝑢2)3 (2𝛾𝑐−𝑔(ℎ1+ℎ2)) .□ Let us introduce an interesting consequence of the third statement in previous theorem. Assuming that all the (approximated and exact) eigenvalues are real numbers, we are able to determine whether  𝜆− ext and  𝜆+ ext are lower and upper bounds of 𝜆− ext and 𝜆+ ext, respectively. It is not possible if complex eigenvalues appears, which will be addressed later. Otherwise, we have the following result: Corollary 1. Let 𝜆± ext, 𝜆± int ∈Rbe the roots of the exact characteristic polynomial 𝑃(𝜆)and  𝜆± ext, 𝜆± int ∈Rapproximations corresponding to the roots of 𝑃(𝑎)(𝜆). Let 𝜆0be the intersection point of 𝑃(𝜆)and 𝑃(𝑎)(𝜆). Then: 1. If 𝜆0∈( 𝜆− ext, 𝜆+ ext), then either  𝜆− ext < 𝜆− ext < 𝜆+ ext < 𝜆+ ext (in this case 𝑃( 𝜆− ext)>0and 𝑃( 𝜆+ ext)<0) or 𝜆− ext < 𝜆− ext < 𝜆+ ext < 𝜆+ ext (in this case 𝑃( 𝜆− ext)<0and 𝑃( 𝜆+ ext)>0). That is, one of the approximations  𝜆± ext is not a bound of the corresponding external eigenvalue. 2. If 𝜆0∉[ 𝜆− ext, 𝜆+ ext], then either  𝜆− ext < 𝜆− ext < 𝜆+ ext < 𝜆+ ext (in this case 𝑃( 𝜆± ext)>0) or 𝜆− ext < 𝜆− ext < 𝜆+ ext < 𝜆+ ext (in this case 𝑃( 𝜆± ext)<0). That is, either both approximations  𝜆± ext are bounds of the corresponding external eigenvalues or none of them verify it. 3. If 𝜆0= 𝜆− ext or 𝜆0= 𝜆+ ext, then either both approximations  𝜆± ext are bounds of the corresponding external eigenvalues or none of them verify it. As mentioned above, this corollary gives us a criteria to know whether the approximated eigenvalues are bounds of the exact ones. In fact, we see that in most of cases, at least one of them is not a bound, which would be desirable to have for a numerical scheme. However, we overcome this problem by redefining the approximation of external eigenvalues. It consists of performing a Newton–Raphson’s iteration starting from an appropriate value. Concretely, we follow these steps: Algorithm 1 (Bounds of external eigenvalues) 1: do 2:  𝜆+ ext ← 𝜆+ ext +𝑐ext 3: while 𝑃′( 𝜆+ ext)≤0 4: do 5:  𝜆− ext ← 𝜆− ext −𝑐ext 6: while 𝑃′( 𝜆− ext)≥0 7:  𝜆± ext ← 𝜆± ext −𝑃( 𝜆± ext)∕𝑃′( 𝜆± ext) Previous algorithm always gives as result a bound of the exact external eigenvalues, which is easy to prove by geometrical considerations. Therefore, in the case that the original approximation is not a bound of the exact eigenvalue, this algorithm corrects it and computes a bound. Otherwise, the algorithm improves the computed bound of the exact eigenvalue. In addition, its convergence is ensured by the fact that 𝑃(𝜆) and 𝑃(𝑎)(𝜆)have the same inflection points (since (𝑃−𝑃(𝑎))(𝜆)is a linear function). The inflection points are 𝛽±=𝑢1+𝑢2 2±√(𝑢1−𝑢2)2 12 +𝑔ℎ1+ℎ2 6, so the approximations  𝜆± ext given by (9) verify  𝜆− ext < 𝛽−< 𝛽+< 𝜆+ ext. In the case where complex eigenvalues arise, a particular numerical treatment will be considered, as explained in Section 3.2. In this case, where the characteristic polynomial (3) has exactly two real eigenvalues, we will need to use some values for the external eigenvalues when defining the numerical scheme. In addition, the approximations  𝜆± ext are not reliable in that case. Then, we will use the bounds of the external eigenvalues 𝜆± 𝑏𝑜𝑢𝑛𝑑 given in Laguerre (1880) for the roots of an arbitrary polynomial, that applied to the characteristic polynomial 𝑃(𝜆), defined by (3), gives 𝜆± 𝑏𝑜𝑢𝑛𝑑 =𝑢1+𝑢2 2± 3 √(𝑢1−𝑢2)2 12 +𝑔ℎ1+ℎ2 6.(11) Notice that 𝜆± 𝑏𝑜𝑢𝑛𝑑 are bounds of the roots of the polynomial only in the case of all real roots, which is not the case. However, they verify that 𝜆− 𝑏𝑜𝑢𝑛𝑑 < 𝛽−< 𝛽+< 𝜆+ 𝑏𝑜𝑢𝑛𝑑, as before. Then, applying Algorithm 1 starting from 𝜆± 𝑏𝑜𝑢𝑛𝑑 we get bounds of the external eigenvalues. Remark 1. Algorithm 1uses a single Newton’s iteration to get the bound of the external eigenvalues. However, one could make more iterations if a more accurate approximation of the eigenvalues is required. 2.2. Hyperbolicity region We study here the hyperbolicity region of the two-layer system (1). The main advantage of this new approach is the fact that it is based on the new approximation of the eigenvalues introduced in previous subsection, so it is not necessary to compute exactly whether we are in the non-hyperbolic region, as it is done in other works, e.g. Sarno et al. (2017) and Krvavica et al. (2018). Advances in Water Resources 182 (2023) 104587 5 M.J. Castro Díaz et al. In the following, we compare the exact and approximated hyperbolicity regions. When approximating this region, we use the eigenvalues approximation introduced in previous subsection and the classical (Schijf and Schöfled,1953). Concerning the later approximation, it estimates that complex eigenvalues arise when (6) is verified. That is, in the |𝑢1−𝑢2|− (ℎ1+ℎ2)plane the approximated hyperbolic region is the one verifying |𝑢1−𝑢2|> 𝐹SS(ℎ1, ℎ2),with 𝐹SS(ℎ1, ℎ2) = √𝑔(1 − 𝑟)(ℎ1+ℎ2). (12) Considering now the eigenvalues approximation (9), we see that 𝛾𝑐 (9d) is a non-negative value, and therefore only 𝑐int may become a complex value. This determines the approximated hyperbolicity region for this approximation of the eigenvalues. Concretely, we obtain complex eigenvalues if 𝑔ℎ1+ℎ2 2+ℎ1ℎ2 (ℎ1+ℎ2)2(𝑢1−𝑢2)2 <√(𝑔ℎ1+ℎ2 2)2 + (𝑟− 1)𝑔2ℎ1ℎ2+2𝑔ℎ1ℎ2 ℎ1+ℎ2 (𝑢1−𝑢2)2. This region cannot be described only in terms of |𝑢1−𝑢2|and (ℎ1+ℎ2), as it also depends explicitly on ℎ1, ℎ2. Thus, the hyperbolicity region in this case is determined by two boundary surfaces, which are defined in the space of axes ℎ1,ℎ2an |𝑢1−𝑢2|by |𝑢1−𝑢2|= (ℎ1+ℎ2)√ √ √ √𝑔(ℎ1+ℎ2) 2ℎ1ℎ2(1 ± √1 − 4(1 − 𝑟)ℎ1ℎ2 (ℎ1+ℎ2)2).(13) Fig. 2 shows a 3D representation of the limits of the hyperbolicity region for several values of 𝑟, namely, 𝑟= 0.3,0.5,0.9,0.99. The space between these two surface is the non-hyperbolic region. In order to compare this approximation with (12), we shall consider some sections or projections into the previous plane. In particular, for any given value ℎ1+ℎ2we consider uni-parametric family of heights given by ℎ1=𝛼(ℎ1+ℎ2), ℎ2= (1 − 𝛼)(ℎ1+ℎ2),for 𝛼∈ [0,1].(14) This allows us to write the boundaries of the hyperbolic region in terms of ||𝑢1−𝑢2||and (ℎ1+ℎ2)for any given value of 𝛼. More explicitly, the non-hyperbolic region in this case is the one determined, in the ||𝑢1−𝑢2||− (ℎ1+ℎ2)plane, by the curves (ℎ1+ℎ2) = 2𝛼(1 − 𝛼)|𝑢1−𝑢2|2 𝑔(1 ± √1 − 4(1 − 𝑟)𝛼(1 − 𝛼)) . That is, the non-hyperbolic region is the one determined by √𝑔(ℎ1+ℎ2)(1 − √1 − 4(1 − 𝑟)𝛼(1 − 𝛼)) 2𝛼(1 − 𝛼)<|𝑢1−𝑢2| <√𝑔(ℎ1+ℎ2)(1 + √1 − 4(1 − 𝑟)𝛼(1 − 𝛼)) 2𝛼(1 − 𝛼). This approximation of the hyperbolic region is analysed in Figs. 3–5, where we represent the exact non-hyperbolic region, the approximations given by the classical approximation (12) and the ones proposed here, for the cases 𝑟= 0.3,0.5,0.9,0.99 and 𝛼= 0.1,0.3,0.5in (14). It is worth mentioning that the cases 𝛼∈ (0.5,1) is symmetric to 𝛼∈ (0,0.5). We also plot in the figures the exact boundaries of the hyperbolic region in the case ℎ1=ℎ2, that are defined by √2𝑔(ℎ1+ℎ2)(1 − √𝑟)<||𝑢1−𝑢2||<√2𝑔(ℎ1+ℎ2)(1 + √𝑟)(15) Notice that the hyperbolic region defined by (15) coincides with the region deduced in Sarno et al. (2021) from an asymptotic analysis in the case ℎ1≈ℎ2. In these figures we see that the exact non-hyperbolic region is delimited by an upper and lower boundaries. However, the classical approximation only gives an upper boundary for the non-hyperbolic region, so all the right-hand side of the boundary is supposed to be non-hyperbolic. Then, the proposed approximation and the one corresponding to ℎ1=ℎ2notably improve the estimation of the non-hyperbolic region. We focus now on each limit of the non-hyperbolic region. We see that the upper limit for the non-hyperbolic region is better approached by (13), whereas the classical boundary overestimates this region (especially for 𝑟= 0.3) and the one in the case ℎ1=ℎ2is not a bound in some cases (see for instance upper figures in Fig. 3). Concerning the lower boundary, we see that the proposed approximation widely overestimates the non-hyperbolic region for some cases (see for instance Fig. 3), whereas the one for the case ℎ1=ℎ2is closer to the exact one. We see that the nearer the proportion between ℎ1and ℎ2is, the better approximation we obtain (compare Figs. 3 and 4). It is also noticeable that in the case 𝛼= 0.5(Fig. 5) the proposed approximation and the one with ℎ1=ℎ2coincide with the exact limits for the non-hyperbolic region, as stated in Theorem 1. As a conclusion from these qualitative analysis, it seems that the upper limit corresponding to the proposed approximation and the lower one corresponding to the case ℎ1=ℎ2bound the non-hyperbolic region (although we have not proven it formally), being a better approximation than taking both boundaries from the approximation of the eigenvalues (9). Then, we propose to consider this approximation for the practical point of view, i.e, 𝐹(𝐿) new(ℎ1, ℎ2)<|𝑢1−𝑢2|< 𝐹(𝑅) ℎ1=ℎ2(ℎ1, ℎ2)(16a) with 𝐹(𝐿∕𝑅) new (ℎ1, ℎ2) = √ √ √ √𝑔(ℎ1+ℎ2)3 2ℎ1ℎ2(1+(−∕+)√1 − 4(1 − 𝑟)ℎ1ℎ2 (ℎ1+ℎ2)2), (16b) and 𝐹(𝐿∕𝑅) ℎ1=ℎ2(ℎ1, ℎ2) = √2𝑔(ℎ1+ℎ2)(1+(−∕+)√𝑟).(16c) We remark that we have empirically checked that the approximated non-hyperbolic region defined by (16) contains the exact one in a wide variety of scenarios and configurations, including small values of the fluid depth (ℎ1+ℎ2). Nevertheless, we are not able to give a theoretical proof of such a result at this moment. We remark that model (1) becomes hyperbolic in the limit 𝑟→0, which is only reproduced by the case ℎ1=ℎ2. Both the classic and the proposed approximations widely overestimate the non-hyperbolic region in this limit case. 3. Numerical discretization of the two-layer system In order to design a numerical scheme for the two-layer system, several aspects have to be taken into account. First, it is well known that solutions of (1) may develop discontinuities, even for initial smooth data. Moreover, this system contains non-conservative products. So, we will consider numerical methods in the framework of path-conservative schemes introduced in Parés (2006). Second, situations with loss of hyperbolicity of system (1) may arise, as commented in previous section. Then, some numerical treatment will be necessary in practice. Here, we propose a second-order finite volume method with flux limiter, which is based on a suitable non-linear combination of the Lax–Wendroff method and a first-order path-conservative PVM scheme (see Castro Díaz and Fernández-Nieto (2012)). These PVM schemes define the numerical viscosity matrix as a polynomial evaluation of the Roe matrix (18). Concretely, we will locally switch between two different PVM methods, depending on where we are with respect to the hyperbolicity region, in such a way that both schemes reduce to the Lax–Wendroff method in case of smooth data. Once the numerical Advances in Water Resources 182 (2023) 104587 6 M.J. Castro Díaz et al. Fig. 2. Boundaries of the non-hyperbolic region defined by the proposed approximation of the eigenvalues on the plane ℎ1-ℎ2-|𝑢1−𝑢2|for 𝑟= 0.3, 0.5, 0.9 and 0.99. scheme is applied, a numerical treatment when complex eigenvalues arise is considered. In next subsections, we develop all the details of the numerical scheme, as well as the numerical treatment for non-hyperbolic situations. In addition, system (1) may be adapted to the submarine avalanche model introduced in Fernández-Nieto et al. (2008). This will further explained in Section 3.4, as well as the numerical discretization of such a model. 3.1. Numerical scheme We consider now the usual framework of finite volume methods: the domain is divided into control volumes 𝑉𝑖= [𝑥𝑖−1∕2, 𝑥𝑖+1∕2], for 𝑖∈. We denote by 𝑥𝑖the centre of each volume and, without loss of generality, a fixed mesh size 𝛥𝑥 =𝑥𝑖+1∕2 −𝑥𝑖−1∕2 is assumed here. In addition, for any variable 𝑓, we denote by 𝑓𝑛its approximation at the fixed time 𝑡𝑛. So, we denote the averaged conservative variable vector 𝒘𝑛 𝑖=1 𝛥𝑥 ∫𝑥𝑖+1∕2 𝑥𝑖−1∕2 𝒘(𝑡𝑛, 𝑥)𝑑𝑥. Finally, the time step 𝛥𝑡𝑛=𝑡𝑛−𝑡𝑛−1 is variable, and it is computed according to the Courant–Friedrichs–Lewy (CFL) condition for the stability of the scheme 𝛥𝑡𝑛≤𝛥𝑥 max𝑖∈|𝜆𝑛 𝑖| where |𝜆𝑖|= max{||𝑆𝐿,𝑖||,||𝑆𝑅,𝑖||}, being 𝑆𝐿,𝑖, 𝑆𝑅,𝑖 the minimum and maximum wave speeds at control volume 𝑉𝑖, respectively. We consider the framework of path-conservative schemes introduced in Parés (2006), which is extended to schemes based on polynomial reconstructions (see Gallardo et al. (2007) and Castro et al. (2017) among many others). First, let us extend the ideas of the hydrostatic reconstruction (Audusse et al.,2004) for the case of the two-layer system. To do so, we define (ℎ1+ℎ2)− 𝑖+1∕2 =𝑏𝑖+ℎ1,𝑖 +ℎ2,𝑖 −𝑏𝑖+1∕2, (ℎ1+ℎ2)+ 𝑖+1∕2 =𝑏𝑖+1 +ℎ1,𝑖+1 +ℎ2,𝑖+1 −𝑏𝑖+1∕2,(17a) and ℎ− 2,𝑖+1∕2 =𝑏𝑖+ℎ2,𝑖 −𝑏𝑖+1∕2 ℎ+ 2,𝑖+1∕2 =𝑏𝑖+1 +ℎ2,𝑖+1 −𝑏𝑖+1∕2,with 𝑏𝑖+1∕2 = max(𝑏𝑖, 𝑏𝑖+1).(17b) The reconstructed values of ℎ1are then defined as ℎ± 1,𝑖+1∕2 = (ℎ1+ℎ2)± 𝑖+1∕2 −ℎ± 2,𝑖+1∕2.(17c) The velocity values at the intercells correspond to the one at the corresponding centre of the cells, that is 𝑢± 𝑘,𝑖+1∕2 =𝑢𝑘,𝑖+1∕2±1∕2 for 𝑘= 1,2. Then, we define the reconstructed values to the left, 𝒘− 𝑖+1∕2, and right 𝒘+ 𝑖+1∕2 of the intercell as 𝒘± 𝑖+1∕2 = (ℎ± 1,𝑖+1∕2, 𝑞± 1,𝑖+1∕2, ℎ± 2,𝑖+1∕2, 𝑞± 2,𝑖+1∕2,)𝑇,(17d) where 𝑞± 𝑘,𝑖+1∕2 =ℎ± 𝑘,𝑖+1∕2𝑢± 𝑘,𝑖+1∕2 =ℎ± 𝑘,𝑖+1∕2𝑢𝑘,𝑖+1 2±1 2 for 𝑘= 1,2.(17e) We remark that the reconstructed states consider a continuous reconstructed bottom at the interfaces, which is typical for the hydrostatic Advances in Water Resources 182 (2023) 104587 7 M.J. Castro Díaz et al. Fig. 3. Projection of the non-hyperbolic region on the plane (ℎ1+ℎ2) − |𝑢1−𝑢2|for the case ℎ1= 0.1(ℎ1+ℎ2),ℎ2= 0.9(ℎ1+ℎ2)and 𝑟= 0.3,0.5,0.9,0.99 (upper left, upper right, lower left and lower right respectively). Comparison of exact region (red zone) and the approximated boundaries given by the classical (12) (dash-dotted green lines) and proposed (13) (dashed blue lines) boundaries for the hyperbolic region. Also the exact boundaries for the case ℎ1=ℎ2(15) (dotted black lines). reconstruction techniques and will result in the cancellation of the geometric source contributions at the interfaces. Now, we first define a generalized path-conservative Roe matrix of 𝐴(𝒘)for the particular choice of segment paths, that is, we set 𝑖+1∕2 =⎛⎜⎜⎜⎜⎝ 0 1 0 0 𝑔ℎ1,𝑖+1∕2 −𝑢2 1,𝑖+1∕2 2𝑢1,𝑖+1∕2 𝑔ℎ1,𝑖+1∕2 0 0 0 0 1 𝑔𝑟ℎ2𝑖+1∕2 0𝑔ℎ2,𝑖+1∕2 −𝑢2 2,𝑖+1∕2 2𝑢2,𝑖+1∕2 ⎞⎟⎟⎟⎟⎠(18) where ℎ𝑙,𝑖+1∕2 = ℎ− 𝑙,𝑖+1∕2 +ℎ+ 𝑙,𝑖+1∕2 2, 𝑢𝑙,𝑖+1∕2 =√ℎ− 𝑙,𝑖+1∕2𝑢− 𝑙,𝑖+1∕2 +√ℎ+ 𝑙,𝑖+1∕2𝑢+ 𝑙,𝑖+1∕2 √ℎ− 𝑙,𝑖+1∕2 +√ℎ+ 𝑙,𝑖+1∕2 ,for 𝑙= 1,2.(19) This matrix satisfies that 𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2) = 𝐹𝑐(𝒘+ 𝑖+1∕2)−𝐹𝑐(𝒘− 𝑖+1∕2)+𝐵𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2), where 𝐹𝑐(𝒘) = ⎛⎜⎜⎜⎜⎝ 𝑞1 𝑞2 1∕ℎ1 𝑞2 𝑞2 2∕ℎ2 ⎞⎟⎟⎟⎟⎠ and 𝐵𝑖+1∕2 =⎛⎜⎜⎜⎜⎝ 0 0 0 0 𝑔ℎ1,𝑖+1∕2 0𝑔ℎ1,𝑖+1∕2 0 0 0 0 0 𝑔𝑟ℎ2𝑖+1∕2 0𝑔ℎ2,𝑖+1∕2 0 ⎞⎟⎟⎟⎟⎠ . Thus, the discrete form of the system (2) will take the form 𝒘𝑛+1 𝑖=𝒘𝑛 𝑖−𝛥𝑡𝑛 𝛥𝑥 (𝑛 𝑖+1∕2 −𝑛 𝑖−1∕2)+(𝛥𝑡𝑛 2𝛥𝑥 )2(𝑛 𝑖−1∕2 +𝑛 𝑖+1∕2).(20) The intercell contributions 𝑖+1∕2 are defined within the framework of PVM schemes (Castro Díaz and Fernández-Nieto,2012) and based on previous reconstructed states: 𝑖+1∕2 =1 2(𝐹𝑐(𝒘+ 𝑖+∕2) + 𝐹𝑐(𝒘− 𝑖+∕2) + 𝐵𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2) − 𝑖+1∕2), where 𝑖+1∕2 =𝛼𝑖+1∕2 0(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2) +𝛼𝑖+1∕2 1(𝐹𝑐(𝒘+ 𝑖+∕2) − 𝐹𝑐(𝒘+ 𝑖+∕2) + 𝐵𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2)) +𝛼𝑖+1∕2 2𝑖+1∕2(𝐹𝑐(𝒘+ 𝑖+∕2) − 𝐹𝑐(𝒘+ 𝑖+∕2) + 𝐵𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2)),(21) with 𝛼𝑖+1∕2 𝑘,𝑘= 0,1,2the coefficients of the selected (second-order) polynomial that define the chosen PVM scheme. Denoting by 𝑖+1∕2 =𝐵𝑖+1∕2(𝒘+ 𝑖+1∕2 −𝒘− 𝑖+1∕2)the nonconservative products, they write 𝑖+1∕2 =⎛⎜⎜⎜⎜⎜⎝ 0 𝑔ℎ1,𝑖+1∕2 (𝛥𝜂1)𝑖+1∕2 0 𝑔ℎ2,𝑖+1∕2(𝑟(𝛥𝜂1)𝑖+1∕2 + (1 − 𝑟)(𝛥𝜂2)𝑖+1∕2)⎞⎟⎟⎟⎟⎟⎠ ,(22) being (𝛥𝜂1)𝑖+1∕2 = (ℎ1+ℎ2)+ 𝑖+1∕2 − (ℎ1+ℎ2)− 𝑖+1∕2, Advances in Water Resources 182 (2023) 104587 8 M.J. Castro Díaz et al. Fig. 4. Projection of the non-hyperbolic region on the plane (ℎ1+ℎ2) − |𝑢1−𝑢2|for the case ℎ1= 0.3(ℎ1+ℎ2),ℎ2= 0.7(ℎ1+ℎ2)and 𝑟= 0.3,0.5,0.9,0.99 (upper left, upper right, lower left and lower right respectively). Comparison of exact region (red zone) and the approximated boundaries given by the classical (12) (dash-dotted green lines) and proposed (13) (dashed blue lines) boundaries for the hyperbolic region. Also the exact boundaries for the case ℎ1=ℎ2(15) (dotted black lines). (𝛥𝜂2)𝑖+1∕2 =ℎ+ 2,𝑖+1∕2 −ℎ− 2,𝑖+1∕2.(23) Finally, 𝑖±1∕2 are second-order corrections corresponding to the Lax– Wendroff method in case of non-conservative systems (see Castro et al. (2009b)), that will be described later (see (30)), being 𝑖±1∕2 = 0 for first order finite volume methods. Remark 2. Notice that the values of (ℎ1+ℎ2)± 𝑖+1∕2,ℎ± 2,𝑖+1∕2 could become negative in wet/dry fronts, as the case of an emerging bottom. With the goal of dealing properly these wet/dry situations, in practice we use (ℎ1+ℎ2)± 𝑖+1∕2 = max ((ℎ1+ℎ2)± 𝑖+1∕2,0), ℎ± 2,𝑖+1∕2 = max (ℎ± 2,𝑖+1∕2,0), which preserves the nonnegativity of (ℎ1+ℎ2)and ℎ2. Let us remark that this assures the nonnegativity of ℎ1as well. Note that the definition of coefficients 𝛼𝑖+1∕2 𝑘will depend on the choice of the PVM numerical solver. These coefficients actually depend on the approximation of the eigenvalues of the system. Here we will use IFCP solver (see Fernández-Nieto et al. (2011)) and PVM-2U (Castro Díaz and Fernández-Nieto,2012) as first order FV solvers. Both schemes are defined in terms of quadratic viscosity matrices. On the one hand, IFCP method uses information on the external and internal waves, producing accurate results for the two-layer shallow water system. On the other hand, PVM-2U method only uses information on the external waves, in a similar way as HLL solver does. Nevertheless, PVM-2U is less diffusive than HLL solver when applied to the two-layer shallow water system (see Castro Díaz and Fernández-Nieto (2012)). Moreover, in order to achieve second-order accuracy, both schemes will be combined with the Lax–Wendroff method through a flux-limiter approach. In what follows we describe both approaches as well as the flux-limiter approach interpreted as PVM solver. 3.1.1. IFCP-FL solver In order to reach second order in space an time in smooth areas, an extension of the IFCP solver based on flux-limiter functions is proposed in this section. The IFCP method introduced in Fernández-Nieto et al. (2011) was specially designed for two-layer systems, and the numerical tests presented there showed that IFCP has an accuracy similar to Roe scheme, being around 3 times faster. Concerning the well-balanced property of the scheme, it is straightforward to check that IFCP scheme is exactly well-balanced for the water at rest solutions 𝑞1=𝑞2= 0, 𝑏 +ℎ1+ℎ2=𝑐𝑠𝑡, and 𝑏+ℎ2=𝑐𝑠𝑡. In addition, this scheme is 𝐿∞-stable with the usual CFL condition (see Fernández-Nieto et al. (2011) for details). It is defined by a parabolic viscosity matrix method that uses the information of the intermediate fields. In particular, this method uses the approximation of the slower (𝑆𝑖+1∕2 𝐿) and faster (𝑆𝑖+1∕2 𝑅) waves of the Riemann problem associated to intercell 𝑥𝑖+1∕2 combined with an intermediate wave speed 𝑆𝑖+1∕2 𝑖𝑛𝑡 . Here we use 𝑆𝑖+1∕2 𝐿= min(𝑆𝐿,𝑖+1∕2, 𝑆− 𝐿,𝑖+1∕2), 𝑆𝑖+1∕2 𝑅= max(𝑆𝑅,𝑖+1∕2, 𝑆+ 𝑅,𝑖+1∕2). being 𝑆𝐿,𝑘, 𝑆𝑅,𝑘 the output of Algorithm 1taking as data the states 𝒘𝑖(𝑘=𝑖),𝒘𝑖+1 (𝑘=𝑖+ 1) or the Roe state (19) (𝑘=𝑖+ 1∕2). The intermediate wave speed is defined as 𝑆𝑖+1∕2 𝑖𝑛𝑡 =𝒔𝑖+1∕2 max (| 𝜆− 𝑖𝑛𝑡,𝑖+1∕2|,| 𝜆+ 𝑖𝑛𝑡,𝑖+1∕2|) Advances in Water Resources 182 (2023) 104587 9 M.J. Castro Díaz et al. Fig. 5. Projection of the non-hyperbolic region on the plane (ℎ1+ℎ2) − |𝑢1−𝑢2|for the case ℎ1=ℎ2= 0.5(ℎ1+ℎ2)and 𝑟= 0.3,0.5,0.9,0.99 (upper left, upper right, lower left and lower right respectively). Comparison of exact region (red zone) and the approximated boundaries given by the classical (12) (dash-dotted green lines) and proposed (13) (dashed blue lines) boundaries for the hyperbolic region. Also the exact boundaries for the case ℎ1=ℎ2(15) (dotted black lines). with 𝒔𝑖+1∕2 ={sign(𝑆𝑖+1∕2 𝐿)if |𝑆𝑖+1∕2 𝐿|≥|𝑆𝑖+1∕2 𝑅|, sign(𝑆𝑖+1∕2 𝑅)otherwise. Then, the viscosity matrix 𝑄𝑖+1∕2 associated to this IFCP method is the evaluation of the matrix 𝑖+1∕2 by the second order polynomial 𝑝IFCP 𝑖+1∕2(𝑥) = 𝛼IFCP 0,𝑖+1∕2 +𝛼IFCP 1,𝑖+1∕2𝑥+𝛼IFCP 2,𝑖+1∕2𝑥2,(24) in such a way that the polynomial 𝑝IFCP 𝑖+1∕2(𝑥)satisfies that 𝑝IFCP 𝑖+1∕2(𝑆𝑖+1∕2 𝐿) = |𝑆𝑖+1∕2 𝐿|, 𝑝IFCP 𝑖+1∕2(𝑆𝑖+1∕2 𝑅) = |𝑆𝑖+1∕2 𝑅| and 𝑝IFCP 𝑖+1∕2(𝑆𝑖+1∕2 𝑖𝑛𝑡 ) = |𝑆𝑖+1∕2 𝑖𝑛𝑡 |. That is, the coefficients 𝛼IFCP 𝑙,𝑙= 0,1,2are the unique solution of the linear system ⎛⎜⎜⎜⎝ 1𝑆𝐿(𝑆𝐿)2 1𝑆𝑅(𝑆𝑅)2 1𝑆𝑖𝑛𝑡 (𝑆𝑖𝑛𝑡)2⎞⎟⎟⎟⎠⎛⎜⎜⎜⎝ 𝛼IFCP 0 𝛼IFCP 1 𝛼IFCP 2 ⎞⎟⎟⎟⎠ =⎛⎜⎜⎜⎝|𝑆𝐿| |𝑆𝑅| |𝑆𝑖𝑛𝑡|⎞⎟⎟⎟⎠ .(25) After some straightforward computations, we get 𝛼IFCP 0=𝛿𝐿𝑆𝑅𝑆𝑖𝑛𝑡 +𝛿𝑅𝑆𝐿𝑆𝑖𝑛𝑡 +𝛿𝑖𝑛𝑡 𝑆𝐿𝑆𝑅, 𝛼IFCP 1= −𝑆𝐿(𝛿𝑅+𝛿𝑖𝑛𝑡) − 𝑆𝑅(𝛿𝐿+𝛿𝑖𝑛𝑡) − 𝑆𝑖𝑛𝑡(𝛿𝐿+𝛿𝑅), 𝛼IFCP 2=𝛿𝐿+𝛿𝑅+𝛿𝑖𝑛𝑡, (26) with 𝛿𝐿=|𝑆𝐿| (𝑆𝐿−𝑆𝑅)(𝑆𝐿−𝑆𝑖𝑛𝑡), 𝛿𝑖𝑛𝑡 =|𝑆𝑖𝑛𝑡| (𝑆𝑖𝑛𝑡 −𝑆𝐿)(𝑆𝑖𝑛𝑡 −𝑆𝑅), 𝛿𝑅=|𝑆𝑅| (𝑆𝑅−𝑆𝐿)(𝑆𝑅−𝑆𝑖𝑛𝑡).(27) Remark 3. In Section 4it will be shown that eventually, for some numerical test cases, a lack of numerical viscosity of the IFCP scheme results into a small stationary artificial shock. This is associated to very small values of the intermediate eigenvalue 𝑆𝑖𝑛𝑡. In order to avoid this problem, one may consider the Harten’s regularization (see Harten (1997)) for the intermediate eigenvalue. Then, |𝑆𝑖𝑛𝑡|is replaced by |𝑆𝑖𝑛𝑡|+1 + sign(𝜖−|𝑆𝑖𝑛𝑡|) 2(𝜖2+𝑆2 𝑖𝑛𝑡 2𝜖−|𝑆𝑖𝑛𝑡|) in the computation of 𝛿𝑖𝑛𝑡 in (27), where 𝜖is a small parameter. In practice, we use 𝜖= 0.25𝑐ext with 𝑐ext defined by (9b). We remark that other entropy-fix techniques could be employed. The goal here is to obtain a second order method as combination of the IFCP method and the Lax–Wendroff method. The resulting scheme will reduce to IFCP method in non-smooth regions and to Lax–Wendroff in smooth areas, achieving second order accuracy in smooth regions. To this aim, first the viscosity matrix of IFCP and Lax–Wendroff schemes are defined by means of a polynomial evaluation: 𝑝IFCP 𝑖+1∕2(𝑥)is given by (24) and the one for the Lax–Wendroff method is 𝑝LW 𝑖+1∕2(𝑥) = 𝛥𝑡 𝛥𝑥 𝑥2. So, both viscosity matrices are the evaluation of the corresponding polynomial on the Roe matrix. To combine them, we consider the polynomial associated to the IFCP-FL method: 𝑝IFCP-FL 𝑖+1∕2 (𝑥) = 𝛼IFCP-FL 0,𝑖+1∕2 +𝛼IFCP-FL 1,𝑖+1∕2 𝑥+𝛼IFCP-FL 2,𝑖+1∕2 𝑥2,(28) Advances in Water Resources 182 (2023) 104587 16 M.J. Castro Díaz et al. Fig. 10. Test 2: Interface 𝜂2between the fluids (left-hand side) and velocity of the lower layer 𝑢2(right-hand side) at times 𝑡= 0.3,0.5,1.5s, computed with the IFCP and PVM-2U methods, and their -FL versions. Inner figures in the right-hand side figures are zooms around 𝑥= 0.5. Fig. 11. Test 2: Interface 𝜂2between the fluids at times 𝑡= 0.5,1.5s, computed with the IFCP-FL method and also considering the Harten’s regularization for the intermediate eigenvalue. that only the IFCP-FL scheme is able to reproduce accurately the twosteps shape of the interface and the velocity near 𝑥= ±0.5at 𝑡= 0.3 s. For larger times, we see that the PVM-2U-FL and IFCP methods also reproduce well the reference solution, being the IFCP-FL method more accurate. Interestingly, the IFCP method produces similar results to the PVM-2U-FL method despite of being first-order accurate, still with qualitatively good accuracy, in the sense that it is able to reproduce the shape of the reference solution, in contrast to the PVM-2U scheme. Here we observe that the ICFP method is much more accurate than the PVM-2U scheme. Notice that in this test we have two shock travelling at a very small velocity. In particular, this makes that the internal wave velocity 𝑆𝑖𝑛𝑡 is very small. In such cases, the IFCP-FL method may develop spurious oscillation as result of the lack of numerical viscosity, similarly to what it is observed for the Roe’s method close to sonic points. To avoid this problem, we could use the Harten’s regularization (see Remark 3) for the computation of the internal eigenvalue. In Fig. 11 we show the results of the IFCP-FL method with/without the regularization. We see some oscillations near the advancing fronts, which does not disappear neither at larger times nor refining the spatial mesh. This behaviour will also be seen in Test 4 (Section 4.2.3). Concerning the efficiency Advances in Water Resources 182 (2023) 104587 17 M.J. Castro Díaz et al. Fig. 12. Test 3: Interfaces 𝜂1, 𝜂2(left-hand side) and velocities 𝑢1, 𝑢2(right-hand side) at times 𝑡= 0.2,0.5s in the case ||𝑢1−𝑢2||= 9 (hyperbolic initial condition), computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT) and (CHMT). of the proposed IFCP-FL method, we have measured the computational time needed by the IFCP and its -FL version for this test. In the current configuration, both methods require less than one second to complete the simulation. Then, in order to compare them, we have set the final time 𝑡𝑓= 10 s, and a finer discretization is used 𝛥𝑥 = 2.5 × 10−3 (domain [−3,3],2400 cells). The IFCP method spends 7.91 s, whereas 8.63 s are needed by the IFCP-FL method. As it is observed, the increase of the computational time when considering the -FL version of the method is 9% approximately. This computational times are similar in Sections Section 4.2.2,4.2.3,4.3.1, and 4.3.2. In what follows, we will focus on the IFCP-FL scheme, unless otherwise specified. 4.2.2. Test 3: hyperbolic conditions with 𝐹(𝑅)<||𝑢1−𝑢2|| As it has been explained, the non-hyperbolic region has a left (𝐹(𝐿)) and a right (𝐹(𝑅)) boundary. One of the main advantages of the proposed approximation with respect to the classic one is the fact that it is possible to approximate 𝐹(𝑅). We consider now the same initial conditions for the heights given in (39) and take 𝑢1(𝑥, 0) = 4, 𝑢2(𝑥, 0) = −5. Then ||𝑢1−𝑢2||= 9 and we are in the hyperbolic region beyond 𝐹(𝑅), that is, in zone (IV) in Fig. 7. Here no treatment is needed since the system is hyperbolic. However, treatment (CFT) applies the friction correction after the first time iteration by reducing the velocity difference to go to zone (I), since it estimates that we are in a nonhyperbolic region. On the contrary, treatments (NFT) and (CHMT) detect that the regime is hyperbolic and no correction is applied. Fig. 12 shows the results in this case for the interfaces 𝜂1, 𝜂2and the velocities 𝑢1, 𝑢2. We see that very different results are obtained when applying the treatment (CFT) whereas no corrections are made with (NFT) and (CHMT). In particular, we see how the velocities are modified to dramatically reduce ||𝑢1−𝑢2||after the first time iteration. It even changes the sign of the velocity 𝑢1, making this treatment fully useless in these situations. 4.2.3. Test 4: non-hyperbolic initial conditions 𝐹(𝐿)<||𝑢1−𝑢2||< 𝐹(𝑅) We consider now the same initial conditions for the heights given in (39), but with a velocity difference ||𝑢1−𝑢2||in order to fall into the non-hyperbolic region. Concretely, we consider ||𝑢1−𝑢2||= 0.5with 𝑢1(𝑥, 0) = 0.2,𝑢2(𝑥, 0) = −0.3. In this case, a correction is needed in order to obtain meaningful results, as shown in Fig. 13. In this figure, we show the solution when applying the treatment for the loss of hyperbolicity explained in Section 3.2, and also without considering any treatment. In the latter case, we see spurious perturbations growing in time, in contrast to the case where we consider corrections (CFT), (NFT) or (CHMT). On the one hand, there is no significant differences between corrections (CFT) and (NFT) in this case, since the left boundaries of the non-hyperbolic region (𝐹𝑆𝑆 and 𝐹(𝐿) new, respectively) almost coincide for 𝑟= 0.99 (see Figs. 3–5). In this case, the friction correction moves the flow to zone (I) after a single time iteration, and these corrections are activated again if the flow goes to zone (II) in any iteration. Usually, they act just in the first time iteration for this flow configuration. In order to see more differences between corrections (CFT) and (NFT), we should consider a lower density ratio, as we will see in Section 4.3. On the other hand, correction (CHMT), which changes to a simplified hyperbolic two-layer system but does not modify the hyperbolicity zone (the flow is in zone (II) during all the simulation), widely differs from the ones applying friction. When looking at the velocities, we see that their magnitudes are not modified, in contrast to corrections based on adding friction. Thus, this correction (CHMT) is an interesting alternative since it avoids the numerical problems related to the loss of hyperbolicity and does not significantly change the magnitude of the variables (specially the velocities). In this test, when the treatments (CFT) or (NFT) are applied, the velocity difference is reduced after a single iteration to fall in the hyperbolic region. Then, we have again a shock travelling at a very slow velocity. In Fig. 14 we observe that it makes appear some spurious oscillations in the shock in all the variables related to the (almost) vanishing intermediate eigenvalue 𝑆𝑖𝑛𝑡. As it has been previously commented in Test 2, it could be avoid by using the Harten’s regularization. We also show how these oscillations does not disappear when reducing the mesh step to 𝛥𝑥∕4. 4.2.4. Test 5: transcritical flow with barotropic tidal forcing We consider now a classical test case consisting on a transcritical flow over a topography. In particular, we consider the configuration Advances in Water Resources 182 (2023) 104587 18 M.J. Castro Díaz et al. Fig. 13. Test 4: Interfaces 𝜂1, 𝜂2(left-hand side) and velocities 𝑢1, 𝑢2(right-hand side) at times 𝑡= 0.5,1.0,1.5s in the case ||𝑢1−𝑢2||= 5 (non-hyperbolic initial condition), computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT) and with no correction (no treat.). Fig. 14. Test 4: 𝜂1, 𝑢1, 𝜂2, 𝑢2at time 𝑡= 1.5s computed with the IFCP-FL method with the treatment (NFT) with mesh steps 𝛥𝑥 and 𝛥𝑥∕4, and with the Harten’s regularization for the internal eigenvalue (mesh steps 𝛥𝑥). Inner figures are zoom. Advances in Water Resources 182 (2023) 104587 19 M.J. Castro Díaz et al. Fig. 15. Test 5: 𝜂2and hyperbolicity zones at times 𝑡= 200,2200,2330 s computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT), and with no correction (no treat.). showed in Castro-Díaz et al. (2010). The computational domain is [0,10] and the initial conditions are 𝑞1(𝑥, 0) = 𝑞2(𝑥, 0) = 0, and 𝑏(𝑥)=0.5𝑒−(𝑥−5)2− 1, ℎ1(𝑥, 0) = {0.5if 𝑥 < 5, 0otherwise; 𝑏(𝑥) + ℎ1(𝑥, 0) + ℎ2(𝑥, 0) = 0. This initial condition corresponds to an internal dam break problem. Free-outflow boundary conditions are assumed. The upper layer then flows from the left to the right (𝑢1>0) while the opposite happens for the lower layer (𝑢2<0). After some time, a steady state is reached, which exhibits an internal hydraulic jump (see Fig. 15 (left)). In this figure, we see that different steady states, with different location for the internal shock, depending on the considered treatment for the loss of hyperbolicity. We see that the solution reached with the treatment (CHMT) is almost identical to the case of no applying a treatment, whereas the treatments based on adding friction ((NFT) and (CFT)) move the location of the shock on the interface. We see that around 𝑥= 4.5the flow goes into zone (II) in Fig. 7, where the treatments are applied. Notice that (NFT) and (CFT) produce identical solutions since the left boundaries (𝐹𝑆𝑆 and 𝐹(𝐿) new) of the non-hyperbolic region are quite close for this case (𝑟= 0.99), as it has been previously commented. Here we also consider a barotropic tidal force at the right boundary for 𝑡 > 200 s. To this end, we consider free-outflow boundary conditions for all the variables, except for ℎ1. For this variable we take ℎ1(10, 𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ free-outflow if 𝑡≤200, ℎ1(10,200) + ℎ2(10,200) +0.005 cos (𝑡∕50 + 𝜋∕4)otherwise; Thus, we obtain an oscillatory behaviour of the solution. In Fig. 15 (right) we show the interface 𝜂2at two different times, together with the hyperbolic zones for each case. Again, the solution when considering the treatment (CHMT) is almost identical to not applying a treatment, in contrast to the solution when a friction treatment is applied. In Fig. 16 we analyse the time evolution of the interface at locations 𝑥= 4.15,4.45,6. We see differences between the results with the friction treatments and (CHMT). We see that the model with (CHMT) is always closer to the solution without any correction. It is even more evident at 𝑥= 4.15. Focusing now on the evolution at 𝑥= 4.45, Fig. 16b shows a zoom, where we see that the model without correction produces small oscillation of the solution, whereas it is smother with the (CHMT). Remark that this spurious oscillations are small and they do not affect neither the stability of the scheme nor the appearance of bigger oscillations. Fig. 17 shows the evolution on time and space of the variable 𝑧=𝐹(𝐿) new∕||𝑢1−𝑢2||. Notice that values 𝑧 > 1represent zones in the hyperbolic region whereas the values 𝑧 < 1correspond to points falling into the non-hyperbolic region. Notice that the plotted cyan and green curves are the isolines corresponding to the points over the boundary 𝐹(𝐿) new =||𝑢1−𝑢2||. That is, the closed region by these curves corresponds to states falling into the non-hyperbolic region. We observe that the results with the correction (CHMT) and with no correction exhibits a similar behaviour of the solution both in the steady and the tidal periods. However, values with 𝐹(𝐿) new∕||𝑢1−𝑢2||≤1do not appear when applying a friction treatment. In general, we can conclude that for this case, although the results with the treatment (CHMT) are much closer to the those results obtained when no treatment is applied, all the treatments for the loss of hyperbolicity could be used, since all the results are similar. We remark that for this test, we have not observed relevant spurious oscillations on the solutions when no treatment is considered. In this case, the scheme described in Section 3.1 is able to deal with this loss of hyperbolicity thanks to the better approximation of the eigenvalues bounds, and the fact of switching to the PVM-2U-FL scheme in zone (II). 4.3. Fluids with low density ratio 𝑟= 0.3 In previous tests, we have shown the relevance of estimating the right boundary of the non-hyperbolic region 𝐹(𝑅)(see Fig. 12). However, there is no significant differences in previous test cases for treatments on the left boundary of (CFT) and (NFT), since the boundaries 𝐹𝑆𝑆 and 𝐹(𝐿) new are very close for 𝑟≈ 1. We perform here some tests consisting of two fluids with 𝑟= 0.3. In this case we expect to have differences between (CFT) and (NFT), where the left boundaries of the non-hyperbolic boundaries differs (see Fig. 5). 4.3.1. Test 6: 𝑟= 0.3, non-hyperbolic initial conditions 𝐹(𝐿)<||𝑢1−𝑢2||< 𝐹(𝑅) Let us consider the initial conditions ℎ2(𝑥, 0) = {0.4if |𝑥|<0.5 0.3otherwise , ℎ1(𝑥, 0) = 0.6 − ℎ2(𝑥, 0), 𝑢1(𝑥, 0) = −1.5, 𝑢2(𝑥, 0) = 1.5. Here we use again 𝛥𝑥 = 0.01 and 𝑏(𝑥)=0. These initial conditions fall into the non-hyperbolic region. Therefore, treatments (CFT) and (NFT) will reduce ||𝑢1−𝑢2||after a single iteration to move the velocity difference to the hyperbolic region configuration. In this sense, the behaviour is similar to previous test, where the treatments (CFT) and (NFT) act only in the first time iteration in order to go to zone (I), whereas treatment (CHMT) does not modify the hyperbolicity zone in which we are and it will be zone (II) for all the simulation time. We see in Fig. 18 that the (CFT) applies a bigger friction force than (NFT) in order to reduce the velocity difference, as excepted. It is clear looking at the layer velocities, and it leads to different results. We also see that meaningless results are obtained if no treatment is considered. A very different solution is observed by changing the model (CHMT), especially for the first layer. We see that the free surface is creating a discontinuous profile in the location of the first step in the interface, which remind us to the hydraulic jump that appears in a transcritical flow over a bottom. In Fig. 19 we see how applying a treatment for the loss of hyperbolicity based on adding friction dramatically change Advances in Water Resources 182 (2023) 104587 20 M.J. Castro Díaz et al. Fig. 16. Test 5: 𝜂2evolution at locations 𝑥𝑖= 4.15,4.5,6m computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT), and with no correction (no treat.). the flow regime in this test. For instance, focusing on the upper layer (𝑖= 1), we see that at initial time the internal Froude number is always greater than 1 (𝐹𝑟1>1) and it changes to 𝐹 𝑟1<1after a single iteration. It makes that the model reproduces a totally different solution. This is an important drawback of adding friction to deal with the loss of hyperbolicity, since we lose these particular solutions. In this sense, the treatment based on changing the model (CHMT) does not change the flow regime, even though some information is lost due to the simplifications made to deal with a hyperbolic model. We want to show now that the solution obtained by changing the model is indeed a reasonable solution. The discontinuity observed at the free surface is related to the almost stationary discontinuity at the interface 𝜂2. It has a similar behaviour to what is obtained by the Shallow Water model with a bottom topography given by 𝑏(𝑥) = 𝜂2(𝑥, 0). In Fig. 20 we see this comparison for 𝑡= 0.3and 𝑡= 1.5s. 4.3.2. Test 7: 𝑟= 0.3, hyperbolic initial conditions 𝐹(𝑅)<||𝑢1−𝑢2|| We consider now a configuration where both internal Froude number are greater than 1 along the domain and hyperbolic initial condition such that 𝐹(𝑅)<||𝑢1−𝑢2||. That is, we are in zone (III) or (IV) according to Fig. 7. To this aim, we consider ℎ2(𝑥, 0) = {0.4if |𝑥|<0.5 0.3otherwise , ℎ1(𝑥, 0) = 1.0 − ℎ2(𝑥, 0), 𝑢1(𝑥, 0) = −3, 𝑢2(𝑥, 0) = 3. Notice that if ||𝑢1−𝑢2||is large enough to be into the hyperbolic region (zones (III) and (IV)) throughout the time iterations, we would obtain similar conclusions to Test 3 4.2.2. That is, no correction is made when considering treatments (NFT) and (CHMT) whereas the velocities are reduced if (CFT) is used. Then, for this test we choose a velocity difference ||𝑢1−𝑢2||= 6 that falls into zones (III) or (IV) (hyperbolic region), but it is not far from the right boundary of the non-hyperbolic region 𝐹(𝑅). Then, when the flow dynamics starts, some of the cells of the domain will come into the non-hyperbolic region and the corresponding correction (NFT or CHMT) will be applied. However, when considering (CFT) the whole velocity field will be modified in order to get ||𝑢1−𝑢2||< 𝐹(𝐿), since the classic treatment does not consider the right boundary of the non-hyperbolic region 𝐹(𝑅). We see this behaviour in Fig. 21, where we plot the free surface/interface, layer velocities, and also the zone in which each cell is according to Fig. 7, at times 𝑡= 0.05 and 0.5 s. Let us describe the results for each treatment of the loss of hyperbolicity: •No treatment: spurious oscillations appear as in previous tests due to the lack of hyperbolicity. •(CFT): after a single iteration, the velocity difference is reduced in the whole domain to achieve ||𝑢1−𝑢2||< 𝐹(𝐿), as excepted. It is done even in those areas where the system is hyperbolic. •(CHMT): In those areas where the system becomes non-hyperbolic (zone (II): 𝐹(𝐿)<||𝑢1−𝑢2||< 𝐹(𝑅)) we change to the hyperbolic model in Section 3.2.2. For instance, for 𝑡= 0.5s it is done just for those cells 𝑉𝑖with 0.37 < 𝑥𝑖<1.04. In the rest of the domain no treatment is needed. •(NFT): Once a cell moves into zone (II), the velocity difference is reduced in this cell and not in the rest of the domain. This creates arbitrary shocks in those areas, so that results with no physical meaning are obtained after a few time iterations. As we observe in Fig. 21, they are even worse than those where no treatment is considered. Finally, notice that similar conclusions for this test would be obtained when considering the procedures and techniques introduced in previous works (Sarno et al.,2017; Krvavica et al.,2018). According to these results, we can conclude that applying a treatment based on adding friction and taking into account the right boundary of the non-hyperbolic region (𝐹(𝑅)) is not appropriate in practical situations where a large velocity difference may occur. Let us remark that situations where the flow goes from zone (I) towards (III,IV), and after some time it returns to zone (I) (through zone (II)) arise in many practical applications of two-layer systems. This is the case, for instance, of a granular mass falling into the water, i.e., submarine avalanches. We deal with this kind of geophysical flows in next subsection. Before, a reasonable question is whether one could consider in the whole domain the hyperbolic simplified model obtained when replacing the second momentum equation in (1) by (34), whose eigenvalues (35) are real values in any situation. That is, instead of considering a treatment for the loss of hyperbolicity, one could use the simplified hyperbolic two-layer model given by (34). In principle, the results should be close as long as the hypothesis 𝜕𝑥𝜂1∼(𝜀)holds. In Fig. 22 Advances in Water Resources 182 (2023) 104587 21 M.J. Castro Díaz et al. Fig. 17. Test 5: 𝑧=𝐹(𝐿) new∕||𝑢1−𝑢2||(𝑧 > 1zone (I), 𝑧 < 1zone (II), in Fig. 7) evolution (𝑡∈ [0,200] left-hand side, 𝑡∈ [2000,5000] right-hand side) computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (no treatment) upper row and (NFT) lower row. Cyan (respectively green) lines are the isolines corresponding to 𝐹(𝐿) new =||𝑢1−𝑢2||computed with no correction (respectively (CHMT)). we show this comparison for times 𝑡= 0.5,1.5s. We see that, in this case, the simplified hyperbolic model produces results that are far from those of the two-layer model with treatment (CHMT). Then, it is not possible to use the simplified model when variations of the free surface are not small. This comparison will be also discussed in the next subsections. 4.4. Application to landslides In this subsection we consider the application of two-layer systems to simulate submarine avalanches. To this aim, we consider the twolayer model described in Section 3.4. We first perform a test consisting of a granular mass falling into a fluid, including a non-flat bathymetry with wet/dry areas, and strong transitions where the changes on the hyperbolicity zone arise. Later, we consider a submarine landslide where the fluid is flowing under a subcritical regime. 4.4.1. Test 8: landslide over a channel In this test, we consider a granular mass that flows downhill of a mountain towards the water. That is, an aerial landslide that falls into a mass of water such as a lake, a river or a coastal area. Actually, this is a realistic simulation where we find transition between the hyperbolicity regions: zone (I) while the material is flowing downhill, and it goes directly to zones (III) and (IV) when the granular mass is going into the water due to the large velocity difference (𝑢2≫ 𝑢1= 0). After that, the velocity difference decreases and it goes to zones (II) and (I). For this problem, treatments based on applying a friction produces unrealistic shape of the water in the entrance of granular material in the water, since they apply a large friction to move ||𝑢1−𝑢2||to the hyperbolic zone, as we will see later. We consider the domain [−22.5,7.5] and a uniform 𝛥𝑥 = 0.02. The bottom topography is 𝑏(𝑥) = max(𝑥+5,0), and the initial conditions (see Fig. 23) ℎ2(𝑥, 0) = {max (12.5 − 𝑏(𝑥),0)if 𝑥 > 5, 0otherwise; ℎ1(𝑥, 0) = max (5 − ℎ2(𝑥, 0) − 𝑏(𝑥),0), and 𝑞1(𝑥, 0) = 𝑞2(𝑥, 0) = 0. The friction coefficient is 𝜇𝑠= tan 30.5◦. Remark 4. Notice that the initial condition depicted in Fig. 23 does not correspond to a two-layer system. In that case, for the small times, we have to different fluids: a water body and an aerial avalanche. This means that the system should be adapted in order to describe this situation. For the ambient water, where the avalanche will fall into, no major modifications are needed. Indeed, system (1) degenerates to a shallow water system when the second layer vanishes and only the upper layer (with index 1) remains. Moreover, the scheme described in Section 3correctly handles such situations. For the case of the aerial avalanche, some simple modifications are needed, namely, 𝑟= 0 must Advances in Water Resources 182 (2023) 104587 22 M.J. Castro Díaz et al. Fig. 18. Test 6: Interfaces 𝜂1, 𝜂2(left-hand side) and velocities 𝑢1, 𝑢2(right-hand side) at times 𝑡= 0.5,1.0,1.5s in the case ||𝑢1−𝑢2||= 3 (non-hyperbolic initial condition), computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT), and with no correction (no treat.). Fig. 19. Test 6: (a) Interfaces 𝜂1, 𝜂2; (b) velocities 𝑢1, 𝑢2; (c,d) composite (𝐺)(4) and internal Froude numbers (𝐹 𝑟𝑖=||𝑢𝑖||∕(𝑔(1 − 𝑟)ℎ𝑖)) at initial time and 𝑡= 0.1s for the model considering treatments (NFT) and (CHMT) for the loss of hyperbolicity. Advances in Water Resources 182 (2023) 104587 23 M.J. Castro Díaz et al. Fig. 20. Test 6: Interfaces 𝜂1, 𝜂2(left-hand side) and velocities 𝑢1, 𝑢2(right-hand side) at times 𝑡= 0.3and 𝑡= 1.5s in the case for the model with treatment (CHMT) and the Shallow Water model with a fixed topography given by 𝜂2(𝑥, 0). Fig. 21. Test 7: Interfaces 𝜂1, 𝜂2(upper row), velocities 𝑢1, 𝑢2(middle row) and hyperbolicity zones (lower row) at times 𝑡= 0.05 (left column) and 𝑡= 0.5(right column) in the case ||𝑢1−𝑢2||= 6 (non-hyperbolic initial condition), computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT) and with no correction (no treat.). Advances in Water Resources 182 (2023) 104587 24 M.J. Castro Díaz et al. Fig. 22. Test 7: Interfaces 𝜂1, 𝜂2(upper row), velocities 𝑢1, 𝑢2(lower row) at times 𝑡= 0.5(left column) and 𝑡= 1.5(right column) in the case ||𝑢1−𝑢2||= 6 (non-hyperbolic initial condition), computed with the IFCP-FL scheme and treatment (CHMT) and the simplified model given by (34) in the whole domain (denoted by Simplified HM). Fig. 23. Test 8: Initial conditions for the bottom, water free surface and granular mass. be assumed when the material is out of the water (see Section 3.4 for more details). Fig. 24 shows the evolution of the water free surfaces and the granular mass for all the treatments, as well as the hyperbolicity zones for the treatments (NFT) and (CHMT). We do not show the hyperbolicity region for treatment (CFT) since it is always (I), as it is expected for this treatment. At the beginning, when the granular mass reaches the water surface, its velocity (𝑢2) is approximately 10 m∕s. This makes that the system is in zone (IV) during the times while the material is going into the water. In such situations, the system is hyperbolic and no treatment is needed. This is indeed the case when applying treatments (NFT) and (CHMT), whereas treatment (CFT) acts to reduce the velocity difference, increasing the magnitude of the water velocity (𝑢1). This makes that the water is strongly dragged, creating then a big wave (see for instance 𝑡= 1,1.5s in Fig. 24). Concerning the results with treatment (NFT), they are good till the moment that the flow goes into the complex region (zone (II)) from zone (III). At this moment, a friction is applied to move to the hyperbolic region and a shock is created at this location, similarly to what has been exposed in Test 7 (see Section 4.3.2). It can be observed looking at Fig. 24 (𝑡= 1 s). Once it happens, the results are qualitatively similar to what it is obtained when considering the classical friction treatment (CFT). However, quite different results are obtained when applying the treatment (CHMT), where the model is modified to its hyperbolic simplification just in those areas in zone (II). We see huge differences between the results obtained when this treatment is considered and those results with (CFT) and (NFT): (i) the generated wave in the water surface has a much smaller amplitude than the big wave generated when a friction treatment is considered; (ii) the granular material position is completely different, since the treatment based on adding friction has notably stopped it, in contrast with treatment (CHMT). We also see that spurious oscillations appear if no treatment is considered because of the non-hyperbolic nature of the system in some parts of the domain. It can be seen, for instance, for 𝑡= 2.5,3.45 s, and also in larger times over the surface of the granular material. In this figure, we also see that the granular material is almost steady with a nonflat surface at time 𝑡= 6.5, although it is not steady yet because of the water movement. It is only possible thanks to the correction of the numerical diffusion accounting for the friction force, as explained in Section 3.4. Actually, the friction makes that the oscillations appearing in the granular surface, when no treatment is applied, do not become smoother and remain steady. It is clear that a treatment is needed in this test, and (NFT) is not appropriate. In Fig. 24 we see that, when considering a friction treatment, a larger water wave is generated, and the water and granular mass are stuck at the front (where water and granular join) at early times, generating a shock in the water free surface. On the contrary, the results computed with (CHMT) are smoother, especially at the front. So, a natural question is treatment, (CFT) or (CHMT), is more appropriate for these applications. Although the system introduced here is a very simplified model for avalanches, taking a look at the laboratory experiments in Viroulet et al. (2014), in particular the entrance of the granular material into the water at early times, we see that the shape of the data at the front position is much closer to the one obtained when using the treatment (CHMT). On the contrary, the shock produced when considering the treatment (CFT) has no physical sense. We remark that we do not expect an accurate and exhaustive comparison with the laboratory results by using the simplified model presented here. To do so, one should consider an appropriate reference system (for instance local/curvilinear), non-hydrostatic effects and a more sophisticated model for the granular media (appropriate rheological law, dilatancy effects, friction between fluid/granular layers, ...), which is out of the scope of this paper. Nevertheless, a qualitative comparison tells us that, for these applications, the treatment (CHMT) is more appropriate than those based on applying a friction. As in previous test, we shall study the possibility of using the simplified hyperbolic model given by (34) in the whole domain. In Advances in Water Resources 182 (2023) 104587 25 M.J. Castro Díaz et al. Fig. 24. Test 8: Interfaces 𝜂1, 𝜂2and hyperbolic zone according to Fig. 7 at times 𝑡= 0,1, 1.5, 2.5, 4.5, 6.5 and 10 s computed with the IFCP-FL scheme and the different treatments for the loss of hyperbolicity: (CFT), (NFT), (CHMT) and with no correction (no treat.). Subfigures (b),(d) are zooms of subfigures (a),(c), respectively. this case, the free surface perturbations produced when considering (CHMT) are small, in contrast to the lower granular layer, except for the beginning when the material reaches the water surface. Actually, at this moment (granular material falling into the fluid) is when complex eigenvalues appear, and therefore in these cells (where larger free surface gradients appear) treatment (CHMT) considers the simplified model. So, in this situation both strategies discretize the same model, whereas these strategies differ where the free surface perturbations are small, which is the range of validity of the simplified hyperbolic model. For this reason, it is reasonable that we obtain close results in both cases for this test. In Fig. 25 we see the relative and absolute differences obtained for the heights ℎ1, ℎ2(lower row) and also, on the upper row, we show on the left-hand side the position of the granular material at time 𝑡= 4.7s (when bigger differences appear) and on the right-hand side the final time. In general we see that the relative errors on ℎ1are small, the bigger differences corresponding to ℎ2. Nevertheless, the relative errors are globally smaller than 7%, being approximately 4.5% when the steady state is reached. As conclusion, in this particular case, it may be reasonable to consider the simplified hyperbolic model in the whole domain, although this is not generally true in other configurations (as for instance in previous test or Test 9 in Section 4.4.2). Finally, let us discuss the computational effort required by each treatment for this test. The cheapest one, computationally speaking, corresponds to the friction treatment based on the classic approximation of the eigenvalues (CFT), which requires 1.7 s for this simulation. The treatments (CHMT) and (NFT) need 1.98 s and 2.14 s, respectively. We see that (CFT) is faster, since only a boundary of the non-hyperbolic region is computed, and Algorithm 1is not considered. Therefore, the approximated eigenvalues are not bounds of the exact eigenvalues, even in the hyperbolic region, as it has been shown in Fig. 8. It is not surprising also that treatment (NFT) requires more computational time than (CHMT) since, once the hyperbolic region is detected, (NFT) needs to solve a linear system to apply the friction correction whereas (CHMT)