ON MATHEMATICAL BALANCING OF A TWO-LAYER SHALLOW FLOW MODEL Wei-Koon Lee 1 , Alistair G.L. Borthwick 1 & Paul H. Taylor 1 1 Department of Engineering Science, University of Oxford, Oxford, U.K. ABSTRACT Two-layer shallow flow models offer an apparently cheap computational means of simulating density-stratified free surface flows. In the present study, we extend an algebraic method for balancing flux gradient and source terms proposed by Rogers et al. [Rogers B.D., Borthwick A.G.L., and Taylor P.H. 2003 Mathematical balancing of flux gradient and source terms prior to using Roe’s approximate Riemann solver, J. of Computational Physics, 192(2): 422- 451] to a two-layer shallow flow of distinct and immiscible liquids. It is found that the resulting scheme is well balanced for the steady state at rest but produces unwanted spurious oscillations in the presence of interface perturbation. We show that the imbalance is associated with the Bernoulli-type functions in the expanded form of the equations of motion after the elimination of the depth evolution terms. INTRODUCTION Owing to their relative simplicity in numerical implementation compared to the full Navier-Stokes equations, the nonlinear shallow water equations (SWEs) have become a popular choice for representing a range of predominantly horizontal free surface flows such as occur in shallow lakes, wide rivers, the coastal zones, estuaries, geophysical domains, oceans and the atmosphere. There exist a number of problems associated with near-horizontal flows that have a more or less pronounced two-layer structure, with each layer having a distinct density. Examples include fresh water intruding upon denser seawater, mudflow or debris flow beneath water, and circulation in a stratified ocean or lake. Physically, there is no exact separation interface between the two layers. Nonetheless, approximation using a two- dimensional dual-layer model can offer a considerable computational advantage in resolving the fluid mechanical behaviour of such flows for practical engineering purposes. In the last ten years, there has been rapid development of layered shallow flow models. Castro et al. (2001) constructed a first-order upwind scheme for one-dimensional flow of two superimposed layers of immiscible shallow liquids in a straight channel with constant rectangular cross-section. Assuming the existence and uniqueness of an entropy solution, the governing equations were linearised and the coupling terms included in the flux such that they may be treated as a standard system of conservation laws. The source term was treated using an approach proposed by Bermúdez & Vázquez-Cendón (1994), which satisfies the conservation property (C-property). Castro et al. (2006) implemented a layered model in two plan dimensions by considering a projected Riemann problem along the normal direction at every cell interface. The idea of the C-property in the context of single-layer SWEs was first introduced by Bermúdez & Vázquez-Cendón (1994), and requires the numerical scheme to calculate exactly the stationary solution corresponding to water at rest at grid nodes. For non-homogeneous hyperbolic systems of conservation laws, such as flows over realistic variable bathymetries, this is of particular importance. Otherwise, large errors in the numerical solution, in particular, the wave speed, may result in unphysical solutions (Roe, 1986). The problem is attributed to the conventional formulation of the SWEs where the surface gradient term within the momentum equations is split into an artificial flux gradient and a source term that includes the effect of the bed slope to ensure hyperbolicity. Greenberg & Leroux (1996) introduced an extension of the C-property to a more general condition by means of the concept of a ‘well-balanced’ scheme that preserves all equilibria of the system at grid nodes. Numerous approaches have since been reported in the literature and continuously improved for single layer models. In particular, the algebraic approach of Rogers et al. (2001 & 2003) uses a deviatoric form of the shallow water equations that inherently balances the set of hyperbolic equations without the need of additional computational effort. For the two-layer shallow flow model, there exists a major difficulty associated with the presence of non-conservative products, which are also common in many other multi-component or multi-phase models involving momentum and energy exchange. Furthermore, owing to nonlinear coupling between the layers, the system eigenvalues of a two-layer shallow water model are not readily accessible. In order to avoid the computational cost of solving the system eigenvalues, some authors have used the strategy of splitting which allows explicit access to the eigenstructure of the modified system, (e.g. Salmon, 2002; Bouchut & DeLuna, 2008; Abgrall & Karni, 2009; Chen & Peng, 2006; Chen et al., 2007). On the other hand, numerical solution of the system eigenvalues is unavoidable if a scheme is used which requires characteristic decomposition, such as the Roe-type method.