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

Fyhn, Eirik Holm; Lervåg, Karl Yngve; Ervik, Åsmund; Wilhelmsen, Øivind

doi:10.1063/1.5126168

The gravity-driven spreading of one fluid in contact with another fluid is of key importance to a range of topics. These phenomena are commonly described by the two-layer shallow-water equations (SWE). When one layer is significantly deeper than the other, it is common to approximate the system with the much simpler one-layer SWE. It has been assumed that this approximation is invalid near shocks, and one has applied additional front conditions to correct 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 capture the behavior of shocks and contact discontinuities without additional closure relations. The result shows that simplification to an effective one-layer model is justified mathematically and can be made without additional knowledge of the shock behavior. The shock speed in the proposed model is consistent with empirical models and identical to front conditions that have been found theoretically by von Kármán and Benjamin. This suggests that the breakdown of the SWE in the vicinity of shocks is less severe than previously thought. We further investigate the applicability of the SW framework to shocks by studying one-dimensional lock-exchange/-release. We derive expressions for the Froude number that are in good agreement with the widely employed expression by Benjamin. The equations are solved numerically to illustrate how quickly the proposed model converges to solutions of the full two-layer SWE. We also compare numerical results from the model with results from experiments and find good agreement.

I. INTRODUCTION

The spreading of two layers of fluids with different densities is of considerable importance. It has been an active field of study since at least 1774, when Franklin, Brownrigg, and Farish¹ 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^2–4 and liquefied gaseous fuels,^5–7 stratified flow inside pipes,⁸ gravity currents particularly in geophysical systems,^9–12 monomolecular layers for evaporation control,¹³ and coalescence in three-phase fluid systems.¹⁴ These applications include immiscible fluids such as oil and water or systems with miscible fluids at a 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,^15,16

F r = \frac{u}{\sqrt{g^{'} h}},

(1)

where u is the velocity, h is the height of the layer that is spreading, and g′ is the effective gravitational acceleration. In two layer spreading, the effective gravitational acceleration is g′ = (1 − ρ₁/ρ₂)g, where ρ₁ and ρ₂ are the two fluid densities and ρ₁ < ρ₂.

An early result for the Froude number of gravity currents was presented by von Kármán.¹⁷ They found that for the edge of a spreading of gravity current at semi-infinite depth, ${F r}_{LE} = \sqrt{2}$ ⁠, where the subscript is short for the “leading edge.” Benjamin¹⁸ later developed a model for Fr_LE for spreading of gravity currents with constant height,

{F r}_{LE}^{2} = \frac{(1 - α) (2 - α)}{(1 + α)},

(2)

where α = h₂/(h₁ + h₂). Here, h₁ and h₂ 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₂ ≪ h₁. More recently, Ungarish¹⁹ extended the result of Benjamin to the spreading of gravity currents into a lighter fluid with an open surface. This result also gives ${F r}_{LE} = \sqrt{2}$ when the spreading fluid becomes 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,²⁰ 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

u_{LE} = β \sqrt{\frac{g^{'} V}{A}},

(3)

where β 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, β represents an effective Froude number where the height at the leading edge is approximated by the average height. The value of β has been discussed in the literature and is commonly set to β = 1.31 in the one-dimensional case and β = 1.41 in the axisymmetric case.^2,3,5,21

A more general approach than the Fay model is the two-layer shallow-water equations (2LSWE) that are derived from the Euler equations by assuming a negligible vertical velocity.^22,23 These equations model the flow of two layers of shallow liquids and may be used to simulate, for instance, gravity currents.²⁴ 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 nonhyperbolic state.²⁵ A breakdown of hyperbolicity causes problems such as ill-posedness and Kelvin-Helmholtz like instabilities.^26–28 Nonhyperbolic equations are generally more difficult to analyze and computationally much more expensive to solve than hyperbolic equations.²⁹ Attempts to amend the nonhyperbolicity of the systems include adding numerical (nonphysical) friction forces,³⁰ operator-splitting approaches,³¹ and introduction of an artificial compressibility.³²

Due to their comparative simplicity, the one-layer shallow-water equations (1LSWE) have often been used to model two-layer phenomena such as 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.^2,21,33 The additional boundary condition at the leading edge has also been used in combination with the 2LSWE.^34,35 In particular, Rottman and Simpson³⁴ 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.^36–38 There is an exhaustive literature on the subject of hyperbolicity of the 1LSWE,^39–44 including the topic of well-balanced formulation, the more general E-balanced schemes, and the identification of resonant vs nonresonant 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 behaviors 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 Fig. 1, which illustrates how solutions to different forms of the 1LSWE compare to the solution of the 2LSWE for a dam-break problem. This figure shows a clear difference between the locally and globally conservative 1LSWE.

FIG. 1.

View large Download slide

An example of how solutions from different formulations of the one-layer shallow-water equations (1LSWE local and 1LSWE global) compare to those from the two-layer shallow-water equations (2LSWE) for a dam-break problem.

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 with previous results from the literature. This indicates that the breakdown of the shallow-water equations in the vicinity of shocks is less severe than previously suggested.

This paper is structured as follows: In Sec. 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 Sec. III, we derive expressions for the Froude number from the full two-layer shallow-water equations. The key result of this paper is presented in Sec. IV, where we show that 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 Sec. V, we define some numerical experiments that are used in Sec. 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 Sec. 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 Fig. 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.

FIG. 2.

View large Download slide

A sketch of a general two-layer shallow-water geometry.

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 traveling wave solutions to obtain physically acceptable solutions of the partial differential equations (PDEs). At the end of this section, we present a necessary energy requirement for the 2LSWE that is used to select correct physical solutions in Sec. IV.

A. The one-layer shallow-water equations

The 1LSWE are typically presented in a globally conservative form where the total momentum is conserved,⁴⁵

\frac{\partial}{\partial t} ρ h + \nabla \cdot (ρ h u) = G_{h},

(4a)

\frac{\partial}{\partial t} (ρ h u) + \nabla \cdot (ρ h u \otimes u) + \nabla (\frac{1}{2} g ρ h^{2}) = G_{h u} - g ρ h \nabla b,

(4b)

where ρ is the density, h is the height, u is the vertically averaged horizontal velocity, ⊗ denotes the tensor product, b is the bottom topography, and G_h and 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 ρ 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,

\frac{\partial}{\partial t} ρ h + \nabla \cdot (ρ h u) = G_{h},

(5a)

\frac{\partial}{\partial t} u + (u \cdot \nabla) u + g \nabla (h + b) = \frac{1}{ρ h} (G_{h u} - u G_{h}) .

(5b)

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 Eqs. (5) and (4) is the admissibility of shocks when the height drops to 0. This will be further discussed in Sec. V A, but the upshot is that such a shock is impossible in Eq. (4), while in Eq. (5), it is possible with a Froude number ${F r}_{LE} = \sqrt{2}$ ⁠. This is exactly the result by von Kármán¹⁷ for two layer spreading with such shocks. In fact, in Sec. IV, we show that the locally conservative form correctly captures the two-layer behavior in certain limits. This result is consistent with previous results, which show that numerical approaches will fail to solve the conservation of global momentum.³¹

B. The two-layer shallow-water equations

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

\frac{\partial}{\partial t} ρ_{1} h_{1} + \nabla \cdot (ρ_{1} h_{1} u_{1}) = G_{h_{1}},

(6a)

\frac{\partial}{\partial t} ρ_{2} h_{2} + \nabla \cdot (ρ_{2} h_{2} u_{2}) = G_{h_{2}},

(6b)

\begin{align} \frac{\partial}{\partial t} (ρ_{1} h_{1} u_{1}) + \nabla \cdot (ρ_{1} h_{1} u_{1} \otimes u_{1}) + \nabla (\frac{1}{2} g ρ_{1} h_{1}^{2}) \\ = G_{h_{1} u_{1}} - g ρ_{1} h_{1} \nabla (h_{2} + b), \end{align}

(6c)

\begin{align} \frac{\partial}{\partial t} (ρ_{2} h_{2} u_{2}) + \nabla \cdot (ρ_{2} h_{2} u_{2} \otimes u_{2}) + \nabla (g ρ_{1} h_{1} h_{2} + \frac{1}{2} g ρ_{2} h_{2}^{2}) \\ = G_{h_{2} u_{2}} + g ρ_{1} h_{1} \nabla (h_{2} + b) - g (ρ_{1} h_{1} + ρ_{2} h_{2}) \nabla b, \end{align}

(6d)

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

This form was originally described by Ovsyannikov²² and is referred to in more recent works as “the conventional two-layer shallow-water model.”³²

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

Conservation of momentum can be considered at three different scales:

The globally conservative form where the total momentum is conserved.
The layerwise conservative form [Eq. (6)] where the momentum in each layer is conserved.
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 Sec. II D. 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,

\nabla \cdot (ρ_{i} h_{i} u_{i} \otimes u_{i}) = u_{i} \nabla \cdot (ρ_{i} h_{i} u_{i}) + ρ_{i} h_{i} (u_{i} \cdot \nabla) u_{i},

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

\frac{\partial}{\partial t} ρ_{1} h_{1} + \nabla \cdot (ρ_{1} h_{1} u_{1}) = G_{h_{1}},

(7a)

\frac{\partial}{\partial t} ρ_{2} h_{2} + \nabla \cdot (ρ_{2} h_{2} u_{2}) = G_{h_{2}},

(7b)

\frac{\partial}{\partial t} u_{1} + (u_{1} \cdot \nabla) u_{1} + \nabla [g (h_{1} + h_{2} + b)] = \frac{1}{ρ_{1} h_{1}} (G_{h_{1} u_{1}} - u_{1} G_{h_{1}}),

(7c)

\frac{\partial}{\partial t} u_{2} + (u_{2} \cdot \nabla) u_{2} + \nabla [g (\frac{ρ_{1}}{ρ_{2}} h_{1} + h_{2} + b)] = \frac{1}{ρ_{2} h_{2}} (G_{h_{2} u_{2}} - u_{2} G_{h_{2}}) .

(7d)

For a comprehensive study of the well-posedness of the locally conservative 2LSWE, see, for instance, Ref. 46.

When conserving the total momentum, the sum of Eqs. (6c) and (6d) is used, which has the advantage of eliminating the interaction between the layers. However, this approach requires an additional conservation law. Ostapenko^47,48 showed that the additional conservation law should be the difference between Eqs. (7d) and (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

\frac{\partial}{\partial t} ρ_{1} h_{1} + \nabla \cdot (ρ_{1} h_{1} u_{1}) = G_{h_{1}},

(8a)

\frac{\partial}{\partial t} ρ_{2} h_{2} + \nabla \cdot (ρ_{2} h_{2} u_{2}) = G_{h_{2}},

(8b)

\begin{align} \frac{\partial}{\partial t} (ρ_{1} h_{1} u_{1} + ρ_{2} h_{2} u_{2}) + \nabla \cdot (ρ_{1} h_{1} u_{1} \otimes u_{1} + ρ_{2} h_{2} u_{2} \otimes u_{2}) \\ + \nabla (\frac{1}{2} g ρ_{1} h_{1}^{2} + ρ_{1} g h_{1} h_{2} + \frac{1}{2} ρ_{2} g h_{2}^{2}) \\ = G_{h_{1} u_{1}} + G_{h_{2} u_{2}} - g (ρ_{1} h_{1} + ρ_{2} h_{2}) \nabla b, \end{align}

(8c)

\frac{\partial}{\partial t} (u_{2} - u_{1}) + (u_{2} \cdot \nabla) u_{2} - (u_{1} \cdot \nabla) u_{1} - \nabla (g δ h_{1}) = J,

(8d)

where

J = \frac{G_{h_{2} u_{2}} - u_{2} G_{h_{2}}}{ρ_{2} h_{2}} - \frac{G_{h_{1} u_{1}} - u_{1} G_{h_{1}}}{ρ_{1} h_{1}}

and where we have defined

δ : = \frac{ρ_{2} - ρ_{1}}{ρ_{2}} .

(9)

D. 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.^49–51 In the 2LSWE, Eqs. (6)–(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.^52–54

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

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

(10)

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

S 〚 u 〛 = \hat{n} \cdot 〚 q 〛 .

(11)

Similarly, if u satisfies a general vector conservation equation,

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

(12)

then if there is some discontinuity in u, we have the result

S 〚 u 〛 = 〚 \hat{n} \cdot (a \otimes b) + q \hat{n} 〛 .

(13)

Equations (11) and (13) can be directly applied to the mass conservation equations and the conservation law for the 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 u · ∇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 ∝h₁∇h₂ 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)],

S 〚 u_{i} 〛 \cdot \hat{n} = 〚 \frac{1}{2} {(\hat{n} \cdot u_{i})}^{2} + g {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} h_{1} + g h_{2} 〛,

(14)

where as before i = 1, 2 denotes the layer. Similarly, for the globally conservative 2LSWE [Eqs. (8c) and (8d)], we find

\begin{matrix} S 〚 ρ_{1} h_{1} u_{1} + ρ_{2} h_{2} u_{2} 〛 & = & 〚 (\hat{n} \cdot u_{1}) ρ_{1} h_{1} u_{1} + (\hat{n} \cdot u_{2}) ρ_{2} h_{2} u_{2} 〛 \\ + 〚 \frac{1}{2} g ρ_{1} h_{1}^{2} + ρ_{1} g h_{1} h_{2} + \frac{1}{2} ρ_{2} g h_{2}^{2} 〛 \hat{n} \end{matrix}

(15)

and

S \hat{n} \cdot 〚 u_{2} - u_{1} 〛 = 〚 \frac{1}{2} [) {(\hat{n} \cdot u_{2})}^{2} - {(\hat{n} \cdot u_{1})}^{2}]) - g δ h_{1} 〛 .

(16)

Finally, we note that in calculations with the Rankine-Hugoniot condition it is useful to observe that $〚 a b 〛 = 〚 a 〛 ⟨b⟩ + ⟨a⟩ 〚 b 〛$ ⁠, where $⟨a⟩ = (a_{l} + a_{r}) / 2$ ⁠.

E. 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.⁵⁰ However, the word entropy is typically restricted to convex functions of the solution variables. As has been showed by Ostapenko,⁴⁷ 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

\begin{matrix} E & = & \frac{1}{2} (ρ_{1} h_{1} {|u_{1}|}^{2} + ρ_{2} h_{2} {|u_{2}|}^{2}) \\ + g [) ρ_{2} h_{2} (\frac{1}{2} h_{2} + b) + (\frac{1}{2} h_{1} + h_{2} + b) ρ_{1} h_{1}]) . \end{matrix}

(17)

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 also be fulfilled by weak, discontinuous solutions. We obtain

\begin{matrix} \frac{\partial E}{\partial t} + \nabla \cdot [q_{1} (g (h_{1} + h_{2} + b) + \frac{1}{2} {|u_{1}|}^{2})) \\ + (q_{2} (g (\frac{ρ_{1}}{ρ_{2}} h_{1} + h_{2} + b) + \frac{1}{2} {|u_{2}|}^{2})] \\ \leq u_{1} \cdot G_{h_{1} u_{1}} + u_{2} \cdot G_{h_{2} u_{2}} - \frac{1}{2} g h_{1}^{2} \frac{\partial ρ_{1}}{\partial t} \\ - g h_{2} (\frac{ρ_{1}}{ρ_{2}} h_{1} + \frac{1}{2} h_{2}) \frac{\partial ρ_{2}}{\partial t} + G_{h_{1}} (g (h_{1} + h_{2} + b)) (- \frac{1}{2} {|u_{1}|}^{2}) \\ + G_{h_{2}} (g (\frac{ρ_{1}}{ρ_{2}} h_{1} + h_{2} + b) - \frac{1}{2} {|u_{2}|}^{2}), \end{matrix}

(18)

where q_i = ρ_ih_iu_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 (Fr_LE) of two-layer systems with the fixed total height.

Shock speeds in two-layer systems with the fixed total height are 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 Fig. 3. Such problems have been studied extensively, and there are a large number of results from laboratory experiments available.^10,34,55 Moreover, much theoretical work has been carried out to model the Froude-number for flows inside rectangular channels,^18,56,57 which means that this is a good candidate for testing the credibility of shock behavior in the 2LSWE.

FIG. 3.

View large Download slide

A sketch of the initial condition for the lock-exchange problem: Two-layer shallow-water flow in a rectangular channel. The gray fluid is lighter than the blue, and the initial shock is the vertical line between blue and gray.

Most previous works have focused on fluids with similar densities such that δ ≪ 1 for δ given by Eq. (9). This is referred to as the Boussinesq case.⁵⁸ The most commonly used front condition applied to such flows is the equation for the Froude number given by Benjamin¹⁸ [Eq. (2)]. For instance, Ungarish⁵⁹ applied the Froude number by Benjamin as a boundary condition when solving the 2LSWE for rectangular geometries. They also generalized this to arbitrary geometries.³⁵

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₁ + h₂ = H. In this case, there are no free surfaces. We therefore add an additional pressure term p₀ 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,

\frac{\partial h_{1}}{\partial t} + \frac{\partial}{\partial x} (h_{1} u_{1}) = 0,

(19a)

\frac{\partial h_{2}}{\partial t} + \frac{\partial}{\partial x} (h_{2} u_{2}) = 0,

(19b)

\frac{\partial u_{1}}{\partial t} + \frac{\partial}{\partial x} (\frac{1}{2} u_{1}^{2} + g h_{2} + \frac{1}{ρ_{1}} p_{0}) = 0,

(19c)

\frac{\partial u_{2}}{\partial t} + \frac{\partial}{\partial x} (\frac{1}{2} u_{2}^{2} + g h_{2} + \frac{1}{ρ_{2}} p_{0}) = 0 .

(19d)

These are the same equations that were used by Rottman and Simpson³⁴ 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 Fr_LE can be obtained from Eq. (19) directly.

With h₁ + h₂ constant, the sum of Eqs. (19a) and (19b) implies that h₁u₁ + h₂u₂ is constant in x. If we assume that the total momentum is 0 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₁u₁ + h₂u₂ = 0.

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

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

where, as before, the subscript l indicates the left side of the shock. We next apply the Rankine-Hugoniot condition to ρ₂(19d) − ρ₁(19c), which gives

S (ρ_{2} S + ρ_{1} \frac{h_{2, l}}{h_{1, l}} S) = \frac{1}{2} (ρ_{2} S^{2} - ρ_{1} \frac{h_{2, l}^{2}}{h_{1, l}^{2}} S^{2}) + ρ_{2} g h_{2, l} - ρ_{1} g h_{2, l} .

(20)

After some algebraic manipulation, we find that

{F r}_{LE}^{2} = \frac{u_{2, l}^{2}}{g δ h_{2, l}} = \frac{2 {(1 - α)}^{2}}{1 - δ α (2 - α)},

(21)

where α = h_2,l/(h_2,l + h_1,l).

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

{F r}_{LE}^{2} = \frac{2 {(1 - α)}^{2} (1 - δ α / 2)}{1 - 2 δ α (1 - α)} .

(22)

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

The Boussinesq approximation is achieved by setting δ = 0, wherever it is not multiplied by g. In this case, Eqs. (21) and (22) coincide and give that

{F r}_{LE} = \sqrt{2} (1 - α) .

(23)

Figure 4 compares our results from the 2LSWE, Eq. (23), to the model by Benjamin,¹⁸ 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 the vicinity of shocks is not so severe as one would think and as has been repeatedly assumed in the literature.^2,21,33,34

FIG. 4.

View large Download slide

Froude numbers calculated from the 2LSWE in the Boussinesq case [Eq. (23)] compared to the equation by Benjamin¹⁸ [Eq. (2)].

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,^60,61 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⁶² found that 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 Fr_LE = 0.90 ± 0.05 for α = 0.37. Inserting α = 0.37 into Eq. (23), we get Fr_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⁵⁶ 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, η ≔ h₁ − h₂ and ϑ ≔ u₁ − u₂, and used η and ηϑ 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 behavior 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.

A. 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). Furthermore, we let $\bar{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 ∈ {1, 2} in a two-layer shallow-water system is source-bounded if there exists

K \in R

such that the source terms satisfy ∀h_i > K,

\frac{\partial}{\partial h_{i}} |\frac{G_{h_{j}}}{ρ_{i} h_{i}}| < 0 and \frac{\partial}{\partial h_{i}} |\frac{G_{h_{j} u_{j}}}{ρ_{i} h_{i}}| < 0

for j = 1 and j = 2.

Lemma 1.

Let (s, d) = (1, 2) or (2, 1),

{D_{k}}_{k \in N}

be a sequence of increasing real numbers, h₀ and f be scalar functions, and q_1,0 and q_2,0 be vector functions. Furthermore, consider a 2LSWE system with initial conditions,

\begin{aligned} h_{d k} (0, x) & = D_{k} + f (x), \\ h_{s k} (0, x) & = h_{0} (x), \\ q_{1 k} (0, x) & = q_{1,0} (x), \\ q_{2 k} (0, x) & = q_{2,0} (x), \end{aligned}

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_{1 k}, h_{2 k}, q_{1 k}, q_{2 k})}_{k \in N}

be physical solutions to the 2LSWE system. If

{(h_{s k}, h_{d k} - D_{k}, q_{s k}, q_{d k} / D_{k})}_{k \in N}

converge and the first and second derivatives are uniformly bounded in the regions where they are well-defined, then

\lim_{k \to \infty} \frac{h_{d k} (t, x)}{D_{k}} = 1 .

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.⁶³ This implies that we can interchange the order of limits and differentiation.⁶⁴ 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 $\bar{h_{i}} = 1$ for i ∈ {1, 2}, the energy is bounded from below by the height and momentum distributions that give E = gρ_i/2. We let $\tilde{E_{k}} = 2 E_{k} / D_{k}^{2}$ be a scaled energy, and it follows by insertion that $\tilde{E_{k}} (0, x) \to g ρ_{d}$ as k → ∞. That is, the scaled energy $\tilde{E_{k}} (0, x)$ approaches the minimal for a system with $\bar{h_{d k} (t, x) ∕ D_{k}} = 1$ in the limit when k → ∞. Source-boundedness of the mass source terms implies that $\bar{h_{d k} ∕ D_{k}} \to 1$ for all $t \in R$ as k → ∞.

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

\frac{\partial \tilde{E_{k}} (0, x)}{\partial t} \leq 0

in the limit when k → ∞. Because all the terms in Eq. are non-negative and because the right-hand side of the scaled version [Eq. ] remains 0 as long as the scaled energy remains minimal, we must have

\lim_{k \to \infty} |\tilde{E_{k}} (t, x) - \tilde{E_{k}} (0, x)| = 0 .

(24)

Assume that ∃ε > 0 and

\forall N \in N

⁠, ∃k > N, such that

|\frac{h_{d k} (t, x)}{D_{k}} - 1| > ε .

This implies that h_dk(t, x)/D_k deviates from 1 by a term which does not vanish in the limit k → ∞. This contradicts Eq. , as discussed above, so h_dk(t, x)/D_k → 1 for all

t \in R

and

x \in R^{2}

as k → ∞.

□

B. 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 → δg with δ 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 ∈ {1, 2} to indicate which layer is shallow relative to the other such that

\frac{\partial}{\partial t} ρ_{s} h_{s} + \nabla \cdot (ρ_{s} h_{s} u_{s}) = G_{h_{s}},

(25a)

\frac{\partial}{\partial t} u_{s} + (u_{s} \cdot \nabla) u_{s} + δ g \nabla (h_{s} + b^{s - 1}) = \frac{1}{ρ_{s} h_{s}} (G_{h_{s} u_{s}} - u_{s} G_{h_{s}}) .

(25b)

Theorem 1.

Let (s, d) = (1, 2) or (2, 1),

{D_{k}}_{k \in N}

be a sequence of increasing real numbers, h₀ and f be scalar functions, and q_1,0 and q_2,0 be vector functions, all defined on

Ω \subseteq R^{n}

. Furthermore, consider 2LSWE in the form of Eq. or Eq. with initial conditions

\begin{aligned} h_{d k} (0, x) & = D_{k} + f (x), \\ h_{s k} (0, x) & = h_{0} (x), \\ q_{1 k} (0, x) & = q_{1,0} (x), \\ q_{2 k} (0, x) & = q_{2,0} (x), \end{aligned}

in which layer d is source-bounded and the density of layer d is constant. Now, let

{(h_{1 k}, h_{2 k}, q_{1 k}, q_{2 k})}_{k \in N}

be physical solutions to the 2LSWE such that q_dk satisfies the boundary condition,

|q_{d k} (t, x)| \leq K f o r x \in \partial Ω,

(26)

with

K \in R

independent of k and where ∂Ω may be at infinity.

If ${(h_{s k}, h_{d k} - D_{k}, q_{s k}, q_{d k} / D_{k})}_{k \in N}$ converge and the first and second derivatives are uniformly bounded in the regions where they are well-defined, then (h_s, u_sk) → (h, u) where (h, u) solves Eq. (25) in the weak sense, u_dk → 0, and $h_{d k} - D_{k} \to (C - ρ_{s}^{d - 1} h_{s} - ρ_{2}^{d - 1} b) / ρ_{d}^{d - 1}$ , where C is constant in space. If the domain on which the solution is defined is infinite in the range or the mass source terms are zero, then C is equal to $C = {\bar{[ρ_{s}^{d - 1} h_{s} + ρ_{d}^{d - 1} f + ρ_{2}^{d - 1} b]}}_{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.⁵⁰

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 u_dk → 0 in the limit k → ∞ by the use of the fundamental theorem of geometric calculus. The reader is referred to Doran and Lasenby⁶⁵ 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

\nabla \cdot (\frac{q_{d k}}{D_{k}}) = \frac{G_{h_{d}}}{D_{k}} - \frac{\partial}{\partial t} \frac{ρ_{d k} h_{d k}}{D_{k}} \overset{k \to \infty}{\to} 0 \Rightarrow \nabla \cdot u_{d k} \overset{k \to \infty}{\to} 0 .

(27)

Next, we show that also the curl of u_dk vanishes in the limit k → ∞. In the following, we use A ∧ B to denote the wedge product, or exterior product, of A and B and A ⊙ B to denote their geometric product. Applying ∇∧, a generalized curl, from the left of the velocity equation of Eq. yields

\frac{\partial w_{d k}}{\partial t} + \nabla \land I^{- 1} ⊙ u_{d k} \land I^{- 1} ⊙ w_{d k} = \frac{\nabla \land (G_{h_{d} u_{d}} - u_{d k} G_{h_{d}})}{ρ_{d k} h_{d k}},

(28)

where

w

_dk ≡ ∇∧u_dk is a bivector which is equal in magnitude to the curl of u_dk but well-defined in any dimension. Here, I⁻¹ = e_m ∧⋯∧e₁, where 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 imply that

w

_dk → 0 in the limit k → ∞.

From the Helmholtz theorem, we know that a vector field defined on a finite domain or which goes sufficiently fast to 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⁶⁵

\oint_{\partial V} L (I_{m - 1} (x^{'})) d^{m - 1} x^{'} \int_{V} \dot{L} (\dot{\nabla} ⊙ I_{m}) d^{m} x^{'},

(29)

where L is any linear function, I_m is the pseudoscalar of the tangent space to V, and I_m−1(x) is the pseudoscalar of the tangent space to ∂V at x. The vector derivative is ∇ ⊙ A = ∇ · A + ∇ ∧ A, and the overdot indicates where it acts. That is, the integrand on the right-hand side of Eq. is ∑_i∂_iL(e_i ⊙ I_m). See, for instance, the textbook by Doran and Lasenby.⁶⁵

Let V be some region where u_dk is differentiable and let

L (A) = G ⊙ A ⊙ u_{d k} = \frac{x^{'} - x}{S_{m - 1} {|x^{'} - x|}^{m}} ⊙ A ⊙ u_{d k},

(30)

where G is the Green’s function for the vector derivative, meaning that ∇ ⊙ G(x′, x) = δ(x − x′) and S_m−1 is the volume of the (m − 1)-sphere. Equation then states that

\begin{matrix} u_{d k} (x) & = & \frac{I_{m}^{- 1}}{S_{m - 1}} ⊙ \{) {(- 1)}^{m} \int_{V} \frac{x^{'} - x}{{|x^{'} - x|}^{m}} ⊙ I_{m} ⊙ [\nabla \cdot u_{d k} (x^{'})) \\ + (w_{d k} (x^{'})] d^{m} x^{'} + \oint_{\partial V} \frac{x^{'} - x}{{|x^{'} - x|}^{m}} ⊙ I_{m - 1} (x^{'}) \\ ⊙ u_{d k} (x^{'}) d^{m - 1} x^{'} \}) . \end{matrix}

(31)

The surface integral can be decomposed into one vector component whose integrand is proportional to $u_{d k} \cdot \hat{n}$ and one triplet-vector component whose integrand is proportional to $u \land \hat{n}$ ⁠. Only the vector component will contribute to u_dk. To show that u_dk → 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} \hat{n} \cdot u_{d k}$ is continuous. By applying Eq. (29) with L given by Eq. (30) on a domain which does not include x, the only contribution comes from surface integral in the limit k → ∞ because lim_k→∞∇ ⋅ u_dk + $w$ _dk = 0 everywhere. I_m−1 has opposite sign on opposite sides on surfaces, so surface integrals from neighboring domains cancel as lim_k→∞u_dk is continuous. Thus, we can extend the integral over ∂V to an integral over ∂Ω by applying the fundamental theorem of geometric calculus in the neighboring domains. From Eq. (26) with Lemma 1, we get that lim_k→∞u_dk must vanish on ∂Ω. Hence, lim_k→∞u_dk = 0 everywhere.

From the momentum equations in Eq. , then

\begin{matrix} \nabla [ρ_{1 k}^{d - 1} h_{1 k} + ρ_{2 k}^{d - 1} (h_{2 k} + b)] & = & () \frac{G_{h_{d} u_{d}} - u_{d k} G_{h_{d}}}{ρ_{d k} h_{d k}} - \frac{\partial u_{d k}}{\partial t} \\ - (u_{d k} \cdot \nabla) u_{d k})) ρ_{d k}^{d - 1} \to 0, \end{matrix}

(32)

and so in the limit k → ∞,

ρ_{s}^{d - 1} h_{s} + ρ_{d}^{d - 1} (h_{d k} - D_{k}) + ρ_{2}^{d - 1} b

is constant in space. Finally, plugging this into the equation for u_s in Eq. or Eq. , we get Eq. . In the regions where the solution is differentiable, the various formulations of the 2LSWE, Eqs. , 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 the Hugoniot locus of Eq. . Let

γ : = 1 / ⟨h_{d}⟩

⁠. In Appendix B, we show that the full set of Rankine-Hugoniot conditions for the 2LSWE [Eqs. and ] may be written as

S 〚 ρ_{s} h_{s} 〛 = \hat{n} \cdot 〚 ρ_{s} h_{s} u_{s} 〛,

(33a)

S \hat{n} \cdot 〚 u_{s} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{s})}^{2} + δ g h_{s} 〛 + g_{1} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}),

(33b)

〚 ρ_{1}^{d - 1} h_{1} + ρ_{2}^{d - 1} h_{2} 〛 = g_{2} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}),

(33c)

S \hat{n} \cdot 〚 u_{d} 〛 = g_{3} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}),

(33d)

where

g_{3} = γ S 〚 h_{d} 〛 (S - ⟨\hat{n} \cdot u_{d}⟩) .

(34)

For Eq. ,

g_{1} = γ 〚 h_{d} 〛 {(S - ⟨\hat{n} \cdot u_{d}⟩)}^{2},

(35)

g_{2} = \frac{ρ_{2}^{d - 1}}{g} g_{1} .

(36)

For Eq. with d = 1,

g_{1} = (γ 〚 h_{1} 〛 {(S - ⟨\hat{n} \cdot u_{1}⟩)}^{2} + δ g g_{2}) \hat{n}

(37)

and

\begin{matrix} g_{2} & = & \frac{γ}{g} () \frac{S}{ρ_{1}} 〚 ρ_{2} h_{2} \hat{n} \cdot u_{2} 〛 - \frac{1}{ρ_{1}} 〚 ρ_{2} h_{2} {(\hat{n} \cdot u_{2})}^{2} 〛 \\ + 〚 h_{1} 〛 (S ⟨\hat{n} \cdot u_{1}⟩ - ⟨{(\hat{n} \cdot u_{1})}^{2}⟩) - g ⟨h_{2}⟩ 〚 h_{1} 〛 - \frac{ρ_{2}}{2 ρ_{1}} 〚 h_{2}^{2} 〛 \\ + 〚 h_{1} 〛 (S - ⟨\hat{n} \cdot u_{1}⟩) (S - 2 ⟨\hat{n} \cdot u_{1}⟩))) . \end{matrix}

(38)

Finally, for Eq. with d = 2,

g_{1} = γ 〚 h_{2} 〛 {(S - ⟨\hat{n} \cdot u_{2}⟩)}^{2}

(39)

and

\begin{matrix} g_{2} & = & \frac{γ}{g} () S 〚 ρ_{1} h_{1} \hat{n} \cdot u_{1} 〛 - 〚 ρ_{1} h_{1} {(\hat{n} \cdot u_{1})}^{2} 〛 + S ρ_{2} 〚 h_{2} 〛 ⟨\hat{n} \cdot u_{2}⟩ \\ - ρ_{2} 〚 h_{2} 〛 ⟨{(\hat{n} \cdot u_{2})}^{2}⟩ - 〚 \frac{1}{2} g ρ_{1} h_{1}^{2} 〛 - ρ_{1} g ⟨h_{1}⟩ 〚 h_{2} 〛 \\ + ρ_{2} 〚 h_{2} 〛 (S - ⟨\hat{n} \cdot u_{2}⟩) (S - 2 ⟨\hat{n} \cdot u_{2}⟩))) . \end{matrix}

(40)

In particular, we note that g₁, g₂, and g₃ vanish for γ = 0 in all cases.

Next, we notice that Eq. (33) with γ = 0 is exactly the Rankine-Hugoniot relations for the locally conservative 1LSWE (25) together with the conditions that $ρ_{1}^{d - 1} h_{1} + ρ_{2}^{d - 1} h_{2}$ is constant and u_d = 0. From Lemma 1, it follows that lim_k→∞γ = 0. Thus, the Hugoniot loci match and this concludes the proof of the theorem. $□$

C. Discussion of the theorem

Theorem 1 shows that we may approximate the thinnest layer of the 2LSWE with the locally conservative 1LSWE where g → (1 − ρ₁/ρ₂)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.

FIG. 5.

View large Download slide

A sketch of how the two layers converge to one-layer cases with increasing D_k for both cases (s, d) = (1, 2) and (s, d) = (2, 1).

The interesting part about Theorem 1 is not that smooth solutions of the 2LSWE 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 have 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}]$ ⁠.^{2,5,18,19,21,33} 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 ${F r}_{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 Fr_LE, have had to track the position of the leading edge and alter the solution.^33,66 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 ${F r}_{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 Ref. 67. Adding viscosity smooths out discontinuities and renders the interaction term h₁∇h₂ 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 traveling 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.⁴⁵ In the two-layer dam-break problem, a lighter fluid of height h₁ spreads on top of a heavier fluid of height h₂ as shown in Fig. 6. The problem has been frequently used in the literature as a benchmark case for spreading models.^68–70

FIG. 6.

View large Download slide

A simple sketch of the one-dimensional dam-break problem.

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 Fig. 7.

FIG. 7.

View large Download slide

A tabular overview of initial conditions for the various test cases. The cases Ra, Rb, and Rc all have reflecting walls to the left and differ in the initial configuration of the fluids. Note that case Rc uses the same initial conditions as Ra with s = 1.

A. 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

\begin{aligned} h_{1} (t = 0, x) & = \{\begin{matrix} h_{0} & if x \leq 0, \\ 0 & if x > 0, \end{matrix}) \\ h_{2} (t = 0, x) & = H - (1 - δ) h_{1}, \\ u_{1} (t = 0, x) & = 0, \\ u_{2} (t = 0, x) & = 0, \end{aligned}

(41)

where h₀ 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,⁷¹

h (x, t) = \{\begin{matrix} h_{0} & if x \leq - c_{0} t, \\ \frac{h_{0}}{9} {(2 - \frac{x}{c_{0} t})}^{2} & if - c_{0} t < x \leq 2 c_{0} t, \\ 0 & if 2 c_{0} t < x, \end{matrix})

(42a)

u (x, t) = \{\begin{matrix} 0 & if x \leq - c_{0} t, \\ \frac{2}{3} (c_{0} + \frac{x}{t}) & if - c_{0} t < x \leq 2 c_{0} t, \\ 0 & if 2 c_{0} t < x, \end{matrix})

(42b)

where $c_{0} = \sqrt{δ g h_{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

h (x, t) = \{\begin{matrix} h_{0} & if x \leq - c_{0} t, \\ \frac{h_{0}}{9} {(2 - \frac{x}{c_{0} t})}^{2} & if - c_{0} t < x \leq \frac{(2 - \sqrt{2}) c_{0} t}{1 + \sqrt{2}}, \\ \frac{4 h_{0}}{{(2 + \sqrt{2})}^{2}} & if \frac{(2 - \sqrt{2}) c_{0} t}{1 + \sqrt{2}} < x \leq \frac{2 c_{0} t}{1 + \sqrt{2}}, \\ 0 & if \frac{2 c_{0} t}{(1 + \sqrt{2})} < x, \end{matrix})

(43a)

u (x, t) = \{\begin{matrix} 0 & if x \leq - c_{0} t, \\ \frac{2}{3} (c_{0} + \frac{x}{t}) & if - c_{0} t < x \leq \frac{(2 - \sqrt{2}) c_{0} t}{1 + \sqrt{2}}, \\ \frac{2 c_{0}}{1 + \sqrt{2}} & if \frac{(2 - \sqrt{2}) c_{0} t}{1 + \sqrt{2}} < x \leq \frac{2 c_{0} t}{1 + \sqrt{2}}, \\ 0 & if \frac{2 c_{0} t}{1 + \sqrt{2}} < x . \end{matrix})

(43b)

A sketch of the two solutions for h is shown in Fig. 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 ${F r}_{LE} = \sqrt{2}$ at the leading edge.

FIG. 8.

View large Download slide

A sketch of the Ritter solution (42) and Eq. (43) for the one-layer dam-break problem.

B. 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 (∂_xh)(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

\begin{aligned} h_{1} (t = 0, x) & = H - h_{2}, \\ h_{2} (t = 0, x) & = \{\begin{matrix} h_{0} & if x \leq 0, \\ 0 & if x > 0, \end{matrix}) \\ u_{1} (t = 0, x) & = 0, \\ u_{2} (t = 0, x) & = 0, \end{aligned}

(44)

where again h₀ is constant.

C. Case Rb: Effect of nonzero depth on both sides of dam

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

\begin{aligned} h_{1} (t = 0, x) & = \{\begin{matrix} h_{0, a} & if x \leq 0, \\ h_{0, b} & if x > 0, \end{matrix}) \\ h_{2} (t = 0, x) & = H - (1 - δ) h_{1}, \\ u_{1} (t = 0, x) & = 0, \\ u_{2} (t = 0, x) & = 0 . \end{aligned}

(45)

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 capture the quantitative behavior of two-layer cases that is not captured by the globally conservative 1LSWE.

D. 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,⁷² who studied the spreading of oil on water. He used a 2.5 m long and 0.62 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 I). The initial depth of water in the experiments was about 30 cm, which is nearly twice the initial heights.

TABLE I.

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₀ (cm)	Width L (cm)	δ
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	2 L methane	2.0	17.3	7	0.746
Chang, Reid, and Fay	2 L nitrogen	2.0	17.3	7	0.34
Chang, Reid, and Fay	0.75 L methane	0.75	6.5	7	0.746
Chang, Reid, and Fay	1 L nitrogen	1.0	8.7	7	0.34

Authors	Experiment	Spill volume (L)	Height h₀ (cm)	Width L (cm)	δ
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	2 L methane	2.0	17.3	7	0.746
Chang, Reid, and Fay	2 L nitrogen	2.0	17.3	7	0.34
Chang, Reid, and Fay	0.75 L methane	0.75	6.5	7	0.746
Chang, Reid, and Fay	1 L nitrogen	1.0	8.7	7	0.34

The second set of experiments that will be considered are those presented by Chang, Reid, and Fay.⁷³ They studied fluids at cryogenic temperatures (cryogens) spreading on water and presented both experimental results and model predictions. In their model, they used the same empirical boundary condition for the spreading rate as discussed in Sec. III.⁷⁴ 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⁷³ deviates from the dam-break case in that the initial reservoir of the cryogen is emptied through a large slit. The spreading then occurred inside a cylinder of length 4 m with an inside diameter of 16.5 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.⁷⁴ 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 I.

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⁷³ and include constant evaporation rates of 0.16 kg m⁻² s⁻¹ for methane and 0.201 kg m⁻² s⁻¹ 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.⁷⁵ 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 kg m⁻³ for nitrogen and 254 kg m⁻³ for methane. Similar to Chang, Reid, and Fay, we will use reduced densities in our simulations as well. Finally, Chang, Reid, and Fay reported 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 Sec. V. The equations are discretized spatially using a finite-volume scheme. We employ the FORCE (first-order centered) flux⁷⁶ and the second-order MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) reconstruction with a minmod limiter⁴⁵ in each finite volume. The solutions are advanced in time with a standard third-order three-stage strong stability-preserving Runge-Kutte method.⁷⁷ 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¹ norm, which for a function $y : Ω \to R$ is defined as

{‖y‖}_{1} \equiv \int_{Ω} |y (x)| d^{n} x .

(46)

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

A. 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₀ = 1 m and δ = 0.7 with 400 grid cells on the domain x/c₀t ∈ (−2, 2). The initial depth, H, was increased stepwise from the initial value, H = h₀, to obtain a larger height difference between the layers. The results presented in Fig. 9 show that the locally conservative 2LSWE converge toward the analytic solution of the locally conservative 1LSWE when H increases. A higher value of H translates into increasing D_k in Theorem 1. This 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.

FIG. 9.

View large Download slide

Comparison of solutions of the local 2LSWE with increasing depths, H (blue lines) to the analytic solution of the locally conservative 1LSWE for case 0. As H increases, the 2LSWE solution approaches the 1LSWE solution.

Figure 10 shows the normalized difference in L¹ for the top-layer height between the analytic solution [Eq. (43)] and the solutions obtained with the 2LSWE, ${\tilde{h}}_{1}$ and h₁, respectively (⁠ $∥ h_{1} - {\tilde{h}}_{1} ∥_{1} / {‖h_{0}‖}_{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₂, u₁, and u₂, were found to exhibit a similar behavior.

FIG. 10.

View large Download slide

The L¹ difference of the top-layer height, h₁, between analytic solutions of case 0 and solutions with the 2LSWE for different initial depths, H. The circles correspond to the locally conservative 2LSWE, and the triangles correspond to the globally conservative 2LSWE. The line indicates a slope of −1.

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

In case Ra, a reflective wall is placed at x = −L (cf. Sec. V B). The case was solved with the parameters h₀ = 1 m, L = 2 m, and δ = 0.6. The domain width was 15 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 s. An illustration of the initial configurations is presented at the top of this 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₂ 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.

FIG. 11.

View large Download slide

Height distribution of the two layers in a dam-break problem at t = 3 s solved with the locally conserved 1LSWE (25) (blue lines), the globally conserved 1LSWE (4) (red lines), and the globally conserved 2LSWE (8) (green lines). The solid lines and dashed lines indicate h₁ + h₂ and h₂, respectively.

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 Fig. 12(a), the L¹ difference of the top-layer height, $∥ h_{1}^{2 LSWE} - h_{1}^{1 LSWE} ∥_{1} / h_{0} L$ ⁠, is plotted for varying times and depths, H/h₀. Figure 12(b) shows the difference of the leading edge position at t = 5 s, |r_2LSWE − r_1LSWE|/h₀, for varying initial depths, H/h₀, and widths, L/h₀. The position of the leading edge is here defined as the smallest x-value where the top layer is thinner than 10⁻⁴ m. Figure 12 shows that the differences in h₁ decrease with time, which is reasonable since the spreading fluid becomes gradually thinner. As expected, the differences decrease with increasing value of H/h₀. Similar to case 0, the errors in the variables h₂, u₁, and u₂ as quantified by the L¹ norm exhibit the same trends as the top layer height (not shown). Furthermore, this figure also 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.

FIG. 12.

View large Download slide

A quantitative comparison of solutions from the globally conservative 2LSWE and the locally conservative 1LSWE for case Ra, showing (a) the L¹ difference in top-layer height and (b) the difference in leading edge position at t = 5 s.

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. In addition, for this case, we observe that the spreading rate from the globally conservative 1LSWE is higher than that from the full 2LSWE.

FIG. 13.

View large Download slide

The position of leading edge as a function of time for case Ra. The blue lines show the result for 2LSWE with various depths.

C. Case Rb: Effect of nonzero depth at both sides of dam

In case Rb, the dam break was initialized according to Eq. (45) with a nonzero depth at both sides of the dam. The case was solved with the parameters δ = 0.6, L = 4 m, H = 50 m, h_0,a = 2 m, and h_0,b = 0.5 m. The width of the domain is 15 m, and the results are again computed with 1000 grid cells.

Figure 14 shows the height distributions at time t = 3 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. This figure shows that the locally conservative 1LSWE should be used for accurate representation of the position of the leading edge.

FIG. 14.

View large Download slide

Height profiles of h₂ and h₁ + h₂ at t = 3 s of solutions to case Rb with the different formulations of the 1LSWE and 2LSWE. “Loc” and “Glob” denote the locally conservative and globally conservative formulation, respectively.

D. 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 Sec. V D (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 β = 1.31 is also included for the Suchon experiments (orange dotted lines).

FIG. 15.

View large Download slide

The spreading distance of oil on water as a function of time. A comparison of results from the locally conservative 1LSWE, the globally conservative 2LSWE, the Fay model, and experimental data from Suchon.⁷²

FIG. 16.

View large Download slide

The spreading distance of liquid methane (CH₄) and nitrogen (N) on water as a function of time. A comparison of results from the locally conservative 1LSWE and experiments by Chang, Reid, and Fay.⁷³

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¹ norm, that is,

d e v (r_{\exp}, r_{sim}) = \frac{‖ r_{\exp} - r_{sim} ‖}{‖ r_{\exp} ‖} .

The deviation in the spreading radius is 4.5% for oil on water, which is significantly better than the Fay model, where the average deviations are 12.6%. For the cryogenic fluids, the deviations in the spreading radius were 10.2% for methane on water and 4.2% 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 is shown in Fig. 12(b) that the difference between the predicted spreading distance from the 2LSWE and the one-layer approximation after 5 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 s that is smaller than 0.2 times the initial height h₀.

A comparison of the green solid lines (1LSWE) and the blue dashed lines (2LSWE) in Fig. 15 confirms a very good agreement between the two formulations. For Fig. 16, and especially for the 2L cases where the initial height is large compared to the water depth, the discrepancy is larger. We find that the deviation in the L¹ norm between the one-layer approximation and 2LSWE results is between 0.7% and 5.5% for all cases. In the 2L cryogen experiments, we observe that the discrepancy is reducing after some time. This is consistent with Fig. 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. However, 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 Fig. 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 is marked by red dots. We see that the 2LSWE capture a more complex behavior, especially in the early phase of the flow, but as expected, the agreement between the profiles improves with time.

FIG. 17.

View large Download slide

A comparison of height profiles produced by the one-layer approximation and the globally conservative 2LSWE at different times for the 2L liquid methane case.

Finally, we note that Chang, Reid, and Fay⁷³ also solved the 1LSWE numerically, but with an imposed boundary condition with Fr_LE = 1.28 at the leading edge. They motivated the use of a constant Froude number at the leading edge by frequent use in the previous literature dealing with spreading of nonboiling fluids. They treated the value of Fr_LE as a parameter that depends on the apparatus and must be determined experimentally and found that Fr_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, Fr_LE approaches to $\sqrt{2}$ from below. Moreover, our analysis shows that ${F r}_{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 that of 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 ${F r}_{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 − ρ₁/ρ₂)g, the one-layer model correctly captures the behavior of shocks and contact discontinuities. In particular, the one-layer model results in ${F r}_{LE} = \sqrt{2}$ ⁠, which is exactly the same as the theoretical predictions by von Kármán,¹⁷ Benjamin,¹⁸ and Ungarish¹⁹ 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 the 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% for oil on water, 10.2% for methane on water, and 4.2% for nitrogen on water. The spreading radius from the one-layer and two-layer descriptions could hardly be distinguished from each other after 10 s 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% 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 Sec. 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 modeling of geophysical phenomena, they represent an attractive possibility for future work.

ACKNOWLEDGMENTS

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 No. 244076/O80).

APPENDIX A: DERIVING RANKINE-HUGONIOT CONDITIONS

To obtain the Rankine-Hugoniot condition for the 2LSWE, we use an approach similar to that of Smoller.⁷⁸ Consider a shock along Γ which is normal to $\hat{n}$ and let ϕ be a test function with compact support D, which lies in the x_nt-plane.

Consider first the locally conservative system. Because the velocity u_i is a weak solution, the normal component must satisfy

\begin{matrix} \int_{D} () ϕ_{t} \hat{n} \cdot u_{i} + ϕ_{n} [) \frac{1}{2} {(\hat{n} \cdot u_{i})}^{2} + {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} g h_{1} \\ + g h_{2} + g b]) - ϕ J_{i})) d x_{n} d t = 0, \end{matrix}

(A1)

where the subscript on ϕ denotes partial differentiations and

J_{i} = u_{i}^{T} \partial_{T} (\hat{n} \cdot u_{i}) - \frac{\hat{n} \cdot (G_{h_{i} u_{i}} - u_{i} G_{h_{i}})}{ρ_{i} h_{i}},

(A2)

where $u_{i}^{T}$ is the tangential component of u_i and ∂_T is differentiation with respect to the tangential direction. The integrand in Eq. (A1) 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₁ and D₂. Because the solution is differentiable inside these regions, we may use Green’s theorem and obtain

\begin{matrix} \int_{D_{j}} (ϕ_{t} \hat{n} \cdot u_{i} + ϕ_{n} [\frac{1}{2} {(\hat{n} \cdot u_{i})}^{2} + {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} g h_{1} + g h_{2} + g b] - ϕ J_{i}) d x_{n} d t \\ = \int_{D_{j}} (\frac{\partial}{\partial t} (ϕ \hat{n} \cdot u_{i}) + \frac{\partial}{\partial x_{n}} [ϕ (\frac{1}{2} {(\hat{n} \cdot u_{i})}^{2}))) \\ (((+ {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} g h_{1} + g h_{2} + g b)]) d x_{n} d t \\ = \pm \lim_{ε \to 0} \int_{Γ \pm ε x_{n}} ϕ ([\frac{1}{2} {(\hat{n} \cdot u_{i})}^{2} + {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} g h_{1})) \\ ((+ g h_{2} + g b] d t - \hat{n} \cdot u_{i} d x_{n}) \end{matrix}

(A3)

because ϕ vanish on the boundary of D. Equation (A1) is obtained by adding Eq. (A3) with j = 1 and j = 2. Because Eq. (A1) must hold for all test functions, we obtain the Rankine-Hugoniot condition for the normal velocity component,

S 〚 \hat{n} \cdot u_{i} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{i})}^{2} + g {(\frac{ρ_{1}}{ρ_{2}})}^{i - 1} h_{1} + g h_{2} 〛 .

(A4)

For the globally conservative system, a similar treatment yields

\begin{matrix} S 〚 ρ_{1} h_{1} u_{1} + ρ_{2} h_{2} u_{2} 〛 & = & 〚 (\hat{n} \cdot u_{1}) ρ_{1} h_{1} u_{1} + (\hat{n} \cdot u_{2}) ρ_{2} h_{2} u_{2} 〛 \\ + 〚 \frac{1}{2} g ρ_{1} h_{1}^{2} + ρ_{1} g h_{1} h_{2} + \frac{1}{2} ρ_{2} g h_{2}^{2} 〛 \hat{n} \end{matrix}

(A5)

and

S \hat{n} \cdot 〚 u_{2} - u_{1} 〛 = 〚 \frac{1}{2} [) {(\hat{n} \cdot u_{2})}^{2} - {(\hat{n} \cdot u_{1})}^{2}]) - g δ h_{1} 〛 .

(A6)

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

Ostapenko⁴⁸ 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 $〚 u_{i}^{T} 〛 = 0$ ⁠. A conservation law for the quantity ∂₁u_i,2 − ∂₂u_i,1 ≡ $w$ _i can be obtained by taking distributional derivatives of the different components of the velocity equation, yielding

\frac{\partial w_{i}}{\partial t} + \nabla \cdot (w_{i} u_{i} + J_{i}^{⊥}) = 0,

(A7)

where

J_{i}^{⊥} = \frac{1}{ρ_{i} h_{i}} (\begin{matrix} G_{h_{i} u_{i}, 2} - u_{i, 2} G_{h_{i}} \\ - G_{h_{i} u_{i}, 1} + u_{i, 1} G_{h_{i}} \end{matrix}) .

(A8)

It may be tempting from Eq. (A7) to conclude that the vorticity must obey the jump condition,

S 〚 w_{i} 〛 = 〚 w_{i} \hat{n} \cdot u_{i} + \hat{n} \cdot J_{i}^{⊥} 〛 .

(A9)

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 cannot 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,

S 〚 ρ_{s} h_{s} 〛 = \hat{n} \cdot 〚 ρ_{s} h_{s} u_{s} 〛,

(B1a)

S \hat{n} \cdot 〚 u_{s} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{s})}^{2} + δ g h_{s} 〛 + g_{1} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}),

(B1b)

〚 ρ_{1}^{d - 1} h_{1} + ρ_{2}^{d - 1} h_{2} 〛 = g_{2} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}),

(B1c)

S \hat{n} \cdot 〚 u_{d} 〛 = g_{3} (γ, S, h_{s}, \hat{n} \cdot u_{s}, \hat{n} \cdot u_{d}) .

(B1d)

Here, g₁ and g₂ differ for Eqs. (7) and (8), while g₃ will be the same. Furthermore, we will show that all of g₁, g₂, and g₃ vanish when γ = 0. Note that Eq. (B1a) follows directly from Eq. (11) applied to the shallowest layer.

We first consider g₃, which can be obtained from mass conservation of layer d. The scalar Rankine-Hugoniot condition (11) immediately yields

S 〚 h_{d} 〛 = \hat{n} \cdot 〚 h_{d} u_{d} 〛 \Rightarrow \hat{n} \cdot 〚 u_{d} 〛 = \frac{〚 h_{d} 〛}{⟨h_{d}⟩} (S - ⟨\hat{n} \cdot u_{d}⟩) = \frac{g_{3}}{S},

where we used that $〚 a b 〛 = 〚 a 〛 ⟨b⟩ + ⟨a⟩ 〚 b 〛$ ⁠. This gives

g_{3} = γ S 〚 h_{d} 〛 (S - ⟨\hat{n} \cdot u_{d}⟩) .

(B2)

Next, we consider the expressions for g₁ and g₂. We first consider the locally conservative momentum equations [(7c) and (7d)]. We apply the Rankine-Hugoniot condition (14) with i = d and insert Eq. (B2) to obtain

〚 ρ_{1}^{d - 1} h_{1} + ρ_{2}^{d - 1} h_{2} 〛 = \frac{ρ_{2}^{d - 1} 〚 h_{d} 〛}{g ⟨h_{d}⟩} {(S - ⟨\hat{n} \cdot u_{d}⟩)}^{2},

(B3)

that is,

g_{2} = \frac{ρ_{2}^{d - 1}}{g} γ 〚 h_{d} 〛 {(S - ⟨\hat{n} \cdot u_{d}⟩)}^{2} .

(B4)

Now consider Eq. (14) with i = s,

S \hat{n} \cdot 〚 u_{s} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{s})}^{2} + g {(\frac{ρ_{1}}{ρ_{2}})}^{s - 1} h_{1} + g h_{2} 〛 .

(B5)

If we consider the cases s = 1 and s = 2 separately and use that d − 1 = 2 − s, we find from Eq. (B3) that

〚 g {(\frac{ρ_{1}}{ρ_{2}})}^{s - 1} h_{1} + g h_{2} 〛 = 〚 δ g h_{s} 〛 + {(\frac{ρ_{1}}{ρ_{2}})}^{s - 1} γ 〚 h_{d} 〛 {(S - ⟨\hat{n} \cdot u_{d}⟩)}^{2},

(B6)

where δ = 1 − ρ₁/ρ₂, as defined in Eq. (9). Inserting into Eq. (B4), we get that

S \hat{n} \cdot 〚 u_{s} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{s})}^{2} + δ g h_{s} 〛 + g_{1},

(B7)

with

g_{1} = {(\frac{ρ_{1}}{ρ_{2}})}^{s - 1} γ 〚 h_{d} 〛 {(S - ⟨\hat{n} \cdot u_{d}⟩)}^{2} .

(B8)

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. (15) and (16)] with $\hat{n}$ and use Eq. (B2), we obtain

\begin{align} 〚 ρ_{1} h_{1} + ρ_{2} h_{2} 〛 & = \frac{1}{g ⟨h_{2}⟩} () S 〚 ρ_{1} h_{1} \hat{n} \cdot u_{1} 〛 - 〚 ρ_{1} h_{1} {(\hat{n} \cdot u_{1})}^{2} 〛 \\ + S ρ_{2} 〚 h_{2} 〛 ⟨\hat{n} \cdot u_{2}⟩ - ρ_{2} 〚 h_{2} 〛 ⟨{(\hat{n} \cdot u_{2})}^{2}⟩ \\ - 〚 \frac{1}{2} g ρ_{1} h_{1}^{2} 〛 - ρ_{1} g ⟨h_{1}⟩ 〚 h_{2} 〛 + ρ_{2} 〚 h_{2} 〛 \\ \times (S - ⟨\hat{n} \cdot u_{2}⟩) (S - 2 ⟨\hat{n} \cdot u_{2}⟩))), \end{align}

(B9a)

S \hat{n} \cdot 〚 u_{1} 〛 = 〚 \frac{1}{2} {(\hat{n} \cdot u_{1})}^{2} + g δ h_{1} 〛 + \frac{〚 h_{2} 〛}{⟨h_{2}⟩} {(S - ⟨\hat{n} \cdot u_{2}⟩)}^{2},

(B9b)

that is,

g_{1} = γ 〚 h_{2} 〛 {(S - ⟨\hat{n} \cdot u_{2}⟩)}^{2}

(B10)

and

\begin{matrix} g_{2} & = & \frac{γ}{g} () S 〚 ρ_{1} h_{1} \hat{n} \cdot u_{1} 〛 - 〚 ρ_{1} h_{1} {(\hat{n} \cdot u_{1})}^{2} 〛 \\ + S ρ_{2} 〚 h_{2} 〛 ⟨\hat{n} \cdot u_{2}⟩ - ρ_{2} 〚 h_{2} 〛 ⟨{(\hat{n} \cdot u_{2})}^{2}⟩ \\ - 〚 \frac{1}{2} g ρ_{1} h_{1}^{2} 〛 - ρ_{1} g ⟨h_{1}⟩ 〚 h_{2} 〛 \\ + ρ_{2} 〚 h_{2} 〛 (S - ⟨\hat{n} \cdot u_{2}⟩) (S - 2 ⟨\hat{n} \cdot u_{2}⟩))) . \end{matrix}

(B11)

Next, we consider the case when d = 1. We may then write Eq. (15) as

\begin{matrix} 〚 h_{1} + h_{2} 〛 & = & g_{2} = \frac{γ}{g} () \frac{S}{ρ_{1}} 〚 ρ_{2} h_{2} \hat{n} \cdot u_{2} 〛 \\ - \frac{1}{ρ_{1}} 〚 ρ_{2} h_{2} {(\hat{n} \cdot u_{2})}^{2} 〛 + 〚 h_{1} 〛 () S ⟨\hat{n} \cdot u_{1}⟩ \\ - ⟨{(\hat{n} \cdot u_{1})}^{2}⟩)) - g ⟨h_{2}⟩ 〚 h_{1} 〛 - \frac{ρ_{2}}{2 ρ_{1}} 〚 h_{2}^{2} 〛 \\ + 〚 h_{1} 〛 (S - ⟨\hat{n} \cdot u_{1}⟩) (S - 2 ⟨\hat{n} \cdot u_{1}⟩))) . \end{matrix}

(B12)

We then insert this into Eq. (16) to get

\begin{matrix} S \hat{n} \cdot 〚 u_{2} 〛 & = & 〚 \frac{1}{2} {(\hat{n} \cdot u_{2})}^{2} + δ g h_{2} 〛 \\ + (γ 〚 h_{1} 〛 {(S - ⟨\hat{n} \cdot u_{1}⟩)}^{2} + δ g g_{2}), \end{matrix}

(B13)

which gives

g_{1} = (γ 〚 h_{1} 〛 {(S - ⟨\hat{n} \cdot u_{1}⟩)}^{2} + δ g g_{2}) .

(B14)

REFERENCES

1.

B.

Franklin

,

W.

Brownrigg

, and

Farish

, “

XLIV. Of the stilling of waves by means of oil. Extracted from sundry letters between Benjamin Franklin, LL. D. F. R. S. William Brownrigg, M. D. F. R. S. and the reverend Mr. Farish

,”

Philos. Trans. R. Soc. London

64

,

445

–

460

(

1774

).

https://doi.org/10.1098/rstl.1774.0044

Google Scholar

Crossref

2.

D. P.

Hoult

, “

Oil spreading on the sea

,”

Annu. Rev. Fluid Mech.

4

,

341

–

368

(

1972

).

https://doi.org/10.1146/annurev.fl.04.010172.002013

Google Scholar

Crossref

3.

J. A.

Fay

, “

Physical processes in the spread of oil on a water surface

,”

Int. Oil Spill Conf. Proc.

1971

,

463

–

467

.

https://doi.org/10.7901/2169-3358-1971-1-463

Crossref

4.

R.

Chebbi

, “

Viscous-gravity spreading of oil on water

,”

AIChE J.

47

,

288

–

294

(

2001

).

https://doi.org/10.1002/aic.690470207

Google Scholar

Crossref

5.

J.

Fay

, “

Spread of large LNG pools on the sea

,”

J. Hazard. Mater.

140

,

541

–

551

(

2007

), part of special issue: LNG Special Issue - Dedicated to Risk Assessment and Consequence Analysis for Liquefied Natural Gas Spills.

https://doi.org/10.1016/j.jhazmat.2006.10.024

Google Scholar

Crossref

PubMed

6.

J.

Fay

, “

Model of spills and fires from LNG and oil tankers

,”

J. Hazard. Mater.

96

,

171

–

188

(

2003

).

https://doi.org/10.1016/s0304-3894(02)00197-8

Google Scholar

Crossref

PubMed

7.

J.

Brandeis

and

D. L.

Ermak

, “

Numerical simulation of liquefied fuel spills: II. Instantaneous and continuous LNG spills on an unconfined water surface

,”

Int. J. Numer. Methods Fluids

3

,

347

–

361

(

1983

).

https://doi.org/10.1002/fld.1650030405

Google Scholar

Crossref

8.

J. F.

Stanislav

,

S.

Kokal

, and

M. K.

Nich

, “

Gas liquid flow in downward and upward inclined pipes

,”

Can. J. Chem. Eng.

64

,

881

–

890

(

1986

).

https://doi.org/10.1002/cjce.5450640601

Google Scholar

Crossref

9.

C.

Adduce

,

G.

Sciortino

, and

S.

Proietti

, “

Gravity currents produced by lock exchanges: Experiments and simulations with a two-layer shallow-water model with entrainment

,”

J. Hydraul. Eng.

138

,

111

–

121

(

2012

).

https://doi.org/10.1061/(asce)hy.1943-7900.0000484

Google Scholar

Crossref

10.

J. O.

Shin

,

S. B.

Dalziel

, and

P. F.

Linden

, “

Gravity currents produced by lock exchange

,”

J. Fluid Mech.

521

,

1

–

34

(

2004

).

https://doi.org/10.1017/s002211200400165x

Google Scholar

Crossref

11.

T.

Moodie

, “

Gravity currents

,”

J. Comput. Appl. Math.

144

,

49

–

83

(

2002

), part of special issue: Selected papers of the International Symposium on Applied Mathematics, August 2000, Dalian, China.

https://doi.org/10.1016/s0377-0427(01)00551-9

Google Scholar

Crossref

12.

C. A.

Mériaux

,

T.

Zemach

,

C. B.

Kurz-Besson

, and

M.

Ungarish

, “

The propagation of particulate gravity currents in a v-shaped triangular cross section channel: Lock-release experiments and shallow-water numerical simulations

,”

Phys. Fluids

28

,

036601

(

2016

).

https://doi.org/10.1063/1.4942241

Google Scholar

Crossref

13.

F.

Stickland

, “

The formation of monomolecular layers by spreading a copper stearate solution

,”

J. Colloid Interface Sci.

40

,

142

–

153

(

1972

).

https://doi.org/10.1016/0021-9797(72)90003-3

Google Scholar

Crossref

14.

A.

Mar

and

S. G.

Mason

, “

Coalescence in three-phase fluid systems

,”

Kolloid-Z. Z. Polym.

224

,

161

–

172

(

1968

).

https://doi.org/10.1007/bf01533973

Google Scholar

Crossref

15.

F. M.

White

,

Fluid Mechanics

, 7th ed. (

McGraw-Hill

,

New York

,

2011

).

Google Scholar

16.

C. L.

Vaughan

and

M. J.

O’Malley

, “

Froude and the contribution of naval architecture to our understanding of bipedal locomotion

,”

Gait Posture

21

,

350

–

362

(

2005

).

https://doi.org/10.1016/j.gaitpost.2004.01.011

Google Scholar

Crossref

PubMed

17.

T.

von Kármán

, “

The engineer grapples with nonlinear problems

,”

Bull. Am. Math. Soc.

46

,

615

–

683

(

1940

).

https://doi.org/10.1090/s0002-9904-1940-07266-0

Google Scholar

Crossref

18.

T. B.

Benjamin

, “

Gravity currents and related phenomena

,”

J. Fluid Mech.

31

,

209

–

248

(

1968

).

https://doi.org/10.1017/s0022112068000133

Google Scholar

Crossref

19.

M.

Ungarish

, “

Benjamin’s gravity current into an ambient fluid with an open surface

,”

J. Fluid Mech.

825

,

R1

(

2017

).

https://doi.org/10.1017/jfm.2017.460

Google Scholar

Crossref

20.

J. A.

Fay

, “

The spread of oil slicks on a calm sea

,” in

Oil on the Sea: Proceedings of a Symposium on the Scientific and Engineering Aspects of Oil Pollution of the Sea, Sponsored by Massachusetts Institute of Technology and Woods Hole Oceanographic Institution and Held at Cambridge, Massachusetts, May 16, 1969

, edited by

D. P.

Hoult

(

Springer US

,

Boston, MA

,

1969

), pp.

53

–

63

.

Google Scholar

Crossref

21.

T. K.

Fannelop

and

G. D.

Waldman

, “

Dynamics of oil slicks

,”

AIAA J.

10

,

506

–

510

(

1972

).

https://doi.org/10.2514/3.50127

Google Scholar

Crossref

22.

L.

Ovsyannikov

, “

Two-layer “shallow water” model

,”

J. Appl. Mech. Tech. Phys.

20

,

127

–

134

(

1979

).

https://doi.org/10.1007/bf00910010

Google Scholar

Crossref

23.

C.

Vreugdenhil

, “

Two-layer shallow-water flow in two dimensions, a numerical study

,”

J. Comput. Phys.

33

,

169

–

184

(

1979

).

https://doi.org/10.1016/0021-9991(79)90014-7

Google Scholar

Crossref

24.

E.

Audusse

,

M.-O.

Bristeau

,

B.

Perthame

, and

J.

Sainte-Marie

, “

A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation

,”

ESAIM: Math. Modell. Numer. Anal.

45

,

169

–

200

(

2011

).

https://doi.org/10.1051/m2an/2010036

Google Scholar

Crossref

25.

P.

Milewski

,

E.

Tabak

,

C.

Turner

,

R.

Rosales

, and

F.

Menzaque

, “

Nonlinear stability of two-layer flows

,”

Commun. Math. Sci.

2

,

427

–

442

(

2004

).

https://doi.org/10.4310/cms.2004.v2.n3.a5

Google Scholar

Crossref

26.

D.

Lannes

and

M.

Ming

, “

The Kelvin-Helmholtz instabilities in two-fluids shallow water models

,” in

Hamiltonian Partial Differential Equations and Applications

, edited by

P.

Guyenne

,

D.

Nicholls

, and

C.

Sulem

(

Springer

,

New York

,

2015

), pp.

185

–

234

.

Google Scholar

Crossref

27.

A. L.

Stewart

and

P. J.

Dellar

, “

Multilayer shallow water equations with complete Coriolis force. Part 3. Hyperbolicity and stability under shear

,”

J. Fluid Mech.

723

,

289

–

317

(

2013

).

https://doi.org/10.1017/jfm.2013.121

Google Scholar

Crossref

28.

M. Y.

Lam

,

M. S.

Ghidaoui

, and

A. A.

Kolyshkin

, “

The roll-up and merging of coherent structures in shallow mixing layers

,”

Phys. Fluids

28

,

094103

(

2016

).

https://doi.org/10.1063/1.4960391

Google Scholar

Crossref

29.

F.

Bouchut

and

T.

Morales de Luna

, “

An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment

,”

ESAIM: Math. Modell. Numer. Anal.

42

,

683

–

698

(

2008

).

https://doi.org/10.1051/m2an:2008019

Google Scholar

Crossref

30.

M. J.

Castro-Díaz

,

E. D.

Fernández-Nieto

,

J. M.

González-Vida

, and

C.

Parés-Madroñal

, “

Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system

,”

J. Sci. Comput.

48

,

16

–

40

(

2011

).

https://doi.org/10.1007/s10915-010-9427-5

Google Scholar

Crossref

31.

F.

Bouchut

and

V.

Zeitlin

, “

A robust well-balanced scheme for multi-layer shallow water equations

,”

Discrete Contin. Dyn. Syst.—Ser. B

13

,

739

(

2010

).

https://doi.org/10.3934/dcdsb.2010.13.739

Google Scholar

Crossref

32.

A.

Chiapolino

and

R.

Saurel

, “

Models and methods for two-layer shallow water flows

,”

J. Comput. Phys.

371

,

1043

–

1066

(

2018

).

https://doi.org/10.1016/j.jcp.2018.05.034

Google Scholar

Crossref

33.

T. M.

Hatcher

and

J. G.

Vasconcelos

, “

Alternatives for flow solution at the leading edge of gravity currents using the shallow water equations

,”

J. Hydraul. Res.

52

,

228

–

240

(

2014

).

https://doi.org/10.1080/00221686.2013.874376

Google Scholar

Crossref

34.

J. W.

Rottman

and

J. E.

Simpson

, “

Gravity currents produced by instantaneous releases of a heavy fluid in a rectangular channel

,”

J. Fluid Mech.

135

,

95

–

110

(

1983

).

https://doi.org/10.1017/s0022112083002979

Google Scholar

Crossref

35.

M.

Ungarish

, “

Two-layer shallow-water dam-break solutions for gravity currents in non-rectangular cross-area channels

,”

J. Fluid Mech.

732

,

537

–

570

(

2013

).

https://doi.org/10.1017/jfm.2013.417

Google Scholar

Crossref

36.

V.

Lombardi

,

C.

Adduce

,

G.

Sciortino

, and

M.

La Rocca

, “

Gravity currents flowing upslope: Laboratory experiments and shallow-water simulations

,”

Phys. Fluids

27

,

016602

(

2015

).

https://doi.org/10.1063/1.4905305

Google Scholar

Crossref

37.

M.

Bjørnestad

and

H.

Kalisch

, “

Shallow water dynamics on linear shear flows and plane beaches

,”

Phys. Fluids

29

,

073602

(

2017

).

https://doi.org/10.1063/1.4994593

Google Scholar

Crossref

38.

T.

Zemach

,

M.

Ungarish

,

A.

Martin

, and

M. E.

Negretti

, “

On gravity currents of fixed volume that encounter a down-slope or up-slope bottom

,”

Phys. Fluids

31

,

096604

(

2019

).

https://doi.org/10.1063/1.5121290

Google Scholar

Crossref

39.

L.

Fraccarollo

and

E. F.

Toro

, “

Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems

,”

J. Hydraul. Res.

33

,

843

–

864

(

1995

).

https://doi.org/10.1080/00221689509498555

Google Scholar

Crossref

40.

J.

Zhou

,

D.

Causon

,

C.

Mingham

, and

D.

Ingram

, “

The surface gradient method for the treatment of source terms in the shallow-water equations

,”

J. Comput. Phys.

168

,

1

–

25

(

2001

).

https://doi.org/10.1006/jcph.2000.6670

Google Scholar

Crossref

41.

P. G.

LeFloch

and

M. D.

Thanh

, “

The Riemann problem for the shallow water equations with discontinuous topography

,”

Commun. Math. Sci.

5

,

865

–

885

(

2007

).

https://doi.org/10.4310/cms.2007.v5.n4.a7

Google Scholar

Crossref

42.

Q.

Liang

and

F.

Marche

, “

Numerical resolution of well-balanced shallow water equations with complex source terms

,”

Adv. Water Resour.

32

,

873

–

884

(

2009

).

https://doi.org/10.1016/j.advwatres.2009.02.010

Google Scholar

Crossref

43.

P. G.

LeFloch

and

M. D.

Thanh

, “

A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime

,”

J. Comput. Phys.

230

,

7631

–

7660

(

2011

).

https://doi.org/10.1016/j.jcp.2011.06.017

Google Scholar

Crossref

44.

J.

Murillo

and

A.

Navas-Montilla

, “

A comprehensive explanation and exercise of the source terms in hyperbolic systems using Roe type solutions. Application to the 1D-2D shallow water equations

,”

Adv. Water Resour.

98

,

70

–

96

(

2016

).

https://doi.org/10.1016/j.advwatres.2016.10.019

Google Scholar

Crossref

45.

R. J.

LeVeque

,

Finite-Volume Methods for Hyperbolic Problems

(

Cambridge University Press

,

2002

).

Google Scholar

Crossref

46.

R.

Monjarret

, “

Local well-posedness of the two-layer shallow water model with free surface

,”

SIAM J. Appl. Math.

75

,

2311

–

2332

(

2015

).

https://doi.org/10.1137/140957020

Google Scholar

Crossref

47.

V. V.

Ostapenko

, “

Complete systems of conservation laws for two-layer shallow water models

,”

J. Appl. Mech. Tech. Phys.

40

,

796

–

804

(

1999

).

https://doi.org/10.1007/bf02468461

Google Scholar

Crossref

48.

V.

Ostapenko

, “

Stable shock waves in two-layer shallow water

,”

J. Appl. Math. Mech.

65

,

89

–

108

(

2001

).

https://doi.org/10.1016/s0021-8928(01)00010-7

Google Scholar

Crossref

49.

G.

Whitham

,

Linear and Nonlinear Waves

(

Wiley

,

1974

).

Google Scholar

50.

H.

Holden

and

N. H.

Risebro

,

Front Tracking for Hyperbolic Conservation Laws

(

Springer

,

2015

), Vol. 152.

Google Scholar

Crossref

51.

D.

Borthwick

, “

Weak solutions

,” in

Introduction to Partial Differential Equations

(

Springer International Publishing

,

Cham

,

2016

), Chap. 10, pp.

177

–

204

.

Google Scholar

Crossref

52.

W. J. M.

Rankine

, “

XV. On the thermodynamic theory of waves of finite longitudinal disturbance

,”

Philos. Trans. R. Soc. London

160

,

277

–

288

(

1870

).

https://doi.org/10.1098/rstl.1870.0015

Google Scholar

Crossref

53.

H.

Hugoniot

, “

Mémoire sur la propagation du mouvement dans un fluide indéfini (Première partie)

,”

J. Math. Pures Appl.

4

,

477

–

492

(

1887

); available at http://sites.mathdoc.fr/JMPA/PDF/JMPA_1887_4_3_A11_0.pdf.

Google Scholar

54.

H.

Hugoniot

, “

Mémoire sur la propagation des mouvements dans les corps et spécialement dans les gas parfaits. Deuxième partie

,”

J. Ec. Polytech.

57

,

1

–

126

(

1889

); available at https://gallica.bnf.fr/ark:/12148/bpt6k4337130/f5.image.

Google Scholar

55.

H. E.

Huppert

and

J. E.

Simpson

, “

The slumping of gravity currents

,”

J. Fluid Mech.

99

,

785

–

799

(

1980

).

https://doi.org/10.1017/s0022112080000894

Google Scholar

Crossref

56.

J.

Priede

, “

Self-contained two-layer shallow water theory of strong internal bores

,” e-print arXiv:1806.06041 [physics.flu-dyn] (

2019

).

Google Scholar

57.

Z.

Borden

and

E.

Meiburg

, “

Circulation based models for Boussinesq gravity currents

,”

Phys. Fluids

25

,

101301

(

2013

).

https://doi.org/10.1063/1.4825035

Google Scholar

Crossref

58.

J.

Boussinesq

,

Théorie Analytique de la Chaleur: Mise en Harmonie avec la Thermodynamique et avec la Théorie Mécanique de la Lumière

(

Gauthier-Villars

,

1903

), Vol. 2.

Google Scholar

59.

M.

Ungarish

, “

Two-layer shallow-water dam-break solutions for non-Boussinesq gravity currents in a wide range of fractional depth

,”

J. Fluid Mech.

675

,

27

–

59

(

2011

).

https://doi.org/10.1017/s0022112010006397

Google Scholar

Crossref

60.

R. J.

Lowe

,

J. W.

Rottman

, and

P. F.

Linden

, “

The non-Boussinesq lock-exchange problem. Part 1. Theory and experiments

,”

J. Fluid Mech.

537

,

101

–

124

(

2005

).

https://doi.org/10.1017/s0022112005005069

Google Scholar

Crossref

61.

V. K.

Birman

,

J. E.

Martin

, and

E.

Meiburg

, “

The non-Boussinesq lock-exchange problem. Part 2. High-resolution simulations

,”

J. Fluid Mech.

537

,

125

–

144

(

2005

).

https://doi.org/10.1017/s0022112005005033

Google Scholar

Crossref

62.

D.

Sher

and

A. W.

Woods

, “

Gravity currents: Entrainment, stratification and self-similarity

,”

J. Fluid Mech.

784

,

130

–

162

(

2015

).

https://doi.org/10.1017/jfm.2015.576

Google Scholar

Crossref

63.

N.

Dunford

and

J. T.

Schwartz

,

Linear Operators, Part 1: General Theory

(

New York Interscience

,

1957

).

Google Scholar

64.

W.

Rudin

,

Principles of Mathematical Analysis

(

McGraw-Hill Professional

,

1976

).

Google Scholar

65.

C.

Doran

and

A.

Lasenby

,

Geometric Algebra for Physicists

(

Cambridge University Press

,

2003

).

Google Scholar

Crossref

66.

T. M.

Hatcher

and

J. G.

Vasconcelos

, “

Finite-volume and shock-capturing shallow water equation model to simulate Boussinesq-type lock-exchange flows

,”

J. Hydraul. Eng.

139

,

1223

–

1233

(

2013

).

https://doi.org/10.1061/(asce)hy.1943-7900.0000775

Google Scholar

Crossref

67.

F.

Marche

, “

Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects

,”

Eur. J. Mech.: B/Fluids

26

,

49

–

63

(

2007

).

https://doi.org/10.1016/j.euromechflu.2006.04.007

Google Scholar

Crossref

68.

V.

Joshi

and

R. K.

Jaiman

, “

An adaptive variational procedure for the conservative and positivity preserving Allen-Cahn phase-field model

,”

J. Comput. Phys.

366

,

478

–

504

(

2018

).

https://doi.org/10.1016/j.jcp.2018.04.022

Google Scholar

Crossref

69.

S.

Soares-Frazão

,

R.

Canelas

,

Z.

Cao

,

L.

Cea

,

H. M.

Chaudhry

,

A. D.

Moran

,

K. E.

Kadi

,

R.

Ferreira

,

I. F.

Cadórniga

,

N.

Gonzalez-Ramirez

,

M.

Greco

,

W.

Huang

,

J.

Imran

,

J. L.

Coz

,

R.

Marsooli

,

A.

Paquier

,

G.

Pender

,

M.

Pontillo

,

J.

Puertas

,

B.

Spinewine

,

C.

Swartenbroekx

,

R.

Tsubaki

,

C.

Villaret

,

W.

Wu

,

Z.

Yue

, and

Y.

Zech

, “

Dam-break flows over mobile beds: Experiments and benchmark tests for numerical models

,”

J. Hydraul. Res.

50

,

364

–

375

(

2012

).

https://doi.org/10.1080/00221686.2012.689682

Google Scholar

Crossref

70.

J. G.

Zhou

,

D. M.

Causon

,

C. G.

Mingham

, and

D. M.

Ingram

, “

Numerical prediction of dam-break flows in general geometries with complex bed topography

,”

J. Hydraul. Eng.

130

,

332

–

340

(

2004

).

https://doi.org/10.1061/(asce)0733-9429(2004)130:4(332)

Google Scholar

Crossref

71.

A.

Ritter

, “

Die fortpflanzung der wasserwellen

,”

Z. Ver. Dtsch. Ing.

36

,

947

–

954

(

1892

).

Google Scholar

72.

W.

Suchon

, “

An experimental investigation of oil spreading over water

,” M.Sc. thesis,

Massachusetts Institute of Technology

,

1970

.

Google Scholar

73.

H.-R.

Chang

,

R.

Reid

, and

J.

Fay

, “

Boiling and spreading of liquid nitrogen and liquid methane on water

,”

Int. Commun. Heat Mass Transfer

10

,

253

–

263

(

1983

).

https://doi.org/10.1016/0735-1933(83)90010-6

Google Scholar

Crossref

74.

H.-R.

Chang

and

R.

Reid

, “

Spreading-boiling model for instantaneous spills of liquefied petroleum gas (LPG) on water

,”

J. Hazard. Mater.

7

,

19

–

35

(

1982

).

https://doi.org/10.1016/0304-3894(82)87002-7

Google Scholar

Crossref

75.

D. S.

Burgess

,

J.

Murphy

, and

M.

Zabetakis

, “

Hazards of LNG spillage in marine transportation

,” Technical Report AD0705078,

Safety Research Center, Bureau of Mines

,

Pittsburgh, PA

,

1970

.

Google Scholar

76.

E. F.

Toro

and

S. J.

Billett

, “

Centred TVD schemes for hyperbolic conservation laws

,”

IMA J. Numer. Anal.

20

,

47

–

79

(

2000

).

https://doi.org/10.1093/imanum/20.1.47

Google Scholar

Crossref

77.

D. I.

Ketcheson

and

A. C.

Robinson

, “

On the practical importance of the SSP property for Runge-Kutta time integrators for some common Godunov-type schemes

,”

Int. J. Numer. Methods Fluids

48

,

271

–

303

(

2005

).

https://doi.org/10.1002/fld.837

Google Scholar

Crossref

78.

J.

Smoller

,

Shock Waves and Reaction-Diffusion Equations

(

Springer-Verlag

,

New York

,

1983

).

Google Scholar

Crossref

2019

Author(s)

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

I. INTRODUCTION

II. THEORY OF THE SHALLOW-WATER EQUATIONS

A. The one-layer shallow-water equations

B. The two-layer shallow-water equations

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

D. The Rankine-Hugoniot condition

E. Physical solutions

III. DERIVATION OF FROUDE NUMBERS FROM THE 2LSWE

IV. REDUCING THE TWO-LAYER SYSTEMS TO EFFECTIVE ONE-LAYER SYSTEMS

A. The constant-height lemma

B. The one-layer approximation theorem

C. Discussion of the theorem

V. CASES FOR THE ONE-DIMENSIONAL DAM-BREAK PROBLEM

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

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

C. Case Rb: Effect of nonzero depth on both sides of dam

D. Case Rc: Comparison to dam-break experiments

VI. NUMERICAL RESULTS

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

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

C. Case Rb: Effect of nonzero depth at both sides of dam

D. Case Rc: Comparison to spreading experiments

VII. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: DERIVING RANKINE-HUGONIOT CONDITIONS

APPENDIX B: RANKINE-HUGONIOT CONDITIONS FOR 2LSWE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

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

I. INTRODUCTION

II. THEORY OF THE SHALLOW-WATER EQUATIONS

A. The one-layer shallow-water equations

B. The two-layer shallow-water equations

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

D. The Rankine-Hugoniot condition

E. Physical solutions

III. DERIVATION OF FROUDE NUMBERS FROM THE 2LSWE

IV. REDUCING THE TWO-LAYER SYSTEMS TO EFFECTIVE ONE-LAYER SYSTEMS

A. The constant-height lemma

B. The one-layer approximation theorem

C. Discussion of the theorem

V. CASES FOR THE ONE-DIMENSIONAL DAM-BREAK PROBLEM

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

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

C. Case Rb: Effect of nonzero depth on both sides of dam

D. Case Rc: Comparison to dam-break experiments

VI. NUMERICAL RESULTS

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

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

C. Case Rb: Effect of nonzero depth at both sides of dam

D. Case Rc: Comparison to spreading experiments

VII. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: DERIVING RANKINE-HUGONIOT CONDITIONS

APPENDIX B: RANKINE-HUGONIOT CONDITIONS FOR 2LSWE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only