A consistent reduction of the two-layer shallow-water equations to an accurate one-layer spreading model

TL;DR

This paper mathematically derives an effective one-layer shallow-water model from the two-layer equations that accurately captures shock behavior without extra closure relations, simplifying analysis of fluid spreading phenomena.

Contribution

It provides a rigorous derivation of a one-layer model that correctly describes shocks in two-layer flows, challenging previous assumptions about the model's limitations near shocks.

Findings

The derived model accurately predicts shock speeds consistent with empirical and theoretical models.

Numerical simulations show rapid convergence of the one-layer model to two-layer solutions.

Predictions from the model agree well with dam break experiments.

Abstract

The gravity-driven spreading of one fluid in contact with another fluid is of key importance to a range of topics. To describe these phenomena, the two-layer shallow-water equations is commonly employed. When one layer is significantly deeper than the other, it is common to approximate the system with the much simpler one-layer shallow water equations. So far, it has been assumed that this approximation is invalid near shocks, and one has applied additional front conditions for the shock speed. In this paper, we prove mathematically that an effective one-layer model can be derived from the two-layer equations that correctly captures the behaviour of shocks and contact discontinuities without any additional closure relations. The proof yields a novel formulation of an effective one-layer shallow water model. The result shows that simplification to an effective one-layer model is well justified mathematically and can be made without additional knowledge of the shock behaviour. The shock speed in the proposed model is consistent with empirical models and identical to the front conditions that have been found theoretically by e.g. von K\'arm\'an and by Benjamin. This suggests that the breakdown of the shallow-water equations in the vicinity of shocks is less severe than previously thought. We further investigate the applicability of the shallow water framework to shocks by studying shocks in one-dimensional lock-exchange/lock-release. We derive expressions for the Froude number that are in good agreement with the widely employed expression by Benjamin. We then solve the equations numerically to illustrate how quickly the proposed model converges to solutions of the full two-layer shallow-water equations. We also compare numerical results using our model with results from dam break experiments. Predictions from the one-layer model are found to be in good agreement with experiments.

Figures

Tables

Table 1. Table 1 : Initial conditions used by the one-dimensional dam-break experiments. The experiments by Suchon are with oil spreading on water, while those from Chang, Reid, and Fay are liquified methane and liquified nitrogen spreading on water.

Authors	Experiment	Spill Volume ($L$)	Height $h_0$ ($\mathrm{cm}$)	Width $L$ ($\mathrm{cm}$)	$\delta$
Suchon	Run 11	$10$	$16.51$	$10.16$	$0.1$
Suchon	Run 14	$7.7$	$16.637$	$7.62$	$0.1$
Suchon	Run 17	$5.1$	$10.9982$	$7.62$	$0.1$
Suchon	Run 18	$5.1$	$16.51$	$5.08$	$0.1$
Chang, Reid, and Fay	$2L$ methane	$2.0$	$17.3$	$7$	$0.746$
Chang, Reid, and Fay	$2L$ nitrogen	$2.0$	$17.3$	$7$	$0.34$
Chang, Reid, and Fay	$0.75L$ methane	$0.75$	$6.5$	$7$	$0.746$
Chang, Reid, and Fay	$1L$ nitrogen	$1.0$	$8.7$	$7$	$0.34$

Equations

$$\mathrm{Fr}=\frac{u}{\sqrt{g^{\prime }h}},$$

$${\mathrm{Fr}}_{\text{LE}}^2=\frac{(1-\alpha )(2-\alpha )}{(1+\alpha )},$$

$$u_{\text{LE}}=\beta \sqrt{\frac{g^{\prime }V}{A}},$$

$$\frac{\partial }{\partial t}\rho h+\nabla \cdot (\rho h\mathbf{u})=G_h,$$

$$\begin{array}{rl}\frac{\partial }{\partial t}(\rho h\mathbf{u})& +\nabla \cdot (\rho h\mathbf{u}\otimes \mathbf{u})\\ & +\nabla \left(\frac{1}{2}g\rho h^2\right)={\mathbf{G}}_{hu}-g\rho h\nabla b.\end{array}$$

$$\frac{\partial }{\partial t}\rho h+\nabla \cdot (\rho h\mathbf{u})=G_h,$$

$$\frac{\partial }{\partial t}\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}+g\nabla (h+b)=\frac{1}{\rho h}({\mathbf{G}}_{hu}-\mathbf{u}{G}_h).$$

$$\frac{\partial }{\partial t}{\rho }_1{h}_1+\nabla \cdot ({\rho }_1{h}_1{\mathbf{u}}_1)$$

$$\frac{\partial }{\partial t}{\rho }_2{h}_2+\nabla \cdot ({\rho }_2{h}_2{\mathbf{u}}_2)$$

$$\frac{\partial }{\partial t}({\rho }_1{h}_1{\mathbf{u}}_1)+\nabla \cdot ({\rho }_1{h}_1{\mathbf{u}}_1\otimes {\mathbf{u}}_1)+\nabla \left(\frac{1}{2}g{\rho }_1{h}_1^2\right)$$

$$\frac{\partial }{\partial t}({\rho }_2{h}_2{\mathbf{u}}_2)+\nabla \cdot ({\rho }_2{h}_2{\mathbf{u}}_2\otimes {\mathbf{u}}_2)+\nabla \left(g{\rho }_1{h}_1{h}_2+\frac{1}{2}g{\rho }_2{h}_2^2\right)$$

$$\nabla \cdot \left({\rho }_i{h}_i{\mathbf{u}}_i\otimes {\mathbf{u}}_i\right)={\mathbf{u}}_i\nabla \cdot ({\rho }_i{h}_i{\mathbf{u}}_i)+{\rho }_i{h}_i\left({\mathbf{u}}_i\cdot \nabla \right){\mathbf{u}}_i,$$

$$\frac{\partial }{\partial t}{\rho }_1{h}_1+\nabla \cdot ({\rho }_1{h}_1{\mathbf{u}}_1)$$

$$\frac{\partial }{\partial t}{\rho }_2{h}_2+\nabla \cdot ({\rho }_2{h}_2{\mathbf{u}}_2)$$

$$\frac{\partial }{\partial t}{\mathbf{u}}_1+({\mathbf{u}}_1\cdot \nabla ){\mathbf{u}}_1+\nabla \left[g(h_1+h_2+b)\right]$$

$$\frac{\partial }{\partial t}{\mathbf{u}}_2+({\mathbf{u}}_2\cdot \nabla ){\mathbf{u}}_2+\nabla \left[g\left(\frac{{\rho }_1}{{\rho }_2}{h}_1+h_2+b\right)\right]$$

$$\frac{\partial }{\partial t}{\rho }_1{h}_1+\nabla \cdot ({\rho }_1{h}_1{\mathbf{u}}_1)=G_{h_1},$$

$$\frac{\partial }{\partial t}{\rho }_2{h}_2+\nabla \cdot ({\rho }_2{h}_2{\mathbf{u}}_2)=G_{h_2},$$

$$\begin{array}{rl}\frac{\partial }{\partial t}\left({\rho }_1{h}_1{\mathbf{u}}_1+{\rho }_2{h}_2{\mathbf{u}}_2\right)& +\nabla \cdot \left({\rho }_1{h}_1{\mathbf{u}}_1\otimes {\mathbf{u}}_1+{\rho }_2{h}_2{\mathbf{u}}_2\otimes {\mathbf{u}}_2\right)\\ & +\nabla \left(\frac{1}{2}g{\rho }_1{h}_1^2+{\rho }_{1}g{h}_1{h}_2+\frac{1}{2}{\rho }_{2}g{h}_2^2\right)={\mathbf{G}}_{h_1{u}_1}+{\mathbf{G}}_{h_2{u}_2}-g({\rho }_1{h}_1+{\rho }_2{h}_2)\nabla b,\end{array}$$

$$\frac{\partial }{\partial t}\left({\mathbf{u}}_2-{\mathbf{u}}_1\right)+({\mathbf{u}}_2\cdot \nabla ){\mathbf{u}}_2-({\mathbf{u}}_1\cdot \nabla ){\mathbf{u}}_1-\nabla \left(g\delta h_1\right)=\mathbf{J},$$

$$\mathbf{J}=\frac{{\mathbf{G}}_{h_2{u}_2}-{\mathbf{u}}_2{G}_{h_2}}{{\rho }_2{h}_2}-\frac{{\mathbf{G}}_{h_1{u}_1}-{\mathbf{u}}_1{G}_{h_1}}{{\rho }_1{h}_1}$$

$$\delta :=\frac{{\rho }_2-{\rho }_1}{{\rho }_2}.$$

$$\frac{\partial }{\partial t}u(t,\mathbf{x})+\nabla \cdot \mathbf{q}(t,\mathbf{x})=J$$

$$S\left\llbracket{u}\right\rrbracket=\boldsymbol{\hat{n}}\cdot\left\llbracket{\boldsymbol{q}}\right\rrbracket.$$

$$\frac{\partial }{\partial t}\mathbf{u}(\mathbf{x},t)+\nabla \cdot (\mathbf{a}\otimes \mathbf{b})+\nabla q(\mathbf{x},t)=\mathbf{J},$$

$$S\left\llbracket{\boldsymbol{u}}\right\rrbracket=\left\llbracket{\boldsymbol{\hat{n}}\cdot(\boldsymbol{a}\otimes\boldsymbol{b})+q\boldsymbol{\hat{n}}}\right\rrbracket.$$

$$S\left\llbracket{\boldsymbol{u}_{i}}\right\rrbracket\cdot\boldsymbol{\hat{n}}=\left\llbracket{\frac{1}{2}(\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{i})^{2}+g\left(\frac{\rho_{1}}{\rho_{2}}\right)^{i-1}h_{1}+gh_{2}}\right\rrbracket,$$

$$S\left\llbracket{\rho_{1}h_{1}\boldsymbol{u}_{1}+\rho_{2}h_{2}\boldsymbol{u}_{2}}\right\rrbracket\\ =\left\llbracket{(\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{1})\rho_{1}h_{1}\boldsymbol{u}_{1}+(\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{2})\rho_{2}h_{2}\boldsymbol{u}_{2}}\right\rrbracket\\ +\left\llbracket{\frac{1}{2}g\rho_{1}h_{1}^{2}+\rho_{1}gh_{1}h_{2}+\frac{1}{2}\rho_{2}gh_{2}^{2}}\right\rrbracket\boldsymbol{\hat{n}}$$

$$S\boldsymbol{\hat{n}}\cdot\left\llbracket{\boldsymbol{u}_{2}-\boldsymbol{u}_{1}}\right\rrbracket=\Biggl{\llbracket}\frac{1}{2}\biggl{[}(\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{2})^{2}\\ -(\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{1})^{2}\biggr{]}-g\delta h_{1}\Biggr{\rrbracket}.$$

$$E=\frac{1}{2}\left(\rho_{1}h_{1}\left\lvert{\boldsymbol{u}_{1}}\right\rvert^{2}+\rho_{2}h_{2}\left\lvert{\boldsymbol{u}_{2}}\right\rvert^{2}\right)\\ +g\Biggl{[}\rho_{2}h_{2}\left(\frac{1}{2}h_{2}+b\right)+\left(\frac{1}{2}h_{1}+h_{2}+b\right)\rho_{1}h_{1}\Biggr{]}.$$

$$\begin{gathered} \frac{\partial E}{\partial t}+\nabla \cdot \left[{\mathbf{q}}_1\left(g\left(h_1+h_2+b\right)+\frac{1}{2}{\left|{\mathbf{u}}_1\right|}^2\right)+{\mathbf{q}}_2\left(g\left(\frac{{\rho }_1}{{\rho }_2}{h}_1+h_2+b\right)+\frac{1}{2}{\left|{\mathbf{u}}_2\right|}^2\right)\right] \\ \le {\mathbf{u}}_1\cdot {\mathbf{G}}_{h_1{u}_1}+{\mathbf{u}}_2\cdot {\mathbf{G}}_{h_2{u}_2}-\frac{1}{2}g{h}_1^2\frac{\partial {\rho }_1}{\partial t}-g{h}_2\left(\frac{{\rho }_1}{{\rho }_2}{h}_1+\frac{1}{2}{h}_2\right)\frac{\partial {\rho }_2}{\partial t} \\ +G_{h_1}\left(g\left(h_1+h_2+b\right)-\frac{1}{2}{\left|{\mathbf{u}}_1\right|}^2\right)+G_{h_2}\left(g\left(\frac{{\rho }_1}{{\rho }_2}{h}_1+h_2+b\right)-\frac{1}{2}{\left|{\mathbf{u}}_2\right|}^2\right), \end{gathered}$$

$$\frac{\partial h_1}{\partial t}+\frac{\partial }{\partial x}(h_1{u}_1)=0,$$

$$\frac{\partial h_2}{\partial t}+\frac{\partial }{\partial x}(h_2{u}_2)=0,$$

$$\frac{\partial u_1}{\partial t}+\frac{\partial }{\partial x}\left(\frac{1}{2}{u}_1^2+g{h}_2+\frac{1}{{\rho }_1}{p}_0\right)=0,$$

$$\frac{\partial u_2}{\partial t}+\frac{\partial }{\partial x}\left(\frac{1}{2}{u}_2^2+g{h}_2+\frac{1}{{\rho }_2}{p}_0\right)=0.$$

$$S=u_{2,l}=-\frac{h_{1,l}}{h_{2,l}}{u}_{1,l},$$

$$\begin{gathered} S\left({\rho }_{2}S+{\rho }_1\frac{h_{2,l}}{h_{1,l}}S\right)=\frac{1}{2}\left({\rho }_2{S}^2-{\rho }_1\frac{h_{2,l}^2}{h_{1,l}^2}{S}^2\right) \\ +{\rho }_{2}g{h}_{2,l}-{\rho }_{1}g{h}_{2,l}. \end{gathered}$$

Full text

I Introduction

The spreading of two layers of fluids with different density is of considerable importance. It has been an active field of study since at least 1774, when Franklin, Brownrigg, and Farish (1774) investigated how oil spreads on water and how this can be used to still waves. Applications where this phenomenon plays an important role include spills of oil (Hoult, 1972; Fay, 1971; Chebbi, 2001) and liquefied gaseous fuels,(Fay, 2007, 2003; Brandeis and Ermak, 1983) stratified flow inside pipes,(Stanislav, Kokal, and Nich, 1986) gravity currents particularly in geophysical systems,(Adduce, Sciortino, and Proietti, 2012; Shin, Dalziel, and Linden, 2004; Moodie, 2002; Mériaux et al., 2016) monomolecular layers for evaporation control,(Stickland, 1972) and coalescence in three-phase fluid systems.(Mar and Mason, 1968) These applications include non-miscible fluids such as oil and water, or systems with miscible fluids at large Richardson number, i.e. where buoyancy dominates mixing effects and ensures separation into layers.

A fundamental property of spreading phenomena is the rate of spreading, or the speed of the leading edge of the spreading fluid. This is typically characterized by the dimensionless Froude number,(White, 2011; Vaughan and O’Malley, 2005)

[TABLE]

where $u$ is the velocity, $h$ is the height of the layer that is spreading, and $g^{\prime}$ is the effective gravitational acceleration. In two layer spreading, the effective gravitational acceleration is $g^{\prime}=(1-\rho_{1}/\rho_{2})g$, where $\rho_{1}$ and $\rho_{2}$ are the two fluid densities and $\rho_{1}<\rho_{2}$.

An early result for the Froude number of gravity currents was presented by von Kármán.von Kármán (1940) They found that for the edge of a spreading gravity current at semi-infinite depth, $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$, where the subscript is short for “leading edge”. Benjamin (1968) later developed a model for $\mathrm{Fr}_{\text{\tiny LE}}$ for spreading of gravity currents with constant height,

[TABLE]

where $\alpha=h_{2}/(h_{1}+h_{2})$. Here $h_{1}$ and $h_{2}$ are the heights of the top and bottom layers, respectively. This model approaches the result by von Kármán when the bottom layer becomes thin, $h_{2}\ll h_{1}$. More recently, Ungarish (2017) extended the result of Benjamin to the spreading of gravity currents into a lighter fluid with an open surface. This result also gives $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$ when the spreading fluid becoms relatively much thinner than the ambient fluid.

The next step beyond characterizing spreading rates is to develop a model that predicts the phenomenon in more detail. An early model was presented by Fay,Fay (1969) who studied the spreading of oil on water. They divided the spreading into three phases; one where inertial forces dominate, one where viscous forces dominate, and one where the surface tension dominates. In the inertial phase, the speed of the front can be written as

[TABLE]

where $\beta$ is an empirical constant and $V$ and $A$ are the volume and area, respectively. Then $V/A$ is the average height of the spreading oil. In this model, $\beta$ represents an effective Froude number where the height at the leading edge is approximated by the average height. The value of $\beta$ has been discussed in the literature and is commonly set to $\beta=1.31$ in the one-dimensional case and $\beta=1.41$ in the axisymmetric case.(Fannelop and Waldman, 1972; Fay, 1971; Hoult, 1972; Fay, 2007)

A more general approach than the Fay model is the two-layer shallow-water equations (2LSWE), which are derived from the Euler equations by assuming a negligible vertical velocity.(Ovsyannikov, 1979; Vreugdenhil, 1979) These equations model the flow of two layers of shallow liquids and may be used to simulate for instance gravity currents.(Audusse et al., 2011) However, internal breaking of waves or large differences in velocities of the two layers can break the hyperbolicity of the equations. Even if the initial conditions are hyperbolic, the system can evolve into a non-hyperbolic state.(Milewski et al., 2004) A breakdown of hyperbolicity causes problems such as ill-posedness and Kelvin-Helmholtz like instabilities.(Lannes and Ming, 2015; Stewart and Dellar, 2013; Lam, Ghidaoui, and Kolyshkin, 2016) Non-hyperbolic equations are generally more difficult to analyse and computationally much more expensive to solve than hyperbolic equations.(Bouchut and Morales de Luna, 2008) Attempts to amend the non-hyperbolicity of the systems include adding numerical (non-physical) friction forces,(Castro-Díaz et al., 2011) operator-splitting approaches,(Bouchut and Zeitlin, 2010) and introduction of an artificial compressibility.(Chiapolino and Saurel, 2018)

Due to their comparative simplicity, the one-layer shallow-water equations (1LSWE) have often been used to model two-layer phenomena like liquid-on-liquid spreading and gravity currents where one assumes that the layers are in a buoyant equilibrium. In this case, a forced constant Froude-number boundary condition at the leading edge of a spreading liquid is used to account for the effect of the missing layer. (Fannelop and Waldman, 1972; Hoult, 1972; Hatcher and Vasconcelos, 2014) The additional boundary condition at the leading edge has also been used in combination with the 2LSWE.(Rottman and Simpson, 1983; Ungarish, 2013) In particular, Rottman and Simpson (1983) argued that a front condition that includes the Froude number is necessary because viscous dissipation and vertical acceleration are too significant to be neglected at the front.

The 1LSWE are always hyperbolic and therefore have fewer challenges than the 2LSWE. However, there are situations where even the 1LSWE are not strictly hyperbolic, meaning that the two eigenvalues of the Jacobian coincide. This situation is found when considering the wet–dry transition, such as the dam break on a dry bottom, or for certain bottom topographies. In particular the case of a gravity current flowing upslope, as in a shallow water wave encountering a beach, is of importance and has seen new developments in recent years.(Lombardi et al., 2015; Bjørnestad and Kalisch, 2017; Zemach et al., 2019) There is an exhaustive literature on the subject of hyperbolicity of the 1LSWE, (Fraccarollo and Toro, 1995; Zhou et al., 2001; LeFloch and Thanh, 2007; Liang and Marche, 2009; LeFloch and Thanh, 2011; Murillo and Navas-Montilla, 2016) including the topic of well-balanced formulation, the more general E-balanced schemes, and the identification of resonant versus non-resonant regimes of flow. These points are mainly of interest for the numerical solution of the equations in specific regimes. As the present paper is focused more on the theoretical developments, a detailed discussion of hyperbolicity is beyond the scope of the present work.

The main results of the present paper are the following. First, we show that the need to impose boundary conditions or empirical closures for the spreading rate when using the 1LSWE instead of the 2LSWE follows from the different shock behaviour of the two formulations. Second, we demonstrate that weak solutions of the 2LSWE converge to weak solutions of a locally conservative form of the one-layer equations. This formulation is different from the standard 1LSWE, and removes the need for front conditions.

This is a strong result as it implies that in many situations, such as when considering liquid spills on water or ocean layers in deep water, one may use the much simpler locally conservative 1LSWE even for two-layer spreading phenomena, without the need for additional boundary conditions or closures. An example is presented in figure 1, which illustrates how solutions to different forms of the 1LSWE compare to the solution of the 2LSWE for a dam-break problem. The figure shows a clear difference between the locally and globally conservative 1LSWE.

We further demonstrate that the constant Froude number at the front of an expanding fluid can be derived directly from the 2LSWE. The Froude numbers obtained from the analysis in this paper are in excellent agreement previous results from the literature. This indicates that the breakdown of the shallow-water equations in vicinity of shocks is less severe than previously suggested.

The paper is structured as following. In section II, we introduce the two-layer shallow-water equations (2LSWE), the one-layer shallow-water equations (1LSWE) and the Rankine-Hugoniot condition for the shock. In section III we derive expressions for the Froude number from the full two-layer shallow-water equations. The key result of the paper is presented in section IV, where we show the 2LSWE can be approximated by a one-layer model when the upper layer is much thicker than the bottom layer, as well as in the opposite situation. In section V we define some numerical experiments that are used in section VI to study how solutions of the 2LSWE approach the one-layer approximations. We show that the results from the simplified model are in good agreement with experimental results. Concluding remarks are provided in section VII.

II Theory of the shallow-water equations

Consider a two-layer system where a fluid of lower density spreads on top of another fluid, as illustrated in figure 2. Assuming that the layers are shallow, the solution of the two-layer shallow-water equations (2LSWE) gives the evolution of height and horizontal velocity of both fluids as a function of position and time.

In the following, we first describe the well-known one-layer shallow-water equations (1LSWE). A straightforward generalization to the 2LSWE is presented next, where we discuss two approaches for reformulating the 2LSWE in a manner that makes them suitable for reduction to an effective one-layer model. We then show how the Rankine-Hugoniot conditions can be used to predict the shock speed. Subsequently we employ the vanishing-viscosity regularization and travelling wave solutions to obtain physically acceptable solutions of the partial differential equations (PDEs). At the end of the section, we present a necessary energy requirement for the 2LSWE that is used to select correct physical solutions in section IV.

II.1 The one-layer shallow-water equations

The 1LSWE are typically presented in a globally conservative form where total momentum is conserved,(LeVeque, 2002)

[TABLE]

where $\rho$ is the density, $h$ is the height, $\boldsymbol{u}$ is the vertically averaged horizontal velocity, $\otimes$ denotes the tensor product, $b$ is the bottom topography, $G_{h}$ and $\boldsymbol{G}_{hu}$ are source functions that may represent external phenomena, such as evaporation, Coriolis forces, wind shear stress, or interfacial shear forces. The bottom topography is assumed to be continuous throughout. The density $\rho$ is assumed constant in space, although it may vary in time.

One may also consider what will be referred to as the locally conservative 1LSWE, that is

[TABLE]

Here the continuity equation (5a) is unchanged. The various forms of the one-layer and two-layer equations all use the same form of the continuity equation.

One particularly striking difference between eq. 5 and eq. 4 is the admissibility of shocks when the height drops to 0. This will be further discussed in section V.1, but the upshot is that such a shock is impossible in eq. 4, while in eq. 5 it is possible with a Froude number $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$. This is exactly the result by von Kármán (1940) for two layer spreading with such shocks. In fact, in section IV, we show that the locally conservative form correctly captures the two-layer behaviour in certain limits. This result is consistent with previous results which show that numerical approaches will fail to solve the conservation of global momentum.(Bouchut and Zeitlin, 2010)

II.2 The two-layer shallow-water equations

The 2LSWE may be written in a general, layerwise form with arbitrary source terms as

[TABLE]

where the subscripts 1 and 2 denote the top and bottom layers respectively. The coupling between the two layers are captured by the last source terms on the right-hand side of the momentum equations.

This form was originally described by Ovsyannikov,Ovsyannikov (1979) and is referred to in more recent works as “the conventional two-layer shallow-water model”.(Chiapolino and Saurel, 2018)

II.3 2LSWE forms that are reducible to one-layer approximations

Conservation of momentum can be considered at three different scales:

The globally conservative form where total momentum is conserved. 2. 2.

The layerwise conservative form (eq. 6) where the momentum in each layer is conserved. 3. 3.

The locally conservative form where the local momentum, or velocity, is conserved.

Although these formulations are equivalent for smooth solutions, they are not generally equivalent, as will be further discussed in section II.4. The layerwise formulation is not easily reducible to a one-layer model. The remaining two approaches can be converted to an effective one-layer approximation, and our analysis will cover both. In the locally conservative form, we combine eqs. 6c and 6d with eqs. 6a and 6b to give equations for velocity rather than momentum. Using the product rule for differentiation,

[TABLE]

we arrive at the set of equations which we shall refer to as the locally conservative version of the 2LSWE,

[TABLE]

For a comprehensive study of the well-posedness of the locally conservative 2LSWE, see for instance.(Monjarret, 2015)

When conserving the total momentum, the sum of eq. 6c and eq. 6d is used, which has the advantage of eliminating the interaction between the layers. However, this approach requires an additional conservation law. Ostapenko (1999, 2001) showed that the additional conservation law should be the difference between eq. 7d and eq. 7c. Ostapenko used these equations in a study of the well-posedness of the 2LSWE. The resulting equations, which we will refer to as the globally conservative version of the 2LSWE, read

[TABLE]

where

[TABLE]

and where we have defined

[TABLE]

II.4 The Rankine-Hugoniot condition

When two sets of equations are equivalent in the classical sense, they may not be equivalent in the weak sense, that is, when interpreted as distributions.(Whitham, 1974; Holden and Risebro, 2015; Borthwick, 2016) In the 2LSWE, eqs. 6, 7 and 8 are equivalent for smooth solutions, but not for weak solutions. In particular, these equations will give different shock velocities. We shall next discuss the mathematical framework used to analyze such discontinuities; the Rankine-Hugoniot condition, named after Rankine and Hugoniot who first introduced it.(Rankine, 1870; Hugoniot, 1887, 1889)

The Rankine-Hugoniot condition states the following. Assume that $u$ satisfies a general scalar conservation equation

[TABLE]

in the weak sense, where $J$ is some source term that does not involve the derivatives of $u$. Further, assume that $u$ has a discontinuity along some curve $\Gamma$. For any function $f$, define the jump across a discontinuity as \left\llbracket{f}\right\rrbracket\equiv f_{r}-f_{l}, where $f_{r}\equiv\lim_{\varepsilon\to 0^{+}}f(\boldsymbol{\xi}+\varepsilon\boldsymbol{\hat{n}})$ and $f_{l}\equiv\lim_{\varepsilon\to 0^{-}}f(\boldsymbol{\xi}+\varepsilon\boldsymbol{\hat{n}})$. The Rankine-Hugoniot condition then states that the discontinuity at any point $\boldsymbol{\xi}\in\Gamma$ propagates along the outward-pointing normal vector $\boldsymbol{\hat{n}}$ with a speed $S$. This speed is called the shock speed and satisfies the relation

[TABLE]

Similarly, if $\boldsymbol{u}$ satisfies a general vector conservation equation,

[TABLE]

then, if there is some discontinuity in $\boldsymbol{u}$, we have the result

[TABLE]

Equations 11 and 13 can be directly applied to the mass conservation equations and the conservation law for total momentum, respectively. In one dimension, the Rankine-Hugoniot conditions can also be applied to the locally conservative momentum equation. In two dimensions, the term $\boldsymbol{u}\cdot\nabla\boldsymbol{u}$ renders the Rankine-Hugoniot condition for the transversal velocity component ill-defined. Nevertheless, for our purposes we do not need the Rankine-Hugoniot condition for the transversal velocity component. See appendix A for a discussion on this. In the layerwise momentum equation, the interaction term $\propto h_{1}\nabla h_{2}$ makes the normal component for the momentum equations ill-defined, which is why we must exclude this formulation of the 2LSWE from the analysis.

We derive the Rankine-Hugoniot conditions in appendix A and find that for the locally conservative 2LSWE (eqs. 7c and 7d),

[TABLE]

where as before $i=1,2$ denotes the layer. Similarly, for the globally conservative 2LSWE eqs. 8c and 8d, we find

[TABLE]

and

[TABLE]

Finally, we note that in calculations with the Rankine-Hugoniot condition it is useful to observe that \left\llbracket{ab}\right\rrbracket=\left\llbracket{a}\right\rrbracket\left\langle{b}\right\rangle+\left\langle{a}\right\rangle\left\llbracket{b}\right\rrbracket where $\left\langle{a}\right\rangle=(a_{l}+a_{r})/2$.

II.5 Physical solutions

When PDEs are considered in the weak sense, it is necessary to impose extra conditions to extract a unique physical solution. Such conditions are called entropy conditions. In this subsection, we will introduce one such condition: the energy requirement. For simplicity, we define a physical solution as one that satisfies the energy requirement.

The energy requirement states that only shocks that dissipate energy are physical. This translates into requiring that the energy of the physical solution does not increase in time except from possible source terms. Energy, in this sense, has the role of a mathematical entropy.Holden and Risebro (2015) However, the word entropy is typically restricted to convex functions of the solution variables. As has been showed by Ostapenko,Ostapenko (1999) energy is indeed a convex function of the globally conservative system, but for the locally conservative system it is convex only for subcritical flow. Because we here cover both cases we use the word energy rather than entropy.

The energy of the 2LSWE reads

[TABLE]

This expression is given in terms of parameters that are already solved for in the 2LSWE. For smooth solutions we may therefore combine the subequations of the 2LSWE to form a conservation law for the energy. By exchanging the equality in this conservation law by an inequality, it can be fulfilled also by weak, discontinuous solutions. We obtain

[TABLE]

where $\boldsymbol{q}_{i}=\rho_{i}h_{i}\boldsymbol{u}_{i}$ for short.

III Derivation of Froude numbers from the 2LSWE

In the following, we briefly illustrate the surprising effectiveness of the 2LSWE to predict shock speeds despite its underlying assumption of negligible vertical acceleration. To do this we apply the Rankine-Hugoniot conditions and the 2LSWE to derive expressions for the leading edge Froude number ($\mathrm{Fr}_{\text{\tiny LE}}$) of two-layer systems with fixed total height.

Shock speeds in two-layer systems with fixed total height is important for instance in lock-exchange and lock-release problems, where a heavy fluid is spreading within a lighter fluid inside a rectangular channel as illustrated in figure 3. Such problems have been studied extensively, and there is a large number of results from laboratory experiments available.(Rottman and Simpson, 1983; Huppert and Simpson, 1980; Shin, Dalziel, and Linden, 2004) Moreover, much theoretical work has been carried out to model the Froude-number for flows inside rectangular channels,(Benjamin, 1968; Priede, 2019; Borden and Meiburg, 2013) which means that this is a good candidate for testing the credibility of shock behaviour in the 2LSWE.

Most previous works have focused on fluids with similar densities such that $\delta\ll 1$ for $\delta$ given by eq. 9. This is referred to as the Boussinesq case.(Boussinesq, 1903) The most commonly used front condition applied to such flows is the equation for the Froude number given by Benjamin,Benjamin (1968) eq. 2. For instance, Ungarish (2011) has applied the Froude number by Benjamin as a boundary condition when solving the 2LSWE for rectangular geometries. They also generalized this to arbitrary geometries.(Ungarish, 2013)

For the particular problem where the two-layer flow is confined inside a rectangular channel, the sum of the layer depths must be constant; $h_{1}+h_{2}=H$. In this case, there are no free surfaces. We therefore add an additional pressure term $p_{0}$ that may vary in time and space but is constant in the vertical direction.

We first consider the locally conservative 2LSWE (7) in one spatial dimension with the added free pressure term,

[TABLE]

These are the same equations that were used by Rottman and Simpson (1983) to study spreading of gravity currents. Rottman and Simpson added eq. 2 as an additional equation for the Froude number at the leading edge, but in the following we will show that a similar expression for $\mathrm{Fr}_{\text{\tiny LE}}$ can be obtained from eq. 19 directly.

With $h_{1}+h_{2}$ constant, the sum of eqs. 19a and 19b implies that $h_{1}u_{1}+h_{2}u_{2}$ is constant in $x$. If we assume that the total momentum is [math] at the boundary, e.g. due to a wall or because the boundary is at infinity and the fluids were initially at rest, we may set $h_{1}u_{1}+h_{2}u_{2}=0$.

By use of the Rankine-Hugoniot condition (11) to eqs. 19a and 19b, we get

[TABLE]

where, as before, the subscript $l$ indicates the left side of the shock. We next apply the Rankine-Hugoniot condition to $\rho_{2}$(19d) $-$ $\rho_{1}$(19c), which gives

[TABLE]

After some algebraic manipulation, we find that

[TABLE]

where $\alpha=h_{2,l}/(h_{2,l}+h_{1,l})$.

We next consider the globally conservative 2LSWE (8). A similar analysis and derivation now gives

[TABLE]

As expected, the different formulation of the 2LSWE leads to a different expressions for the Froude number.

The Boussinesq approximation is achieved by setting $\delta=0$ wherever it is not multiplied by $g$. In this case eqs. 21 and 22 coincide and gives that

[TABLE]

Figure 4 compares our results from the 2LSWE, eq. 23, to the model by Benjamin,Benjamin (1968) eq. 2. As can be seen, the difference is small. Equation 2 is obtained by balancing forces and does not rely on any assumptions regarding negligible vertical velocities. The similarity of eqs. 2 and 23 therefore indicates that the breakdown of the shallow-water equations in vicinity of shocks is not so severe as one would think and as has been repeatedly assumed in the literature.(Hatcher and Vasconcelos, 2014; Rottman and Simpson, 1983; Hoult, 1972; Fannelop and Waldman, 1972)

One advantage of the 2LSWE is that it does not use the Boussinesq approximation. The non-Boussinesq case has more recently received attention in the literature,Lowe, Rottman, and Linden (2005); Birman, Martin, and Meiburg (2005) and eqs. 21 and 22 could be of interest in this regard.

The treatment presented here is under the assumption of negligible mixing between the layers. In systems with mixing, Sher and Woods (2015) has found the spreading is slower because the density difference at the leading edge, and hence the effective gravity, is reduced with time. With their time-dependent reduced gravity, they found experimentally that $\mathrm{Fr}_{\text{\tiny LE}}=0.90\pm 0.05$ for $\alpha=0.37$. Inserting $\alpha=0.37$ into eq. 23 we get $\mathrm{Fr}_{\text{\tiny LE}}=0.89$. That is, if mixing is taken into account in the shallow water framework by introducing a slowly varying time-dependent density difference and possibly some source terms that do not affect the Rankine-Hugoniot condition, the resulting Froude number at the leading edge is in good agreement with the observed value.

Finally, we note that Priede (2019) has also found an expression for the Froude number in the 2LSWE with constant height. They restricted the analysis to the Boussinesq case and got a result which differs slightly from eq. 23. The reason for the deviation is that they rewrote the equations in terms of new variables, $\eta\vcentcolon=h_{1}-h_{2}$ and $\vartheta\vcentcolon=u_{1}-u_{2}$, and used $\eta$ and $\eta\vartheta$ as conserved quantities before they applied the Rankine-Hugoniot condition. This changes the weak solutions and hence the shock speed.

IV Reducing the two-layer systems to effective one-layer systems

In this section, we present a theorem with a constructive proof that demonstrates that it is possible to reduce the 2LSWE into an effective one-layer model while preserving the correct behaviour of shocks and contact discontinuities. The theorem shows that this decoupling is possible when the depth of one layer becomes large compared to the other layer. We show that additional closures for the shock velocity are not needed, which differs from previous reductions to one-layer models presented in the literature.

IV.1 The constant-height lemma

In the following, we denote by $s$ and $d$ the relatively shallow and deep layers, respectively. This means that with $(s,d)=(1,2)$, the top layer is shallow relative to the bottom layer, and vice versa for $(s,d)=(2,1)$. Further, we let $\overline{f}$ denote the average of $f$ over the region in which it is defined.

In order to state and prove the theorem, we will use a concept we call source-boundedness. We will also use a lemma that states that in the indicated limits of the theorem, the relative height of the deepest layer does not change with time.

Definition 1** (Source-boundedness).**

Layer $i\in\{1,2\}$ in a two-layer shallow-water system is source-bounded if there exists $K\in\mathds{R}$ such that the source terms satisfy $\forall h_{i}>K$,

[TABLE]

for $j=1$ and $j=2$.

Lemma 1**.**

Let $(s,d)=(1,2)$ or $(2,1)$, $\{D_{k}\}_{k\in\mathds{N}}$ be a sequence of increasing real numbers, $h_{0}$ and $f$ be scalar functions, and $\boldsymbol{q}_{1,0}$ and $\boldsymbol{q}_{2,0}$ be vector functions. Further, consider a 2LSWE system with initial conditions

[TABLE]

where layer $d$ is source-bounded and where both layer $d$ and the bottom layer (these are the same if $d=2$) have constant average density. Now let $\{(h_{1k},h_{2k},\boldsymbol{q}_{1k},\boldsymbol{q}_{2k})\}_{k\in\mathds{N}}$ be physical solutions to the 2LSWE system. If $\{(h_{sk},h_{dk}-D_{k},\boldsymbol{q}_{sk},\boldsymbol{q}_{dk}/D_{k})\}_{k\in\mathds{N}}$ converge and the first and second derivatives are uniformly bounded in the regions where they are well-defined, then

[TABLE]

Proof.

First, note that since the second derivatives are uniformly bounded, the mean value theorem implies that the first derivatives are equicontinuous. Then, since the first derivatives are also bounded, the Arzelà-Ascoli theorem gives that there is a subsequence where the first derivatives are uniformly convergent.(Dunford and Schwartz, 1957) This implies that we can interchange the order of limits and differentiation.(Rudin, 1976) From the definition of the energy in eq. 17 it is clear that all terms are non-negative. This, in addition to the fact that the energy is a convex function of the heights, means that for a system with constant bottom topography and $\overline{h_{i}}=1$ for $i\in\{1,2\}$, the energy is bounded from below by the height and momentum distributions that give $E=g\rho_{i}/2$. We let $\widetilde{E_{k}}=2E_{k}/D_{k}^{2}$ be a scaled energy, and it follows by insertion that $\widetilde{E_{k}}(0,\boldsymbol{x})\to g\rho_{d}$ as $k\to\infty$. That is, the scaled energy $\widetilde{E_{k}}(0,\boldsymbol{x})$ approaches the minimal for a system with $\overline{h_{dk}(t,\boldsymbol{x})/D_{k}}=1$ in the limit when $k\to\infty$. Source-boundedness of the mass source terms implies that $\overline{h_{dk}/D_{k}}\to 1$ for all $t\in\mathds{R}$ as $k\to\infty$.

Further, since a physical solution must satisfy the energy conservation (18), it similarly follows by use of the source-boundedness that

[TABLE]

in the limit when $k\to\infty$. Because all the terms in eq. 17 are non-negative and because the right hand side of the scaled version eq. 18 remain 0 as long as the scaled energy remains minimal, we must have

[TABLE]

Assume that $\exists\varepsilon>0$ and $\forall N\in\mathds{N}$, $\exists k>N$, such that

[TABLE]

This implies that $h_{dk}(t,\boldsymbol{x})/D_{k}$ deviates from 1 by a term which does not vanish in the limit $k\to\infty$. This contradicts eq. 24, as discussed above, so $h_{dk}(t,\boldsymbol{x})/D_{k}\to 1$ for all $t\in\mathds{R}$ and $\boldsymbol{x}\in\mathds{R}^{2}$ as $k\to\infty$. ∎

IV.2 The one-layer approximation theorem

In the following theorem, we show that in the similar limits as above, the 2LSWE may be reduced to the locally conservative 1LSWE (5) with a reduced gravity, $g\to\delta g$ with $\delta$ as defined in eq. 9. In the case where the top layer is shallow relative to the bottom layer, the bottom topography term drops out of the equation governing the top layer. As before, we use $s\in\{1,2\}$ to indicate which layer is shallow relative to the other, such that

[TABLE]

Theorem 1**.**

Let $(s,d)=(1,2)$ or $(2,1)$, $\{D_{k}\}_{k\in\mathds{N}}$ be a sequence of increasing real numbers, $h_{0}$ and $f$ be scalar functions, and $\boldsymbol{q}_{1,0}$ and $\boldsymbol{q}_{2,0}$ be vector functions, all defined on $\Omega\subseteq\mathds{R}^{n}$. Further, consider a 2LSWE in the form of eq. 7 or eq. 8 with initial conditions

[TABLE]

in which layer $d$ is source-bounded and the density of layer $d$ is constant. Now let $\{(h_{1k},h_{2k},\boldsymbol{q}_{1k},\boldsymbol{q}_{2k})\}_{k\in\mathds{N}}$ be physical solutions to the 2LSWE such that $\boldsymbol{q}_{dk}$ satisfies the boundary condition

[TABLE]

with $K\in\mathds{R}$ independent of $k$ and where $\partial\Omega$ may be at infinity.

If $\{(h_{sk},h_{dk}-D_{k},\boldsymbol{q}_{sk},\boldsymbol{q}_{dk}/D_{k})\}_{k\in\mathds{N}}$ converge and the first and second derivatives are uniformly bounded in the regions where they are well-defined, then $(h_{s},\boldsymbol{u}_{sk})\to(h,\boldsymbol{u})$ where $(h,\boldsymbol{u})$ solves eq. 25 in the weak sense, $\boldsymbol{u}_{dk}\to\boldsymbol{0}$, and $h_{dk}-D_{k}\to\left(C-\rho_{s}^{d-1}h_{s}-\rho_{2}^{d-1}b\right)/\rho_{d}^{d-1}$, where $C$ is constant in space. If the domain on which the solution is defined is infinite in range or the mass source terms are zero, then $C$ is equal to $C=\overline{\left[\rho_{s}^{d-1}h_{s}+\rho_{d}^{d-1}f+\rho_{2}^{d-1}b\right]}_{t=0}$.

Proof.

First, we note that weak solutions of eqs. 7 and 8 will be piecewise differentiable and their states on both sides of a discontinuity are connected by a Hugoniot locus. A Hugoniot locus at some location in phase space is defined as all those states for which there is a shock speed that satisfies the Rankine-Hugoniot condition.(Holden and Risebro, 2015)

To prove the theorem, it is therefore sufficient to show i) local convergence for regions where the solution is differentiable and ii) that the states that are allowed by the Hugoniot loci of the 2LSWE ((7) and (8)) converge to those of the 1LSWE (25). As before, we may interchange the order of limits and differentiation since the second derivatives are uniformly bounded.

We will first prove i). This will be done by proving that $\boldsymbol{u}_{dk}\to\boldsymbol{0}$ in the limit $k\to\infty$ by the use of the fundamental theorem of geometric calculus. The reader is referred to Doran and Lasenby (2003) for an overview of this branch of mathematics. The purpose of using this theorem is to give a way to explicitly express a vector quantity in terms of its divergence, curl and boundary conditions.

From conservation of mass and through source-boundedness and lemma 1, we get that

[TABLE]

Next, we show that also the curl of $\boldsymbol{u}_{dk}$ vanish in the limit $k\to\infty$. In the following, we use $A\wedge B$ to denote the wedge product, or exterior product, of $A$ and $B$, and $A\odot B$ to denote their geometric product. Applying $\nabla\wedge$, a generalized curl, from the left of the velocity equation of eq. 7 yields

[TABLE]

where $w_{dk}\equiv\nabla\wedge\boldsymbol{u}_{dk}$ is a bivector which is equal in magnitude to the curl of $\boldsymbol{u}_{dk}$ but well-defined in any dimension. Here $I^{-1}=\boldsymbol{e}_{m}\wedge\dotsb\wedge\boldsymbol{e}_{1}$, where $\boldsymbol{e}_{i}$ is the unit vector in direction $i$ and $m$ is the number of dimensions. Source-boundedness and the fact that $w_{dk}=0$ at $t=0$ implies that $w_{dk}\to\boldsymbol{0}$ in the limit $k\to\infty$.

From the Helmholtz theorem, we know that a vector field defined on a finite domain or which goes sufficiently fast to $\boldsymbol{0}$ is uniquely specified by its boundary condition, curl and divergence. Using techniques from geometric calculus, we can give an analytic expression. The fundamental theorem of geometric calculus states that(Doran and Lasenby, 2003)

[TABLE]

where $L$ is any linear function, $I_{m}$ is the pseudoscalar of the tangent space to $V$ and $I_{m-1}(\boldsymbol{x})$ is the pseudoscalar of the tangent space to $\partial V$ at $\boldsymbol{x}$. The vector derivative is $\nabla\odot A=\nabla\cdot A+\nabla\wedge A$, and the overdot indicates where it acts. That is, the integrand on the right hand side of eq. 29 is $\sum_{i}\partial_{i}L(\boldsymbol{e}_{i}\odot I_{m})$. See for instance the textbook by Doran and Lasenby.Doran and Lasenby (2003)

Let $V$ be some region where $u_{dk}$ is differentiable and let

[TABLE]

where $G$ is the Green’s function for the vector derivative, meaning that $\nabla\odot G(x^{\prime},x)=\delta(x-x^{\prime})$, and $S_{m-1}$ is the volume of the $(m-1)$-sphere. Equation 29 then states that

[TABLE]

The surface integral can be decomposed into one vector component whose integrand is proportional to $\boldsymbol{u}_{dk}\cdot\boldsymbol{\hat{n}}$ and one triplet-vector component whose integrand is proportional to $\boldsymbol{u}\wedge\boldsymbol{\hat{n}}$. Only the vector component will contribute to $\boldsymbol{u}_{dk}$. To show that $\boldsymbol{u}_{dk}\to 0$, it remains only to show that the last integral in eq. 31 goes to zero. This is proved in part ii) by showing that $\lim_{k\to\infty}\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{dk}$ is continuous. By applying eq. 29 with $L$ given by eq. 30 on a domain which does not include $\boldsymbol{x}$ the only contribution comes from surface integral in the limit $k\to\infty$, because $\lim_{k\to\infty}\nabla\cdot\boldsymbol{u}_{dk}+w_{dk}=0$ everywhere. $I_{m-1}$ has opposite sign on opposite sides on surfaces, so surface integrals from neighbouring domains cancel as $\lim_{k\to\infty}\boldsymbol{u}_{dk}$ is continuous. Thus, we can extend the integral over $\partial V$ to an integral over $\partial\Omega$ by applying the fundamental theorem of geometric calculus in the neighboring domains. From eq. 26 with lemma 1, we get that $\lim_{k\to\infty}\boldsymbol{u}_{dk}$ must vanish on $\partial\Omega$. Hence, $\lim_{k\to\infty}\boldsymbol{u}_{dk}=\boldsymbol{0}$ everywhere.

From the momentum equations in eq. 7, then,

[TABLE]

and so in the limit $k\to\infty$, $\rho_{s}^{d-1}h_{s}+\rho_{d}^{d-1}(h_{dk}-D_{k})+\rho_{2}^{d-1}b$ is constant in space. Finally, plugging this into the equation for $\boldsymbol{u}_{s}$ in eq. 7 or eq. 8, we get eq. 25. In the regions where the solution is differentiable, the various formulations of the 2LSWE, eqs. 6, 7 and 8, are equivalent. This completes the proof of i).

For the proof of ii), we will compare the Hugoniot loci of the 2LSWE in the limit $k\to\infty$ to the Hugoniot locus of eq. 25. Let $\gamma\vcentcolon=1/\left\langle{h_{d}}\right\rangle$. In appendix B, we show that the full set of Rankine-Hugoniot conditions for the 2LSWE eqs. 7 and 8 may be written as

[TABLE]

where

[TABLE]

For eq. 7,

[TABLE]

For eq. 8 with $d=1$,

[TABLE]

and

[TABLE]

And finally, for eq. 8 with $d=2$,

[TABLE]

and

[TABLE]

In particular, we note that $g_{1},g_{2}$, and $g_{3}$ vanish for $\gamma=0$ in all cases.

Next, we notice that eq. 33 with $\gamma=0$ is exactly the Rankine-Hugoniot relations for the locally conservative 1LSWE (25) together with the conditions that $\rho_{1}^{d-1}h_{1}+\rho_{2}^{d-1}h_{2}$ is constant and $\boldsymbol{u}_{d}=\boldsymbol{0}$. From lemma 1, it follows that $\lim_{k\to\infty}\gamma=0$. Thus the Hugoniot loci match, and this concludes the proof of the theorem. ∎

IV.3 Discussion of the theorem

Theorem 1 shows that we may approximate the thinnest layer of the 2LSWE with the locally conservative 1LSWE where $g\to(1-\rho_{1}/\rho_{2})g$ according to eq. 25. The approximation becomes more accurate when the depth of the deepest layer is increased without increasing momentum or other key properties. Figure 5 shows a sketch of how the two-layer cases converge to one-layer cases when we increase the “depth”, $D_{k}$.

The interesting part about theorem 1 is not that smooth solutions of the 2LSWEs can be approximated by solutions of the 1LSWE. It is rather that a particular form of the 1LSWE, the locally conservative form, also captures weak solutions, meaning that it gives the correct shock speeds and relations between height- and velocity-distributions at either sides of discontinuities. This is important, because while 1LSWE has been used to model two-layer spreading before, it has always been under the assumption that one must use additional equations at discontinuities in order to account for the effects of the additional layer.

A surprising implication of this result is that it suggests that the shallow water framework works better to describe shocks than one would anticipate from the assumption of negligible vertical acceleration. By analytical and experimental considerations not related to the shallow water framework, it has been found that the Froude number at the leading edge of a spreading fluid in a two-layer system lies in the range $[1,\sqrt{2}]$.(Benjamin, 1968; Ungarish, 2017; Fay, 2007; Fannelop and Waldman, 1972; Hoult, 1972; Hatcher and Vasconcelos, 2014) Theorem 1 implies that this is also true in the shallow water model. Using the Rankine-Hugoniot condition of the locally conservative 1LSWE we get $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$.

From a practical standpoint, the result presented in this paper makes it more straightforward to use the shallow-water framework to model two-layer flow with discontinuous distributions, such as oil-spills. Previous numerical schemes which have been created to ensure that the height- and velocity-distributions satisfy front conditions, which typically involves $\mathrm{Fr}_{\text{\tiny LE}}$, have had to track the position of the leading edge and alter the solution.(Hatcher and Vasconcelos, 2014, 2013) In contrast, when using the 1LSWE, which correctly captures shocks of 2LSWE, one automatically obtains numerical solutions that satisfy the Rankine-Hugoniot conditions and hence satisfies the front-condition $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$.

Finally, we remark that the mathematical tools used to prove theorem 1 are not directly applicable to the layerwise formulation of the 2LSWE (6). One way to possibly find if there is a one-layer model also for the layerwise 2LSWE is to viscously regularize the equations by an added viscosity. How viscosity looks in the shallow water framework is for instance given in.(Marche, 2007) Adding viscosity smooths out discontinuities and renders the interaction term $h_{1}\nabla h_{2}$ well-defined. The equations can then be investigated numerically by studying how shocks emerge when the viscosity coefficient is reduced. They can also be investigated analytically by looking at travelling wave solutions inside the emerging shocks.

V Cases for the one-dimensional dam-break problem

In this section we present the cases that will be used to investigate theorem 1 numerically. The cases represent variations of the one-dimensional dam-break problem.(LeVeque, 2002) In the two-layer dam-break problem, a lighter fluid of height $h_{1}$ spreads on top of a heavier fluid of height $h_{2}$ as shown in figure 6. The problem has been frequently used in the literature as a benchmark case for spreading models.(Joshi and Jaiman, 2018; Soares-Frazão et al., 2012; Zhou et al., 2004)

In the following, we first consider the dam-break problem in an unrestricted spatial domain (“Case 0”). This case will be used for convergence analyses. We next consider the dam-break problem with a reflective wall boundary-condition (“Case R” for “reflective”), which is used both to compare qualitative differences between the forms of the 1LSWE and 2LSWE and to compare results with experimental data on two-layer spreading. An overview of the cases is provided in figure 7.

V.1 Case 0: Dam-break in an unrestricted spatial domain

The initial conditions for the standard, one-dimensional dam-break problem that is not restricted in the flow-direction are

[TABLE]

where $h_{0}$ is constant.

In this particular case, the corresponding one-layer problem has self-similar analytic solutions for both variants of the 1LSWE. With the standard 1LSWE (4), there is the well-known Ritter solution,(Ritter, 1892)

[TABLE]

where $c_{0}=\sqrt{\delta gh_{0}}$. This solution is obtained from the assumption that eq. 4 is valid across discontinuities, as is normally the case when working with the 1LSWE. For the locally conservative form (5), the analytic solution is

[TABLE]

A sketch of the two solutions for $h$ is shown in figure 8. One can see that the Ritter solution expands more than $2.4$ times faster than the solution of the locally conservative form. The latter solution is the only one with a discontinuous height profile, and it has a constant Froude number of $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$ at the leading edge.

V.2 Case Ra: Quantify inaccuracies in the one-layer approximation

In Case Ra, the initial conditions are the same as for Case 0 (eq. 41). However, a reflective wall is placed to the left of the dam at position $x=-L$ with boundary conditions $(\partial_{x}h)(x=-L,t)=0$ and $u(x=-L,t)=0$. The reflective wall removes the self-similarity of the solution, which enables a study of how the accuracy of the one-layer approximation from theorem 1 evolves in time.

We also consider a variant of this case where the top layer becomes deep, i.e. $s=2$ in theorem 1. Here the initial conditions become

[TABLE]

where again $h_{0}$ is constant.

V.3 Case Rb: Effect of non-zero depth on both sides of dam

Case Rb is a variant of case Ra where the initial conditions are relaxed to allow a non-zero depth to the right of the dam, that is,

[TABLE]

In this case, the difference between the solutions of the locally and globally conservative 1LSWE will be less notable, because both give shocks. This case will be used to show that the locally conservative 1LSWE captures quantitative behaviour of two-layer cases that is not captured by the globally conservative 1LSWE.

V.4 Case Rc: Comparison to dam-break experiments

Finally, in case Rc we compare numerical results of the dam-break case with experimental results for liquid-on-liquid spreading. In particular, we compare the spreading radius predicted by the one-layer approximation from theorem 1 (eq. 25) and by the two-layer equations to two sets of experimental results. We use the same initial conditions as in case Ra, that is, eq. 41 with a reflective wall at $x=-L$.

The first set of experiments is from Suchon,Suchon (1970) who studied the spreading of oil on water. He used a $2.5\text{\,}\mathrm{m}$ long and $0.62\text{\,}\mathrm{m}$ wide channel with glass walls. The initial dam was controlled by a thin aluminum plate that was manually removed to start the experiment. We use initial conditions corresponding to 4 different runs by Suchon, see table 1. The initial depth of water in the experiments was about $30\text{\,}\mathrm{cm}$, which is nearly twice the initial heights.

The second set of experiments that will be considered are those presented by Chang, Reid, and Fay.Chang, Reid, and Fay (1983) They studied fluids at cryogenic temperatures (cryogens) spreading on water and presented both experimental results as well as model predictions. In their model, they used the same empirical boundary condition for the spreading rate as discussed in section III.(Chang and Reid, 1982) We will demonstrate that their experimental results can be reproduced to a high accuracy without any empirical boundary condition or model for the spreading rate.

It should be noted that the experimental setup by Chang, Reid, and Fay (1983) deviates from the dam-break case in that the initial reservoir of cryogen is emptied through a large slit. The spreading then occurred inside a cylinder of length $4\text{\,}\mathrm{m}$ with an inside diameter of $16.5\text{\,}\mathrm{cm}$ where half of the volume was filled with water. However, the case should be well approximated by a dam-break since the slit height is of the same order of magnitude as that of the leading edge of the spreading liquid. To the best of our knowledge, the numerical predictions by Chang, Reid, and Fay were also based on the dam-break case. Chang, Reid, and Fay do not list the initial height and width of the released cryogens, only the initial volumes. The initial conditions are therefore estimated based on the description of the apparatus given by Chang and Reid.Chang and Reid (1982) We assume that the spreading occurs in a channel of the same width as that of the experiment. We then estimated the area of the release tank and used this to find an estimate for the initial height and width from a given initial volume. The initial conditions used are listed in table 1.

To accurately represent the spreading of cryogens, it is necessary to account for evaporation due to heat flow from the water and surrounding air. The evaporation gives a source term in the mass conservation laws. We follow Chang, Reid, and Fay (1983) and include constant evaporation rates of $0.16\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ for methane and $0.201\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ for nitrogen in the mass balances. These values are the same as those used by Chang, Reid, and Fay, which are based on experimental studies of the relevant substances.(Burgess, Murphy, and Zabetakis, 1970) Evaporation leads to the formation of bubbles in the liquid, which reduces its density. Chang, Reid, and Fay called this reduced density the effective cryogenic density, and they estimated it based on experimental results to be $660\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ for nitrogen and $254\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ for methane. Similar to Chang, Reid, and Fay, we will use reduced densities in our simulations as well. Finally, Chang, Reid, and Fay report the formation of some ice on the water surface downstream of the cryogen distributor. This effect is not accounted for in the models, although it is also not expected to have a large impact on the spreading rates.

VI Numerical results

In this section, we will discuss results from the cases described in section V. The equations are discretized spatially using a finite-volume scheme. We employ the FORCE (first-order centered) flux (Toro and Billett, 2000) and the second-order MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) reconstruction with a minmod limiter (LeVeque, 2002) in each finite volume. The solutions are advanced in time with a standard third-order three-stage strong stability-preserving Runge-Kutte method.(Ketcheson and Robinson, 2005) Although some of the equations have terms that are in general not conservative, e.g. the convective term in the locally conservative 2LSWE (7), it should be noted that these become conservative when the equations are restricted to a single spatial dimension. The Courant–Friedrichs–Lewy (CFL)-number is $0.9$ for all cases.

For quantification of errors, we use the $L^{1}$ norm, which for a function $y:\Omega\to\mathds{R}$ is defined as

[TABLE]

We present numerical results for the cases described in section V for $\delta$ of $0.6$ and $0.7$, which are similar to those of cryogenic spills on water (see for instance table 1). To find the height and velocity distributions in the one-layer model for other values of $\delta$, a temporal scaling is all that is needed. This is because the equations are invariant under the transformations $\delta g\to\lambda\delta g$, $t\to\sqrt{\lambda}^{-1}t$ and $u_{s}\to\sqrt{\lambda}u_{s}$ for all $\lambda$. We find that the convergence is also similar for other values of $\delta$.

VI.1 Case 0: Dam-break in an unrestricted spatial domain

In Case 0, the dam-break occurs in a one-dimensional, spatially unrestricted domain. We solved this case with the parameters $h_{0}=$1\text{,}\mathrm{m}$$ and $\delta=0.7$ with 400 grid cells on the domain $x/c_{0}t\in(-2,2)$. The initial depth, $H$, was increased stepwise from the initial value, $H=h_{0}$, to obtain a larger height difference between the layers. The results presented in Figure 9 show that the locally conservative 2LSWE converge towards the analytic solution of the locally conservative 1LSWE when $H$ increases. A higher value of $H$ translates into increasing $D_{k}$ in theorem 1. The figure demonstrates a general trend found for the agreement between the 1LSWE and the 2LSWE, namely that the locally conserved 1LSWE become an increasingly good approximation to the complete 2LSWE with increasing $H$.

Figure 10 shows the normalized difference in $L^{1}$ for the top-layer height between the analytic solution (eq. 43), and the solutions obtained with the 2LSWE, $\tilde{h}_{1}$ and $h_{1}$, respectively ($\lVert h_{1}-\tilde{h}_{1}\rVert_{1}/\left\lVert{h_{0}}\right\rVert_{1}$). The differences are shown as a function of the initial depth, $H$. The circles correspond to the locally conservative 2LSWE (7), the triangles correspond to the globally conservative 2LSWE (8), and the solid line indicates a slope of -1. The plot shows that the difference between the globally and the locally conservative 2LSWE is small as expected. Other relevant variables such as $h_{2}$, $u_{1}$, and $u_{2}$, were found to exhibit a similar behaviour.

VI.2 Case Ra: Quantify inaccuracies in the one-layer approximation

In Case Ra, a reflective wall is placed at $x=-L$, cf. section V.2. The case was solved with the parameters $h_{0}=$1\text{,}\mathrm{m}, $L=$2\text{\,}\mathrm{m}, and $\delta=0.6$. The domain width was $15\text{\,}\mathrm{m}$ and 1000 grid cells were used.

We first compare the two situations $s=1$ (the top layer is shallow) and $s=2$ (the bottom layer is shallow), cf. theorem 1. Figure 11 shows a dam break in the left column ($s=1$) and the cross section of a gravity current in the right column ($s=2$) at $t=$3\text{,}\mathrm{s}$$. An illustration of the initial configurations is presented at the top of the figure. The globally conservative 2LSWE (green lines) are compared with the locally conservative 1LSWE (green lines) and the globally conservative 1LSWE (red lines). For the dam break case (right column), we initialize the bottom layer in a perturbed state where the depth $h_{2}$ is constant. This is done to show that a perturbation of the initial solution of the relatively deep layer does not prevent the 2LSWE to converge to the one-layer approximation when the depth increases. As the depths increase, we see that the solutions of the 2LSWE converge toward the locally conservative 1LSWE as predicted by the theorem.

To quantify how the solutions of the 2LSWE converge to those of the locally conservative 1LSWE, we will compare the solutions at various times, initial depths, $H$, and initial widths, $L$. We consider solutions of the globally conservative 2LSWE and the locally conservative 1LSWE and evaluate two quantifiable differences. In figure 12(a), the $L^{1}$ difference of the top-layer height, $\lVert h_{1}^{\mathrm{2LSWE}}-h_{1}^{\mathrm{1LSWE}}\rVert_{1}/h_{0}L$ is plotted for varying times and depths, $H/h_{0}$. Figure 12(b) shows the difference of the leading edge position at $t=$5\text{,}\mathrm{s}$$, $|r_{\mathrm{2LSWE}}-r_{\mathrm{1LSWE}}|/h_{0}$, for varying initial depths, $H/h_{0}$, and widths, $L/h_{0}$. The position of the leading edge is here defined as the smallest $x$-value where the top layer is thinner than ${10}^{-4}\text{\,}\mathrm{m}$. Figure 12 shows that the differences in $h_{1}$ decrease with time, which is reasonable since the spreading fluid becomes gradually thinner. As expected, the differences decrease with increasing value of $H/h_{0}$. Similar to Case 0, the errors in the variables $h_{2}$, $u_{1}$, and $u_{2}$ as quantified by the $L^{1}$ norm exhibit the same trends as the top layer height (not shown). Further, the figure shows that the difference of the leading edge position decreases with decreasing width, which is reasonable because the spreading fluid becomes thinner as the initial volume decreases.

It is also interesting to see how the rate of spreading evolves with increasing depth, $H$. Figure 13 shows the leading edge position as a function of time for the two variants of the 1LSWE and the globally conservative 2LSWE for different depths $H$. Again we observe a rapid convergence of the 2LSWE to the locally conservative 1LSWE. Also for this case, we observe that the spreading rate from the globally conservative 1LSWE is higher that that from the full 2LSWE.

VI.3 Case Rb: Effect of non-zero depth at both sides of dam

In Case Rb, the dam-break was initialized according to eq. 45 with a non-zero depth at both sides of the dam. The case was solved with the parameters $\delta=0.6$, $L=$4\text{,}\mathrm{m}, $H=$50\text{\,}\mathrm{m}, $h_{0,a}=$2\text{,}\mathrm{m}, and $h_{0,b}=$0.5\text{\,}\mathrm{m}. The width of the domain is $15\text{\,}\mathrm{m}$ and the results are again computed with 1000 grid cells.

Figure 14 shows the height distributions at time $t=$3\text{,}\mathrm{s}$$ for both the globally and locally conserved 1LSWE and for both versions of the 2LSWE, eqs. 7 and 8 at different depths $H$. The solutions of both formulations of the 1LSWE have shocks, but the shock velocities differ. As expected, the height profiles of the locally and globally conservative 2LSWE are similar, however, they are only in agreement with the locally conservative 1LSWE. The figure shows that the locally conservative 1LSWE should be used for accurate representation of the position of the leading edge.

VI.4 Case Rc: Comparison to spreading experiments

Figures 15 and 16 present a comparison of the spreading rates predicted from the one-layer approximation from theorem 1 (green solid line) with the experiments described in section V.4 (symbols) and with the full 2LSWE (blue dashed lines). All cases were solved with 2000 grid cells. A comparison to the simpler Fay model (eq. 3) with $\beta=1.31$ is also included for the Suchon experiments (orange dotted lines).

We find good agreement between the one-layer approximation and available experimental data. The deviation in the spreading radius was calculated as the relative difference in the $L^{1}$ norm, that is,

[TABLE]

The deviation in the spreading radius is $4.5\text{\,}\mathrm{\char 37\relax}$ for oil on water, which is significantly better than the Fay model, where the average deviations are $12.6\text{\,}\mathrm{\char 37\relax}$. For the cryogenic fluids, the deviations in the spreading radius were $10.2\text{\,}\mathrm{\char 37\relax}$ for methane on water, and $4.2\text{\,}\mathrm{\char 37\relax}$ for nitrogen on water.

We remark that in the experiments, the depth of the water is not much larger than the initial depth of the spreading liquid. In the experiments by Suchon, it is about twice the initial height of the oil, and in the cryogen experiments by Chang, Reid, and Fay, it is about the same as the initial cryogen height. However, it was shown in figure 12(b) that the difference between the predicted spreading distance from the 2LSWE and the one-layer approximation after $5\text{\,}\mathrm{s}$ is still small, even at these initial depths. In all cases, the initial width $L$ is less than half the initial width, and so we expect a deviation at $5\text{\,}\mathrm{s}$ that is smaller than 0.2 times the initial height $h_{0}$.

A comparison of the green solid lines (1LSWE) and the blue dashed lines (2LSWE) in figure 15 confirms a very good agreement between the two formulations. For figure 16, and especially for the $2\text{\,}\mathrm{L}$ cases where the initial height is large compared to the water depth, the discrepancy is larger. We find that the deviation in the $L^{1}$ norm between the one-layer approximation and 2LSWE results is between $0.7\text{\,}\mathrm{\char 37\relax}$ and $5.5\text{\,}\mathrm{\char 37\relax}$ for all cases. In the $2\text{\,}\mathrm{L}$ cryogen experiments, we observe that the discrepancy is reducing after some time. This is consistent with figure 12(a). The evaporation that occurs during the spreading of cryogenic fluids likely accelerates the decrease in error.

The results indicate that the proposed one-layer approximation may be used as an approximation to the 2LSWE even for cases where the depth ratio is small, as long as the main interest is to predict the spreading distance. Although, in these cases one should not expect that the one-layer approximation captures all of the qualitative flow patterns that are captured by the 2LSWE. In figure 17, we compare the profiles of the top layer at three different times for the 1LSWE (green solid lines) and the 2LSWE (blue dashed lines). The observed spreading distance from the corresponding experiments are marked by red dots. We see that the 2LSWE captures a more complex behaviour, especially in the early phase of the flow, but as expected, the agreement between the profiles improves with time.

Finally, we note that Chang, Reid, and Fay (1983) also solved the 1LSWE numerically, but with an imposed boundary condition with $\mathrm{Fr}_{\text{\tiny LE}}=1.28$ at the leading edge. They motivated the use of a constant Froude number at the leading edge by frequent use in previous literature dealing with spreading of non-boiling fluids. They treated the value of $\mathrm{Fr}_{\text{\tiny LE}}$ as a parameter that depends on the apparatus and must be determined experimentally and found that $\mathrm{Fr}_{\text{\tiny LE}}=1.28$ worked best for their apparatus. The analysis in this paper shows that the constant Froude number can in fact be derived from the 2LSWE. That is, as the relative depth of the bottom layer increases, $\mathrm{Fr}_{\text{\tiny LE}}$ approaches $\sqrt{2}$ from below. Moreover, our analysis shows that $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$ is true also for boiling liquids and even when including other relevant source terms in the 2LSWE model.

VII Conclusions

We have presented a comprehensive study of two-layer spreading where the depth of one layer is significantly larger than the other. The main result is that the two-layer shallow-water equations can be approximated by an effective one-layer model with an effective gravitational constant as described in theorem 1. In the literature, the globally conservative one-layer momentum equations are frequently used. We have demonstrated both analytically and numerically that the locally conservative momentum equations should be used instead for a precise representation with an effective one-layer model.

Earlier works in the literature have made use of an additional boundary condition for the speed of the leading edge as a closure relation for the one-layer spreading model. The speed is typically represented in terms of a constant Froude number which is adjusted to match experimental data so that $\mathrm{Fr}_{\text{\tiny LE}}\in[1,\sqrt{2}]$. We have shown that this boundary condition can in fact be derived from the full two-layer shallow water equations. By using the locally conservative version of the one-layer shallow water equations with the effective gravitational constant $(1-\rho_{1}/\rho_{2})g$, the one-layer model correctly captures the behaviour of shocks and contact discontinuities. In particular, the one-layer model results in $\mathrm{Fr}_{\text{\tiny LE}}=\sqrt{2}$, which is exactly the same as the theoretical predictions by von Kármán (1940); Benjamin (1968); Ungarish (2017) in the limit captured by theorem 1.

By using the same mathematical tools that were used to derive the one-layer approximation in theorem 1, we derived an expression for the Froude number at the front of a spreading fluid inside a rectangular cavity from the full two-layer shallow water equations. The expression that we obtained from the analysis of the shallow-water equations is in good agreement with the expression by Benjamin. The agreement between these expressions suggests that the validity breakdown of the shallow-water equations in vicinity of shocks is less severe than previously suggested.

We compared to available experimental data for one-dimensional dam break experiments and found good agreement between the one-layer model derived in this work and experiments, where the mean relative deviation in the spreading radius was $4.5\text{\,}\mathrm{\char 37\relax}$ for oil on water, $10.2\text{\,}\mathrm{\char 37\relax}$ for methane on water, and $4.2\text{\,}\mathrm{\char 37\relax}$ for nitrogen on water. The spreading radius from the one-layer and two-layer descriptions could hardly be distinguished from each other after 10 seconds of spreading, but the fluid profiles from the two formulations differed at short times. In comparison, the mean relative deviation in the spreading radius of the Fay model was $12.6\text{\,}\mathrm{\char 37\relax}$ for oil on water.

The treatment in this paper has also included source terms, as long as they are source-bounded. Source terms representing Coriolis forces are not source-bounded as defined in section IV because they are proportional to the depth. Thus they are not covered by the present analysis. It should be possible to include Coriolis-like source terms in the analysis, because although they are proportional to the depth, they are also proportional to the flow velocity which vanishes with increasing depth in the deep layer. Since Coriolis forces are relevant particularly for the modelling of geophysical phenomena, it represents an attractive possibility for future work.

VIII Acknowledgements

The authors wish to acknowledge fruitful discussions with Hans Langva Skarsvåg and Svend Tollak Munkejord. This work was undertaken as part of the research project “Predicting the risk of rapid phase-transition events in LNG spills (Predict-RPT)”, and the authors would like to acknowledge the financial support of the Research Council of Norway under the MAROFF programme (Grant 244076/O80).

Appendix A Deriving Rankine-Hugoniot Conditions

To obtain Rankine-Hugoniot condition for the 2LSWE we use an approach similar to that of Smoller.Smoller (1983) Consider a shock along $\Gamma$ which is normal to $\boldsymbol{\hat{n}}$ and let $\phi$ be a test function with compact support $D$ which lies in the $x_{n}t$-plane.

Consider first the locally conservative system. Because the velocity $\boldsymbol{u}_{i}$ is a weak solution, the normal component must satisfy

[TABLE]

where the subscript on $\phi$ denote partial differentiations and

[TABLE]

$u_{i}^{T}$ is the tangential component of $\boldsymbol{u}_{i}$ and $\partial_{T}$ is differentiation with respect to the tangential direction. The integrand in eq. 47 is a normal function and hence we can seperate the integral into two, one for each region where the velocities and heights are differentiable. Call these regions $D_{1}$ and $D_{2}$. Because the solution is differentiable inside these regions we may use Green’s theorem and obtain

[TABLE]

because $\phi$ vanish on the boundary of $D$. Equation 47 is obtained by adding eq. 49 with $j=1$ and $j=2$. Because eq. 47 must hold for all test functions we obtain the Rankine-Hugoniot condition for the normal velocity component,

[TABLE]

For the globally conservative system a similar treatment yields

[TABLE]

and

[TABLE]

The treatment presented above is not applicable for the tangential component of the velocity equations because they involve a term on the form $\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{i}\partial_{n}u_{i}^{T}$. Note, that the tangential velocity components does not enter any of the other Rankine-Hugoniot conditions and that the equations are consistent if the tangential component are continuous across shocks. Requiring \left\llbracket{u_{i}^{T}}\right\rrbracket=0 has the additional advantage of making $\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{i}\partial_{n}u_{i}^{T}$ well defined as the product of $\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{i}$ and $\partial_{n}u_{i}^{T}$. Otherwise the distribution $\boldsymbol{\hat{n}}\cdot\boldsymbol{u}_{i}\partial_{n}u_{i}^{T}$ can not be decomposed without relying on some mollification scheme.

Ostapenko (2001) proposed that in the two-dimensional case one can use a Rankine-Hugoniot condition for the vorticity instead of the tangential velocity component. Unfortunately, the proposed equation works only if one assumes \left\llbracket{u_{i}^{T}}\right\rrbracket=0. A conservation law for the quantity $\partial_{1}u_{i,2}-\partial_{2}u_{i,1}\equiv w_{i}$ can be obtained by taking distributional derivatives of the different components of the velocity equation, yielding

[TABLE]

where

[TABLE]

It may be tempting from eq. 53 to conclude that the vorticity must obey the jump condition

[TABLE]

However, this is only true if the vorticity $w_{i}$ can be interpreted as a normal function. If the tangential velocity component is discontinuous across the shock, the vorticity would have a contribution similar to a delta distribution at the shock. In that case one can not seperate the integral into two as was done in the derivation above, and the Rankine-Hugoniot condition would gain an additional contribution from the delta-like term.

For the purposes of this paper we can ignore the tangential velocity components across jumps. The relevant equation in the one layer system is equal in both the globally conservative two-layer system and in the locally conservative two-layer system in the relevant limits. Solutions of the two-layer systems are therefore also solutions of the locally conservative one-layer system.

Appendix B Rankine-Hugoniot conditions for 2LSWE

In this appendix, we will show that the Rankine-Hugoniot conditions for the 2LSWE may be written as eq. 33, repeated here for convenience,

[TABLE]

Here $g_{1}$ and $g_{2}$ differ for eq. 7 and eq. 8, while $g_{3}$ will be the same. Further, we will show that all of $g_{1},g_{2}$, and $g_{3}$ vanish when $\gamma=0$. Note that eq. 56a follows directly from eq. 11 applied to the shallowest layer.

We first consider $g_{3}$, which can be obtained from mass conservation of layer $d$. The scalar Rankine-Hugoniot condition (11) immediately yields

[TABLE]

where we used that \left\llbracket{ab}\right\rrbracket=\left\llbracket{a}\right\rrbracket\left\langle{b}\right\rangle+\left\langle{a}\right\rangle\left\llbracket{b}\right\rrbracket. This gives

[TABLE]

Next we consider the expressions for $g_{1}$ and $g_{2}$. We first consider the locally conservative momentum equations ((7c) and (7d)). We apply the Rankine-Hugoniot condition (14) with $i=d$ and insert eq. 58 to obtain

[TABLE]

that is,

[TABLE]

Now consider eq. 14 with $i=s$,

[TABLE]

If we consider the cases $s=1$ and $s=2$ separately and use that $d-1=2-s$, we find from eq. 59 that

[TABLE]

where $\delta=1-\rho_{1}/\rho_{2}$, as defined in eq. 9. Inserting into eq. 60 we get that

[TABLE]

with

[TABLE]

Finally, we consider the globally conservative momentum equations ((8c) and (8d)). We first consider the case where $d=2$. If we take the scalar product of the Rankine-Hugoniot conditions eqs. 16 and 15 with $\boldsymbol{\hat{n}}$ and use eq. 58, we obtain

[TABLE]

that is,

[TABLE]

and

[TABLE]

Next we consider the case when $d=1$. We may then write eq. 15 as

[TABLE]

We then insert this into eq. 16 to get

[TABLE]

which gives

[TABLE]