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^{1} 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,^{8} gravity currents particularly in geophysical systems,^{9–12} monomolecular layers for evaporation control,^{13} and coalescence in three-phase fluid systems.^{14} 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}

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 − *ρ*_{1}/*ρ*_{2})*g*, where *ρ*_{1} and *ρ*_{2} are the two fluid densities and *ρ*_{1} < *ρ*_{2}.

An early result for the Froude number of gravity currents was presented by von Kármán.^{17} They found that for the edge of a spreading of gravity current at semi-infinite depth, $FrLE=2$, where the subscript is short for the “leading edge.” Benjamin^{18} later developed a model for Fr_{LE} for spreading of gravity currents with constant height,

where *α* = *h*_{2}/(*h*_{1} + *h*_{2}). Here, *h*_{1} and *h*_{2} are the heights of the top and bottom layers, respectively. This model approaches the result by von Kármán when the bottom layer becomes thin, *h*_{2} ≪ *h*_{1}. More recently, Ungarish^{19} extended the result of Benjamin to the spreading of gravity currents into a lighter fluid with an open surface. This result also gives $FrLE=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,^{20} 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

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.^{24} 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.^{25} 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.^{29} Attempts to amend the nonhyperbolicity of the systems include adding numerical (nonphysical) friction forces,^{30} operator-splitting approaches,^{31} and introduction of an artificial compressibility.^{32}

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^{34} 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.

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.

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,^{45}

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,

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 $FrLE=2$. This is exactly the result by von Kármán^{17} 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.^{31}

### B. The two-layer shallow-water equations

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

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.

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

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

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

where

and where we have defined

### 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

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 $\u301af\u301b\u2261fr\u2212fl$, where $fr\u2261lim\epsilon \u21920+f(\xi +\epsilon n^)$ and $fl\u2261lim\epsilon \u21920\u2212f(\xi +\epsilon n^)$. The Rankine-Hugoniot condition then states that the discontinuity at any point $\xi $ ∈ Γ propagates along the outward-pointing normal vector $n^$ with a speed *S*. This speed is called the shock speed and satisfies the relation

Similarly, if * u* satisfies a general vector conservation equation,

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

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* · ∇

*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 ∝*

**u***h*

_{1}∇

*h*

_{2}makes the normal component for the momentum equations ill-defined, which is why we must exclude this formulation of the 2LSWE from the analysis.

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

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

and

Finally, we note that in calculations with the Rankine-Hugoniot condition it is useful to observe that $\u301aab\u301b=\u301aa\u301bb+a\u301ab\u301b$, where $a=(al+ar)/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.^{50} However, the word *entropy* is typically restricted to convex functions of the solution variables. As has been showed by Ostapenko,^{47} 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

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

where **q**_{i} = *ρ*_{i}*h*_{i}**u**_{i} for short.

## III. DERIVATION OF FROUDE NUMBERS FROM THE 2LSWE

In the following, we briefly illustrate the surprising effectiveness of the 2LSWE to predict shock speeds despite its underlying assumption of negligible vertical acceleration. To do this, we apply the Rankine-Hugoniot conditions and the 2LSWE to derive expressions for the leading edge Froude number (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.

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.^{58} The most commonly used front condition applied to such flows is the equation for the Froude number given by Benjamin^{18} [Eq. (2)]. For instance, Ungarish^{59} applied the Froude number by Benjamin as a boundary condition when solving the 2LSWE for rectangular geometries. They also generalized this to arbitrary geometries.^{35}

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

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

These are the same equations that were used by Rottman and Simpson^{34} 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*_{1} + *h*_{2} constant, the sum of Eqs. (19a) and (19b) implies that *h*_{1}*u*_{1} + *h*_{2}*u*_{2} is constant in *x*. If we assume that the total momentum is 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*_{1}*u*_{1} + *h*_{2}*u*_{2} = 0.

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

After some algebraic manipulation, we find that

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

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

Figure 4 compares our results from the 2LSWE, Eq. (23), to the model by Benjamin,^{18} 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}

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^{62} 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^{56} 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*_{1} − *h*_{2} and *ϑ* ≔ *u*_{1} − *u*_{2}, 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 $f\xaf$ 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*

*i*∈ {1, 2} in a two-layer shallow-water system is source-bounded if there exists $K\u2208R$ such that the source terms satisfy ∀

*h*

_{i}>

*K*,

*j*= 1 and

*j*= 2.

*Lemma 1.*

*Let*(

*s*,

*d*) = (1, 2)

*or*(2, 1), ${Dk}k\u2208N$

*be a sequence of increasing real numbers, h*

_{0}

*and f be scalar functions, and*

*q*_{1,0}

*and*

*q*_{2,0}

*be vector functions. Furthermore, consider a 2LSWE system with initial conditions,*

*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*${(h1k,h2k,q1k,q2k)}k\u2208N$

*be physical solutions to the 2LSWE system. If*${(hsk,hdk\u2212Dk,qsk,qdk/Dk)}k\u2208N$

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

*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.^{63} This implies that we can interchange the order of limits and differentiation.^{64} 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 $hi\xaf=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 $Ek\u0303=2Ek/Dk2$ be a scaled energy, and it follows by insertion that $Ek\u0303(0,x)\u2192g\rho d$ as *k* → ∞. That is, the scaled energy $Ek\u0303(0,x)$ approaches the minimal for a system with $hdk(t,x)\u2215Dk\xaf=1$ in the limit when *k* → ∞. Source-boundedness of the mass source terms implies that $hdk\u2215Dk\xaf\u21921$ for all $t\u2208R$ as *k* → ∞.

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

*ε*> 0 and $\u2200N\u2208N$, ∃

*k*>

*N*, such that

*h*

_{dk}(

*t*,

*x*)/

*D*

_{k}deviates from 1 by a term which does not vanish in the limit

*k*→ ∞. This contradicts Eq. (24), as discussed above, so

*h*

_{dk}(

*t*,

*)/*

**x***D*

_{k}→ 1 for all $t\u2208R$ and $x\u2208R2$ as

*k*→ ∞.$\u25a1$

### 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

*Let*(

*s*,

*d*) = (1, 2)

*or*(2, 1)

*,*${Dk}k\u2208N$

*be a sequence of increasing real numbers, h*

_{0}

*and f be scalar functions, and*

*q*_{1,0}

*and*

*q*_{2,0}

*be vector functions, all defined on*$\Omega \u2286Rn$

*. Furthermore, consider 2LSWE in the form of*Eq. (7)

*or*Eq. (8)

*with initial conditions*

*in which layer d is source-bounded and the density of layer d is constant. Now, let*${(h1k,h2k,q1k,q2k)}k\u2208N$

*be physical solutions to the 2LSWE such that q*

_{dk}

*satisfies the boundary condition,*

*with*$K\u2208R$

*independent of k and where ∂*Ω

*may be at infinity*.

*If* ${(hsk,hdk\u2212Dk,qsk,qdk/Dk)}k\u2208N$ *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*$hdk\u2212Dk\u2192C\u2212\rho sd\u22121hs\u2212\rho 2d\u22121b/\rho dd\u22121$

*, 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=\rho sd\u22121hs+\rho dd\u22121f+\rho 2d\u22121b\xaft=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.^{50}

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^{65} 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.

**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. (7) yields

_{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*

^{−1}=

**e**_{m}∧⋯∧

**e**_{1}, 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*→ ∞.

**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

^{65}

*L*is any linear function,

*I*

_{m}is the pseudoscalar of the tangent space to

*V*, and

*I*

_{m−1}(

*) is the pseudoscalar of the tangent space to*

**x***∂V*at

*. The vector derivative is ∇ ⊙*

**x***A*= ∇ ·

*A*+ ∇ ∧

*A*, and the overdot indicates where it acts. That is, the integrand on the right-hand side of Eq. (29) is ∑

_{i}

*∂*

_{i}

*L*(

*e*

_{i}⊙

*I*

_{m}). See, for instance, the textbook by Doran and Lasenby.

^{65}

*V*be some region where

*u*

_{dk}is differentiable and let

*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 (29) then states that

The surface integral can be decomposed into one vector component whose integrand is proportional to $udk\u22c5n^$ and one triplet-vector component whose integrand is proportional to $u\u2227n^$. 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 $limk\u2192\u221en^\u22c5udk$ 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.

*k*→ ∞, $\rho sd\u22121hs+\rho dd\u22121(hdk\u2212Dk)+\rho 2d\u22121b$ is constant in space. Finally, plugging this into the equation for

**u**_{s}in Eq. (7) or Eq. (8), we get Eq. (25). In the regions where the solution is differentiable, the various formulations of the 2LSWE, Eqs. (6)–(8), are equivalent. This completes the proof of (i).

*k*→ ∞ to the Hugoniot locus of Eq. (25). Let $\gamma :=1/hd$. In Appendix B, we show that the full set of Rankine-Hugoniot conditions for the 2LSWE [Eqs. (7) and (8)] may be written as

*d*= 1,

*d*= 2,

*g*

_{1},

*g*

_{2}, and

*g*

_{3}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 $\rho 1d\u22121h1+\rho 2d\u22121h2$ 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.$\u25a1$

### 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 − *ρ*_{1}/*ρ*_{2})*g* according to Eq. (25). The approximation becomes more accurate when the depth of the deepest layer is increased without increasing momentum or other key properties. Figure 5 shows a sketch of how the two-layer cases converge to one-layer cases when we increase the “depth,” *D*_{k}.

The interesting part about Theorem 1 is not that smooth solutions of the 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,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 $FrLE=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 $FrLE=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*_{1}∇*h*_{2} well-defined. The equations can then be investigated numerically by studying how shocks emerge when the viscosity coefficient is reduced. They can also be investigated analytically by looking at 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*.^{45} In the *two-layer* dam-break problem, a lighter fluid of height *h*_{1} spreads on top of a heavier fluid of height *h*_{2} as shown in Fig. 6. The problem has been frequently used in the literature as a benchmark case for spreading models.^{68–70}

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.

### 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

where *h*_{0} is constant.

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

where $c0=\delta gh0$. 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

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 $FrLE=2$ at the leading edge.

### 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 (*∂*_{x}*h*)(*x* = −*L*, *t*) = 0 and *u*(*x* = −*L*, *t*) = 0. The reflective wall removes the self-similarity of the solution, which enables a study of how the accuracy of the one-layer approximation from Theorem 1 evolves in time.

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

where again *h*_{0} 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,

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,^{72} 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.

Authors . | Experiment . | Spill volume (L) . | Height h_{0} (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_{0} (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.^{73} 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.^{74} 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^{73} 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.^{74} 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^{73} and include constant evaporation rates of 0.16 kg m^{−2} s^{−1} for methane and 0.201 kg m^{−2} s^{−1} for nitrogen in the mass balances. These values are the same as those used by Chang, Reid, and Fay, which are based on experimental studies of the relevant substances.^{75} 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^{−3} for nitrogen and 254 kg m^{−3} for methane. Similar to Chang, Reid, and Fay, we will use reduced densities in our simulations as well. Finally, Chang, Reid, and Fay 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^{76} and the second-order MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) reconstruction with a minmod limiter^{45} in each finite volume. The solutions are advanced in time with a standard third-order three-stage strong stability-preserving Runge-Kutte method.^{77} Although some of the equations have terms that are in general not conservative, e.g., the convective term in the locally conservative 2LSWE (7), it should be noted that these become conservative when the equations are restricted to a single spatial dimension. The Courant-Friedrichs-Lewy (CFL)-number is 0.9 for all cases.

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

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\u2192\lambda \u22121t$ and $us\u2192\lambda us$ 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*_{0} = 1 m and *δ* = 0.7 with 400 grid cells on the domain *x*/*c*_{0}*t* ∈ (−2, 2). The initial depth, *H*, was increased stepwise from the initial value, *H* = *h*_{0}, to obtain a larger height difference between the layers. The results presented in 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*.

Figure 10 shows the normalized difference in *L*^{1} for the top-layer height between the analytic solution [Eq. (43)] and the solutions obtained with the 2LSWE, $h\u03031$ and *h*_{1}, respectively ($\u2225h1\u2212h\u03031\u22251/h01$). The differences are shown as a function of the initial depth, *H*. The circles correspond to the locally conservative 2LSWE (7), the triangles correspond to the globally conservative 2LSWE (8), and the solid line indicates a slope of −1. The plot shows that the difference between the globally and the locally conservative 2LSWE is small as expected. Other relevant variables, such as *h*_{2}, *u*_{1}, and *u*_{2}, were found to exhibit a similar behavior.

### 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*_{0} = 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*_{2} is constant. This is done to show that a perturbation of the initial solution of the relatively deep layer does not prevent the 2LSWE to converge to the one-layer approximation when the depth increases. As the depths increase, we see that the solutions of the 2LSWE converge toward the locally conservative 1LSWE as predicted by the theorem.

To quantify how the solutions of the 2LSWE converge to those of the locally conservative 1LSWE, we will compare the solutions at various times, initial depths, *H*, and initial widths, *L*. We consider solutions of the globally conservative 2LSWE and the locally conservative 1LSWE and evaluate two quantifiable differences. In Fig. 12(a), the *L*^{1} difference of the top-layer height, $\u2225h12LSWE\u2212h11LSWE\u22251/h0L$, is plotted for varying times and depths, *H*/*h*_{0}. Figure 12(b) shows the difference of the *leading edge* position at *t* = 5 s, |*r*_{2LSWE} − *r*_{1LSWE}|/*h*_{0}, for varying initial depths, *H*/*h*_{0}, and widths, *L*/*h*_{0}. The position of the leading edge is here defined as the smallest *x*-value where the top layer is thinner than 10^{−4} m. Figure 12 shows that the differences in *h*_{1} decrease with time, which is reasonable since the spreading fluid becomes gradually thinner. As expected, the differences decrease with increasing value of *H*/*h*_{0}. Similar to case 0, the errors in the variables *h*_{2}, *u*_{1}, and *u*_{2} as quantified by the *L*^{1} norm exhibit the same trends as the top layer height (not shown). 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.

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.

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

### 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).

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

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*_{0}.

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*^{1} 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.

Finally, we note that Chang, Reid, and Fay^{73} 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 $2$ from below. Moreover, our analysis shows that $FrLE=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 $FrLE\u2208[1,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 − *ρ*_{1}/*ρ*_{2})*g*, the one-layer model correctly captures the behavior of shocks and contact discontinuities. In particular, the one-layer model results in $FrLE=2$, which is exactly the same as the theoretical predictions by von Kármán,^{17} Benjamin,^{18} and Ungarish^{19} 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.^{78} Consider a shock along Γ which is normal to $n^$ and let *ϕ* be a test function with compact support *D*, which lies in the *x*_{n}*t*-plane.

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

where the subscript on *ϕ* denotes partial differentiations and

where $uiT$ 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*_{1} and *D*_{2}. Because the solution is differentiable inside these regions, we may use Green’s theorem and obtain

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,

For the globally conservative system, a similar treatment yields

and

The treatment presented above is not applicable for the tangential component of the velocity equations because they involve a term on the form $n^\u22c5ui\u2202nuiT$. 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 $\u301auiT\u301b=0$ has the additional advantage of making $n^\u22c5ui\u2202nuiT$ well defined as the product of $n^\u22c5ui$ and $\u2202nuiT$. Otherwise the distribution $n^\u22c5ui\u2202nuiT$ cannot be decomposed without relying on some mollification scheme.

Ostapenko^{48} 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 $\u301auiT\u301b=0$. A conservation law for the quantity *∂*_{1}*u*_{i,2} − *∂*_{2}*u*_{i,1} ≡ $w$_{i} can be obtained by taking distributional derivatives of the different components of the velocity equation, yielding

where

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

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,

Here, *g*_{1} and *g*_{2} differ for Eqs. (7) and (8), while *g*_{3} will be the same. Furthermore, we will show that all of *g*_{1}, *g*_{2}, and *g*_{3} vanish when *γ* = 0. Note that Eq. (B1a) follows directly from Eq. (11) applied to the shallowest layer.

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

where we used that $\u301aab\u301b=\u301aa\u301bb+a\u301ab\u301b$. This gives

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

that is,

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

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

with

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 $n^$ and use Eq. (B2), we obtain

that is,

and

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

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

which gives