Hydrodynamic models and confinement effects by horizontal boundaries

TL;DR

This paper investigates how rigid horizontal boundaries influence ideal fluid dynamics, revealing that contact with boundaries can cause finite-time singularities in the models, with analysis supported by numerical simulations.

Contribution

It analyzes the effects of boundary contact on long-wave models and Euler systems, highlighting conditions leading to singularities and providing numerical validation.

Findings

Contact with boundaries can cause finite-time singularities.

Models predict qualitative and quantitative effects of confinement.

Numerical simulations support theoretical predictions.

Abstract

Confinement effects by rigid boundaries in the dynamics of ideal fluids are considered from the perspective of long-wave models and their parent Euler systems, with the focus on the consequences of establishing contacts of material surfaces with the confining boundaries. When contact happens, we show that the model evolution can lead to the dependent variables developing singularities in finite time. The conditions and the nature of these singularities are illustrated in several cases, progressing from a single layer homogeneous fluid with a constant pressure free surface and flat bottom, to the case of a two-fluid system contained between two horizontal rigid plates, and finally, through numerical simulations, to the full Euler stratified system. These demonstrate the qualitative and quantitative predictions of the models within a set of examples chosen to illustrate the theoretical results.

Figures

Equations

$${\mathbf{v}}_t+\mathbf{v}\cdot \nabla \mathbf{v}=-\frac{\nabla p}{\rho }-g\mathbf{k},\nabla \cdot \mathbf{v}=0,{\rho }_t+\mathbf{v}\cdot \nabla \rho =0.$$

$$\frac{\partial p}{\partial z}\to -g\rho ,\text{and}\rho \to {\rho }_b(z),\text{ for }|x|\to \infty ,w(x,0)=w(x,h)=0\text{for all x}.$$

$$\begin{array}{rl}& {{\eta }_i}_t+({\overline{u}}_i{\eta }_i{)}_x=0,i=1,2\\ & {{\overline{u}}_1}_t+{\overline{u}}_1{{\overline{u}}_1}_x-g{{\eta }_1}_x+\frac{P_x}{{\rho }_1}+D_1=0,\\ & {{\overline{u}}_2}_t+{\overline{u}}_2{{\overline{u}}_2}_x+g{{\eta }_2}_x+\frac{P_x}{{\rho }_2}+D_2=0,\\ & {\eta }_1+{\eta }_2=h,({\eta }_1{\overline{u}}_1+{\eta }_2{\overline{u}}_2{)}_x=0.\end{array}$$

$${\overline{u}}_1(x,t)=\frac{1}{{\eta }_1(x,t)}{\int }_{h-{\eta }_1(x,t)}^h{u}_1(x,z,t)\mathrm{d}z,{\overline{u}}_2(x,t)=\frac{1}{{\eta }_2(x,t)}{\int }_0^{{\eta }_2(x,t)}{u}_2(x,z,t)\mathrm{d}z,$$

$$D_i=\frac{1}{3{\eta }_i}{\partial }_x[{\eta }_i^3({{\overline{u}}_i}_{xt}+{\overline{u}}_i{{\overline{u}}_i}_{xx}-({{\overline{u}}_i}_x{)}^2)]+\dots ,i=1,2,$$

$$\epsilon \equiv \underset{(x,t)}{\max }[{\eta }_i/L]\ll 1,i=1,2$$

$$\sigma (x,t)\equiv \frac{1}{h}\left({\rho }_2{\overline{u}}_2-{\rho }_1{\overline{u}}_1\right).$$

$$\begin{array}[]{l}{\displaystyle{\eta_{t}=}}\,{\displaystyle{-\partial_{x}\left({{h\eta(\eta-h)\sigma\over(\rho_{2}-\rho_{1})\eta-\rho_{2}h}}\right)}}\,,\vspace*{0.3cm}\\ {\displaystyle{{\sigma}_{t}}}={\displaystyle{-\partial_{x}\left({h\over 2}\,{\rho_{2}(h-\eta)^{2}-\rho_{1}\eta^{2}\over((\rho_{2}-\rho_{1})\eta-\rho_{2}h)^{2}}\,\sigma^{2}+\frac{g(\rho_{2}-\rho_{1})\eta}{h}\right)}}\,.\end{array}$$

$$\begin{array}[]{l}\eta_{t}=-\partial_{x}(\eta\,{{\overline{u}}}_{2})\,,\vspace*{0.3cm}\\ {{{\overline{u}}}_{2}}_{t}=-\partial_{x}\left({1\over 2}\,{{{\overline{u}}}_{2}}^{2}+g\eta\right)\,.\end{array}$$

$$\dot{x}(t)=u_2(x(t),\eta (x(t),t),t),x(0)=x_0,$$

$$\begin{array}{rl}{\dot{Z}}_1(t)& ={\eta }_x(x_{\kappa }(t),t){\dot{x}}_{\kappa }(t)+{\eta }_t(x_{\kappa }(t),t)={\eta }_x(x_{\kappa }(t),t)u_2(x_{\kappa }(t),\eta (x_{\kappa }(t),t),t)+{\eta }_t(x_{\kappa }(t),t)\\ & =w_2(x_{\kappa }(t),\eta (x_{\kappa }(t),t),t)=w_2(x_{\kappa }(t),Z_1(t),t).\end{array}$$

$$\dot{z}(t)=w_2(x_{\kappa }(t),z(t),t),z(0)=0.$$

$$\eta \to \eta /h,{\overline{u}}_2\to \frac{{\overline{u}}_2}{\sqrt{gh}},x\to \frac{x}{h},t\to t\sqrt{\frac{g}{h}},$$

$${\eta }_t+(\eta u{)}_x=0,u_t+u{u}_x+{\eta }_x=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\eta (x(t),t)={\eta }_x\frac{\mathrm{d}x}{\mathrm{d}t}+{\eta }_t={\eta }_{x}u-(\eta u{)}_x=-u_x\eta ,$$

$${\dot{\eta }}_{\kappa }(t)=-{\upsilon }_{\kappa }(t){\eta }_{\kappa }(t),{\eta }_k(0)=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}u(x_{\kappa }(t),t)=-{\eta }_x(x_{\kappa }(t),t)=0,$$

$$\eta_{t}(x_{0},t_{0})=-(u_{x}\eta+u\eta_{x})\Big{|}_{x=x_{0},t=t_{0}}=-u_{x}(x_{0},t_{0})\epsilon=O(\epsilon)$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}={\lambda }_{\pm }=u\pm \sqrt{\eta },x(0)=x_0,$$

$$\left(\begin{array}[]{cc}u&\eta\\ 1&u\end{array}\right)$$

$$R_{\pm }=u\pm 2\sqrt{\eta }$$

$$R_{-}=u-2\sqrt{\eta },$$

$${\eta }_t+(3\sqrt{\eta }+c){\eta }_x=0,\eta (x,0)={\eta }_0(x).$$

$$x-(3\sqrt{\eta }+c)t=F(\eta ),$$

$$x-(3\sqrt{\eta }+c)t=F(\eta ),-\frac{3}{2\sqrt{\eta }}t=F^{\prime }(\eta ),$$

$$\begin{array}[]{l}x(\eta)=F(\eta)-2\,\eta\,F^{\prime}(\eta)-{\displaystyle\frac{2c}{3}}\sqrt{\eta}\,\,F^{\prime}(\eta)\,,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr t(\eta)=-{\displaystyle\frac{2}{3}}\sqrt{\eta}\,\,F^{\prime}(\eta)\,.\end{array}$$

$$F^{\prime }({\eta }_c)+2F^{\prime \prime }({\eta }_c){\eta }_c=0$$

$$x_c=F({\eta }_c)-2{\eta }_c{F}^{\prime }({\eta }_c)-\frac{2c}{3}\sqrt{{\eta }_c}{F}^{\prime }({\eta }_c),t_c=-\frac{2}{3}\sqrt{{\eta }_c}{F}^{\prime }({\eta }_c).$$

$$\eta(x,0)=\left\{\begin{array}[]{cccl}M&&&x<-M^{1/\alpha}\\ (-x)^{\alpha}&&&-M^{1/\alpha}\leq x<0\\ 0&&&x\geq 0\end{array}\right.\,$$

$$x-(3\sqrt{\eta }-2\sqrt{M})t=-{\eta }^{1/\alpha }.$$

Full text

1 Introduction

In investigations of the behavior of stratified fluids under gravity with ideal fluid models, it is well known that setups for which constant density level sets (isopycnals) intersect non-vertical rigid boundaries can cause several difficulties that need to be addressed, either by careful considerations of the mathematical properties of the governing equations or by modifying them to include additional physical effects (see e.g., [6, 9], or Carrier & Greenspan [11] and Phillips [33] for the case of sloped boundaries). With this in mind, under these circumstances, it seems important to assess whether the simplicity afforded by the governing equations for ideal fluids can continue to be physically relevant, or whether the assumptions underlying the foundation of these equations need to be modified, and more complex systems become essential.

To this end, in previous work [8, 7], we have focussed on the consequences that isopycnals intersecting horizontal confining surfaces can have on the structure of the incompressible, variable density Euler system, such as its conservation laws and Hamiltonian formalism. In this work, we want to shift our attention to issues that concern more specifically the time evolution of data in this class. Much of the difficulties here stem from the mathematical assumptions underlying the foundation of the governing equations, specifically the property that “fluid particles on a wall stay on the wall ” ([13], p.22). When this assumption is accepted, the shrinking or opening of “islands” of density variations along the confining boundaries becomes a rather subtle issue, even in limiting cases such as layered fluids separated by sharp interfaces (see figure 1 for a schematic). The latter setup in turn admits the well-known water-wave, or air-water, limit, in which the density of the lighter fluid, air, is assumed to be negligible and enters the dynamics only as means of maintaining constant pressure on the water’s free surface.

There are some similarities between the “splash” and “splat” singularities studied in [12, 21, 22] (or their backward-time versions) and the focus of this work on the interaction of a fluid’s density isoline with rigid boundaries. For instance, the splash example provided in [21], due to the reflectional symmetry across a vertical line, can be viewed as interaction of the free surface with a vertical wall. In this respect, it is worth mentioning the exact free-surface solutions for an Euler fluid evolving under gravity pioneered by John [24], by assuming a simple quadratic dependence on the spatial variables of the velocity potential. This leads to parabolic interface shapes that collapse self-similarly onto their axis for certain time limits. It is remarkable that in our investigation we find a similar idea and behavior playing a prominent role for models of interfaces. In fact, rather than with the full Euler system we will be mostly concerned in this work with models that can extract the essential elements of confining-boundary/interface interactions. While a direct investigation of Euler equations could in principle be extended to include boundary issues for similar setups, much insight can be gained by the study of such reduced models that incorporate boundary conditions and physically relevant properties either exactly or approximately, e.g., through layer averaging of the horizontal velocity field. Such is the case for the Airy (shallow water) system for the water wave problem or its strongly nonlinear counterpart for the layered stratifications [29, 20, 4, 5].

Of course, as already remarked earlier, a description of fluids neglecting phenomena such as viscosity, diffusivity, and surface tension is realistically limited from a physical perspective. However, under suitable assumptions on the magnitude of the corresponding parameters, and on the scales of initial data and evolution times, omitting these effects can be expected to capture correctly the qualitative behavior of the fluid, and in special cases even to reproduce it quantitatively. Similar considerations can be made within classes of simplified models, for instance with regard to whether neglecting dispersion in wave propagation can be justified. For example, in many setups the evolution supported by the Airy (nondispersive) shallow-water system can be viewed as providing a “skeleton” for the relevant physics, shedding light on where (and how) singularity generation in the form of shocks may arise, which in turn would signal the need for additional physical mechanisms to continue the evolution beyond limiting events; such may be the situation encountered in the collapsing of dry spots or of density islands. In our opinion, Stoker’s words ([35], page XIX) about the relevance of approximate models such as the shallow water equations “for clarifying some of the mysteries concerning the dynamical causes” still ring true when applied to the present context.

This paper is essentially organized in two parts. The first concentrates on the dynamics of a single layer fluid (under gravity) and on the interaction of its free surface with a flat bottom boundary. The second part focuses mainly on the case of stratified fluids filling a two-dimensional domain between two flat plates. Both parts present a combination of analytical and numerical approaches to the dynamical questions arising in the study of the time evolution for certain classes of initial conditions.

Specifically, after a brief review in §2 of the mathematical and physical setting of a stratified, incompressible Euler fluid between two horizontal boundaries, we state the fundamental kinematic implications of the appropriate no-flux boundary conditions. The peculiar features of shock formation for the Airy’s model when the free surface connects to a “dry spot” at the bottom of fluid layer is investigated in §3 by means of several families of initial data. We first focus on the simple wave case in §3.1. With the aim of extending this simple wave study to more general setups, we then examine, in §3.2, a class of initial data leading to exact solutions that are the counterpart of the well-known self-similar solution of the Hopf equation. These are also a convenient starting point for a perturbative approach to the study of contacts with the upper lid of the lower layer fluid in the (nondispersive) two-fluid system when the upper fluid density is small.

Next, in §4, we look at how the kinematic properties for the full Euler system transfer to generic quasi-linear systems whose structure is representative of long-wave models, such as the general two-layer systems. The study of these properties is carried out in particular by focussing on their behavior near the hyperbolic-elliptic transition curves in the hodograph plane.

In the two-fluid context of §4, general results about the behavior of contact points with hyperbolic-elliptic transition curves are illustrated by the persistence of an interface contact to a rigid boundary, for as long as singularities do not develop in the system’s evolution.

The general features of our study are then illustrated with a suite of numerical simulations in our last section, §5. The simulations of both the models, when closed form solutions are not available, and the parent, stratified Euler equations focus on a class of initial conditions and time scales that isolate the consequences of contact phenomena for density isolines.

2 Physical setup and equations of motion

We consider the Euler equations for an ideal incompressible and inhomogeneous fluid, subject to gravity, in two dimensions:

[TABLE]

Here ${\boldsymbol{v}}=(u,w)$ is the velocity field, $\rho$ and $p$ are the density and pressure fields, respectively, $g$ is the constant gravity acceleration, $\{{\boldsymbol{i}},{\boldsymbol{k}}\}$ are the unit directional vectors for, respectively, horizontal $x$ and vertical $z$ coordinates, and all physical variables depend on spatial coordinates $(x,z)$ and time $t$. The fluid domains we shall consider here are slabs in the $(x,z)$ plane rigidly confined by horizontal plates of infinite extent located at $z=0$ and $z=h$. Our study will focus on two-dimensional dynamics, though it can be generalized to fully three-dimensional cases. The Euler equations (2.1) are then supplemented by initial data ${\boldsymbol{v}}(x,z,0)={\boldsymbol{v}}_{0}(x,z)$, $\rho(x,z,0)=\rho_{0}(x,z)$, for some functions $({\boldsymbol{v}}_{0},\rho_{0})$, and by boundary conditions

[TABLE]

The first two conditions in (2.2) impose that the fluid at the far ends of the channel be in hydrostatic equilibrium.

A simplification of system (2.1) which retains the essential properties of stratification can be obtained by considering a system of two fluids of homogeneous densities whose interface is described by a smooth function $z=\eta(x,t)=\eta_{2}(x,t)$, see figure 2. As well known, taking layer-averages in this case offers a convenient starting point to derive simplified models. This results in the (non-closed) system

[TABLE]

Here ${{\overline{u}}}_{1,2}$ are the layer-mean velocities, defined as

[TABLE]

and $P(x,t)$ is the interfacial pressure. The terms $D_{1}$, $D_{2}$ at the right hand side of system (2.3) are

[TABLE]

where dots represent terms with nonlocal dependence on the averaged velocities. These terms collect the non-hydrostatic correction to the pressure field, and make the evolution of system (2.3) dispersive, as can be readily seen by linearizing the equations around a flat interface with just the explicit terms in (2.4). When an asymptotic expansion based on the long-wave assumption, e.g.,

[TABLE]

(with $L$ denoting a typical horizontal scale of the initial data), is carried out, expressions (2.4) explicitly define the leading order dispersive terms in the small parameter $\varepsilon$; truncating equations (2.3) at this order results in the strongly nonlinear system studied in [30, 14], which is a two-layer version of the so-called Serre-Green-Naghdi model (see, e.g., [34, 36, 23, 26]) for a single layer with a free surface.

The constraints in the last line of (2.3) can be used to reduce this model to a system of only two unknowns (see, e.g., [7]), which simplifies considerably if the dispersive terms $D_{1}$ and $D_{2}$ are neglected altogether. These non-dispersive equations of motion are most naturally expressed in terms of the thickness $\eta(x,t)\equiv\eta_{2}(x,t)$ of the lower layer and the layer-averaged weighted vorticity

[TABLE]

In these variables the resulting equations read (see e.g., [7]):

[TABLE]

When the limit of vanishing upper fluid density is further considered, together with the nondispersive limit obtained by dropping the dispersive terms $D_{2}$, system (2.3) reduces to the well known shallow water Airy’s model, as $\rho_{1}=0$ implies $P_{x}=0$ (consistently with the condition of constant pressure at the free surface in the classical water wave problem) and the lower layer $i=2$ equations read simply

[TABLE]

In general, the nondispersive approximations above can be expected to lose accuracy as time progresses when steep gradients in the dependent variables become relevant. As we shall see, for the classes of initial data and time scales we study the relative effects may typically be considered negligible when tested against numerical solutions of the parent Euler equations. An exception occurs when quantitatively accurate predictions of thresholds for critical phenomena become necessary, as seen later in §5.2.1.

Having defined the physical setup and motion equations, as well as their models in the appropriate limits, we next look at what information can be extracted in general from the parent system by simple kinematic considerations.

2.1 Kinematic effects of the boundary-interface contact

We consider the full Euler equations in a two-layer setup. The case of free surface under constant pressure can be obtained as the appropriate limiting case from this, as $\rho_{1}\to 0$. Suppose that the interface, assumed to be a graph of class $C^{1}$, initially touches the bottom at one point $x_{0}$, that is, $\eta(x_{0},0)=0$; then there is a moving contact point $x=x_{\kappa}(t)$, with $x_{\kappa}(0)=x_{0}$, at the bottom such that $\eta(x_{\kappa}(t),t)=0$ for all times $t$. It is not difficult to see that $x_{\kappa}(t)$ is the solution of the Cauchy problem

[TABLE]

where $u_{2}$ is the horizontal component of the velocity field in the lower layer. Indeed, let $Z_{1}(t)=\eta(x_{\kappa}(t),t)$. Then, using the interface velocity condition $\eta_{t}+\eta_{x}u_{2}=w_{2}$, we obtain

[TABLE]

Hence $Z_{1}(t)$ is a solution of the Cauchy problem

[TABLE]

Due to the boundary condition $w_{2}(x,0,t)=0$, the identically zero function is also a solution of the same Cauchy problem. If $(t,z)\mapsto w_{2}(x_{\kappa}(t),z,t)$ is a Lipschitz function with respect to $z$ and uniformly continuous with respect to $t$, then the solution is unique and we can conclude that $Z_{1}(t)=0$ for all $t$, as claimed.

Suppose now that $\eta(x_{0},0)=0$ for all $x_{0}$ in an interval $I_{0}$. Then we can apply the previous result to the end points of the interval (respectively for their left and right neighborhoods) as this evolves with time into $I_{t}$, and conclude that the end points never detach from the boundary. Hence, no lower layer fluid can come into contact with points within the interval $I_{t}$ as long as the interface remains a regular, single-valued function of the form $z=\eta(x,t)$. A similar argument holds for the contact of lower, denser fluid with the top boundary plate. In particular, when fluid velocities are initially zero in this case, contact is maintained during the slumping that can be expected in the evolution of the interface for as long as this remains a regular graph. Hence, to break the contact a different mechanism, one that relies on loss of regularity, is needed to detach the lower fluid from the top plate.

In order to shed some light on these issues, we now turn our attention to models of Euler fluid dynamics, which under appropriate assumptions can isolate relevant dynamics in simplified settings more amenable to analysis.

3 A single fluid over a flat surface: Airy’s model

We now consider the simplest model of a homogeneous fluid with free surface under a constant pressure, with a horizontal flat bottom placed at $z=0$. As summarized in § 2, such a model is provided by the long-wave approximation at leading order, which yields the classic Airy system (2.8). Written in non-dimensional form by the rescaling

[TABLE]

system (2.8) becomes

[TABLE]

where we have abused notation a little by not differentiating between scaled and unscaled variables; further, unless explicitly indicated, by $u=u(x,t)$ we will denote the layer averaged-horizontal velocity instead of the corresponding variable in the Euler equations, and from here on we will drop bars and layer subscripts, restoring the full notation of § 2 whenever necessary to avoid confusion. (Note that in this notation the gravity acceleration is normalized to have unit magnitude, which coincides with the convention of [35], section 2.3).

Our first aim is to show that the persistence of contact points observed above for a two-fluid system governed by the full Euler equations continues to hold for this greatly simplified model, as long as the regularity of the initial data is maintained by its evolution. We take $(\eta(x,t),u(x,t))$ to be a solution of (3.2) such that $\eta(x_{0},0)=0$. This obviously means that initially there is (at least) one dry point on the floor. Our claim is that the bottom will remain dry for some time and will not get immediately wet; more precisely, once again we have to find a (moving) point $x_{\kappa}(t)$ such that $\eta(x_{\kappa}(t),t)=0$ for all $t$ (for as long as this function is sufficiently regular before shocks form — this caveat will often be left implicit in the following). It is easy to show that the function $x_{\kappa}(t)$ we are looking for is the solution of the Cauchy problem $\dot{x}(t)=u(x(t),t)$, $x(0)=x_{0}$. Indeed,

[TABLE]

so that the function $\eta_{\kappa}(t)\equiv\eta(x_{\kappa}(t),t)$ is a solution of the Cauchy problem

[TABLE]

where the velocity gradient at the contact point, $\upsilon_{\kappa}(t)\equiv u_{x}(x_{\kappa}(t),t)$, can be viewed as a given function of time. Loss of regularity in finite times of the function $u$ typically occurs in the form of a shock, and hence leads to unbounded $\upsilon_{\kappa}$’s. Due to uniqueness, we can conclude that $\eta_{\kappa}(t)=0$ for all times before a shock forms. If the interface is tangent to the bottom somewhere, then it remains tangent, even if the tangency point can move along the bottom. Indeed, if we assume that initially $\eta_{x}(x_{0},0)=0$, we can introduce the function $\gamma_{\kappa}(t)\equiv\eta_{x}(x_{\kappa}(t),t)$, the slope of the surface at the contact point, and easily show that $(\eta_{\kappa}(t),\gamma_{\kappa}(t))$ is a solution of a Cauchy problem for a homogeneous first-order linear system, with initial data $(\eta_{\kappa}(0),\gamma_{\kappa}(0))=(0,0)$. Again from uniqueness, we can conclude that $\eta(x_{\kappa}(t),t)=0$ and $\eta_{x}(x_{\kappa}(t),t)=0$ for all times before the possible shock.

We remark that

[TABLE]

meaning that the velocity of the dry point is constant. Hence, with the additional hypothesis $u(x_{0},0)=0$ we see that $x_{\kappa}(t)=x_{0}$, that is, a dry point does not move if the initial velocity is zero at that point. (This fact was already noticed in [31].) In particular, if the velocity field of the fluid initially vanishes, we have two consequences: (i) suppose that there is only one dry point $x_{0}$, then this point does not move (even if the fluid configuration is not symmetric with respect to the line $x=x_{0}$, which would imply this as a consequence of preservation of symmetry); (ii) suppose that there is a whole interval of dry points, then these points will remain dry at all times so long as a shock does not develop. We remark that working with the “sound speed” variable (see e.g., [35], Chapter 10) $c_{s}=\sqrt{\eta}$, both consequences can be seen as corollaries of the vanishing of this speed at the contact point. (As will be seen in section 3.2, consequence (ii) is in fact false if the interface does not touch the bottom with sufficient regularity: an example is given by the initial data $\eta(x,0)=(1-x^{2})\chi_{I}(x)$, $u(x,0)=0$, where $\chi$ is the characteristic function of the interval $I=(-1,1)$.) Viceversa, if the interface does not touch the bottom at time $t=0$, then the uniqueness hypothesis for solutions of equation (3.4) in backward times shows that the interface will never reach the bottom, as long as only classical solutions of the Airy system are considered (an explicit example will be shown in section 3.2). A qualitative description of the role of shocks can be also given: suppose a smooth interface has a small-amplitude minimum at a given time $t_{0}$, say $\eta(x_{0},t_{0})=\epsilon\ll 1$ and $\eta_{x}(x_{0},t_{0})=0$. Then

[TABLE]

for as long as the velocity gradient $u_{x}(x_{0},t_{0})$ remains bounded. Thus, the argument would break down if a large velocity gradient were to develop; this would imply a fast change of the interface which could eventually make contact of the interface with the bottom possible.

It is worthwhile to add a comment on the curve $t\mapsto(x_{\kappa}(t),t)$. Since $\eta(x_{\kappa}(t),t)=0$, a posteriori we can say that it belongs to both solution families (the characteristics) of

[TABLE]

where $\lambda_{\pm}=u\pm\sqrt{\eta}$ are the eigenvalues of the matrix

[TABLE]

of the quasilinear system (3.2). Thus, another interpretation of the curve $(x_{\kappa}(t),t)$ is that the point $(\eta(x_{\kappa}(t),t),u(x_{\kappa}(t),t))$ stays on the “sonic” line $\eta=0$ where the eigenvalues coincide.

3.1 Airy’s system shocks for simple waves

When the free surface contacts the bottom, the process of shock formation for Airy’s system (3.2) solutions deserves a closer look: as we have seen, wetting can only happen with loss of smoothness, i.e., the free surface cannot detach from the bottom before a shock forms. For a class of special solutions of system (3.2), the so-called “simple waves,” we will prove that the vertical location of the shock position depends on the initial steepness of the interface near the bottom.

As well known, along curves defined by (3.7) the quantities $R_{\pm}$ (the Riemann invariants)

[TABLE]

are constant; when such constant is independent of the initial condition $x_{0}\in{\cal{I}}$ for some open interval $\cal{I}\subset\hbox{I\hskip-2.15277ptR}$, the corresponding solution $(\eta(x,t),u(x,t))$ is commonly referred to as simple wave (see, e.g,  [37], p.164). Let us focus on the $R_{-}$ invariant,

[TABLE]

and let $c$ denote its constant value. Thus, since $u(x,t)=2\sqrt{\eta(x,t)}+c$, the Cauchy problem for the Airy system reduces to

[TABLE]

We briefly review the shock formation in this case. The generic (hodograph) solution is

[TABLE]

where $F$ is the (local) inverse function of the initial condition $\eta_{0}$. The shock point $(x_{c},t_{c})$ is the lowest point, in the $(x,t)$ plane, of the shock curve

[TABLE]

where the second equation comes from ${{\mathrm{d}}x}/{{\mathrm{d}}\eta}=0$. Therefore, the parametric equations of the shock curve are

[TABLE]

Thus, the shock value $\eta_{c}$ of the field $\eta$ is the minimum of the function $t(\eta)$, where ${{\mathrm{d}}t}/{{\mathrm{d}}\eta}=0$. Hence we have

[TABLE]

and

[TABLE]

Next, we examine the shock’s position dependence on the regularity of the interface function $\eta$ at the intersection (dry) point at the bottom boundary. We consider the family of initial interfaces

[TABLE]

with $M>0$ and $\alpha>1$. We choose the constant value of the Riemann invariant $R_{-}$ to be $c=-2\sqrt{M}$, so that $u=2\sqrt{\eta}-2\sqrt{M}$ and the fluid is at rest at infinity. The hodograph solution is

[TABLE]

For these initial data, we can see that a shock occurs at finite times. Indeed, at the points $(x_{L}(t),M)$ and $(x_{R}(t),0)$, depicted in figure 3, from (3.18) it follows that ${\dot{x}}_{L}=\sqrt{M}$ and ${\dot{x}}_{R}=-\sqrt{2M}$.

The shock curve is given by

[TABLE]

When $\alpha\neq 2$ the second component $t(\eta)$ is a monotonic function of $\eta$, and the shock can only appear at the limiting values of the interface transition region,

[TABLE]

and the shock time is

[TABLE]

Thus, the interface shock formation depends on whether the exponent $\alpha$ is smaller or greater than the threshold value $\alpha=2$. For sufficiently steep interfaces the shock appears far from the bottom, while for flatter interfaces the shock appears right at the bottom at $t=0^{+}$. For $\alpha=2$, the shock time is $t_{c}=1/3$, and the shock starts as a finite amplitude jump from $\eta=0$ to $\eta=M$, since the evolution is nothing but the time-reversal version of the classical dam-break problem (see, e.g., [35]). At this critical case, higher-order terms in the expansion of the initial datum near the bottom determine the evolution. Hence, while still in the simple-wave approximation, it is instructive to study an initial interface of the form

[TABLE]

where now $\alpha>2$ and $x_{M}^{2}+\epsilon\,x_{M}^{\alpha}=M$. For sake of simplicity we let $|\epsilon|\ll 1$ and consider the leading order terms in this small parameter. For the branch that corresponds to $x<0$, the asymptotic expansion for $\epsilon\to 0$ of the inverse of the initial condition is

[TABLE]

with an error $o(\epsilon)$. The hodograph solution is therefore approximately given by

[TABLE]

so that (assuming the asymptotic relation holds for derivatives as well)

[TABLE]

or

[TABLE]

Thus, the critical time $t_{c}$ continues to be a monotonic function of the critical interface height $\eta_{c}$. The shock position depends uniquely on the sign of $\epsilon$ and not on the exponent $\alpha$ (see figure 4). If $\epsilon$ is positive the minimum of (3.26) is attained when $\eta_{c}=M$, meaning that the shock appears far from the bottom. When $\epsilon$ is negative the shock appears at the contact point with the fluid’s bottom boundary. The main difference with respect to the qualitatively similar case analyzed above with the local interface $\eta(x,0)=(-x)^{\alpha}$ for $1<\alpha<2$ is given by the critical time which is now $O(\epsilon)$-close to $t=1/3$ and does not develop immediately at $t=0^{+}$ as in the previous case.

The fact that the interface touches the bottom boundary plays a fundamental role in the shock position. To illustrate this, let us consider a vertically translated version of the local initial data (3.17),

[TABLE]

where $\delta<M$ and $1<\alpha<2$. The relation between the velocity and $\eta$ is taken as previously as in the simple-wave case, i.e., $u=2\sqrt{\eta}-2\sqrt{M}$. The hodograph solution is

[TABLE]

which implies that the time component of the shock curve is

[TABLE]

The minimum of this function is given by $\eta_{c}=\alpha\delta/(2-\alpha)$, which is strictly greater than $\delta$.

While we do not provide here a proof, numerical simulations show that the above results may apply to more general classes of initial data, whose local behavior is illustrated by the special cases considered here. In particular, in Appendix B we will provide an example of how the smoothing of the interface at $\eta=M$ does not affect the qualitative shock formation displayed in this section.

3.2 Self-similar parabolic solutions

All solutions of the Airy system (3.2) we have considered so far belong to the simple wave class, that is, one of the two Riemann invariants is constant across characteristics. This gave rise to the examples we have considered for shock formation when the free surface initially contacts the bottom boundary with $\eta(x,0)=0$ on the half-line $x\geq 0$.

Inspired by the well-known linear scaling solution of the Hopf equation $u_{t}+uu_{x}=0$ with $u=x/t$, we now study a particular class of local similarity solutions, already briefly mentioned in [32], for which the simple wave restriction can be relaxed, while maintaining continuity (in most cases) for finite times. These similarity solutions will also allow us to study setups where the free surface can come arbitrarily close to the bottom, but does not contact it. Furthermore, this class of solutions can be used as an advantageous starting point to move into two-fluid setups, and illustrate the difference between evolution from initial conditions where the interface touches the upper (as opposed to the bottom) boundary in the air-water case of two-fluid dynamics. We remark that self-similar solutions of polynomial form for free surface flows are known for the parent Euler equations themselves (as opposed to models) [24, 27, 28], albeit for fluids of infinite depth.

3.2.1 Local behavior of parabolic interfaces

We consider solutions of the Airy system evolving from the class of initial data locally defined by

[TABLE]

with $\mu_{0}>0$ and $\gamma_{0}>0$. Contact of the free surface with the bottom is clearly achieved at $x=0$ if $\mu_{0}=0$. It is immediate to see that the second order polynomial dependence of the right hand side of the Airy system (3.2) is maintained by these data, and similarly to the polynomial balance argument for the linear scaling solution of the Hopf equation $u=x/t$, whereby nonlinear terms balance with space and time derivative to yield cancellations, we make the ansatz

[TABLE]

which we assume holds locally for some finite time $0\leq t<T$. Substitution into the Airy’s system (3.2) then yields the following ODE’s for the scaling parameters,

[TABLE]

with

[TABLE]

Being nonlinear, these ODE’s will in general support a movable singular time $T>0$ depending on the initial data. Since the first equation is of Riccati type, it can be linearized and the system solved by quadratures (note that the $\nu$ and $\gamma$ equations are uncoupled from the $\mu$ evolution). The general solution in parametric form follows from

[TABLE]

where the constant $K_{1}$ and $K_{2}$ are determined by the initial conditions

[TABLE]

The minimum $\eta(0,t)=\mu(t)$ of the parabola displays the persistence of the contact point with the bottom, as the second equation in system (3.32) shows that $\mu=0$ is an invariant manifold for these ODE’s. The $\mu$-evolution by eliminating $\nu$ and $\gamma$ can be written as

[TABLE]

Thus, if $\mu_{0}>0$, $\nu_{0}\neq 0$ and the negative branch of the square root is taken so that $\dot{\mu}(0)<0$, then $\mu(t)$ is a decreasing function of time. Nonetheless, by solution uniqueness the bottom can never be reached, unless a singularity develops and Lipschitz continuity of (3.32) or (3.36) is lost. Indeed, if there exists a time $\bar{t}$ such that $\mu(\bar{t})=0$, then from (3.36) the quadrature solution leads to a divergent integral

[TABLE]

Elimination of $\mu$ and $\nu$ shows that the curvature parameter $\gamma(t)$ satisfies the differential equation

[TABLE]

The auxiliary variable $\tau=(\gamma/\gamma_{0})^{1/3}$ simplifies this evolution equation to

[TABLE]

and is useful for the parametrization of the quadrature solution. (The variable $\mu=K_{1}\gamma^{1/3}$ would play a similar role, however it is inconvenient in general because $\mu$ should be allowed to vanish, which leads to $K_{1}=0$ and thus to $\mu$ vanishing for all times.) The sign of $\dot{\gamma}$ in the solution of (3.38) is determined by the initial data through the third equation in (3.32), as $\dot{\gamma}(0)=-3\nu_{0}\gamma_{0}$. The explicit solution of (3.38) splits into different regimes, according to the sign of $\nu_{0}$ and $\gamma_{0}$.

For the purpose of comparing model predictions with the dynamics of the full Euler fluids, we now focus on initial conditions for which the fluid is initially at rest, i.e., $\nu_{0}=0$ and the parabolic surface has positive concavity, $\gamma_{0}>0$. Therefore $\dot{\gamma}(0)=0$, and the correct sign in (3.38) is chosen, when $\nu_{0}=0$, by observing that

[TABLE]

so that the curvature parameter $\gamma$ is monotonically increasing at time $t=0^{+}$. We have

[TABLE]

whose implicit solution is

[TABLE]

with respect to the auxiliary time-like variable $\tau$. As functions of this variable, the coefficients in (3.31) are, explicitly,

[TABLE]

The curvature diverges in finite time $t=t_{s}$, where

[TABLE]

At this time the velocity coefficient $\nu$ diverges to negative infinity and, if the fluid surface is initially away from the channel bottom boundary, or $\mu_{0}>0$, then the height coefficient $\mu$ diverges as well. However, if initially the surface is in contact with the channel bottom, or $\mu_{0}=0$, then the minimum of the parabolic surface will remain at zero, $\mu(t)\equiv 0$, at all times before the singularity of infinite curvature develops at $t=t_{s}$.

We remark that when $\gamma_{0}<0$, i.e., the parabola has negative concavity initially, the formulae above must be modified to account for the correct branches of the square roots. We report the appropriate, explicit relations here for completeness:

[TABLE]

3.2.2 Local parabolic solutions continuously connected to uniform depth state

The local parabolic solutions can be globally extended by “gluing” them together with a constant background $\eta(x,t)=M>0$. We will report the full details of the analysis of the time evolution of such solutions elsewhere, as these deserve a separate study to analyze all relevant scenarios. Here, we simply report results for the initial condition

[TABLE]

where $\chi_{\cal I}(x)$ is the characteristic function of the interval ${\cal I}=(a,b)\subset\hbox{I\hskip-2.15277ptR}$, that is

[TABLE]

with $H$ the Heaviside step function, and ${\cal I}^{c}$ is the complement set in IR; here ${\cal A}_{0}=(-a_{0},a_{0})$ is a symmetric interval about the origin $x=0$. If we choose $M>\mu_{0}>0$ and $a_{0}=\sqrt{(M-\mu_{0})/\gamma_{0}}$ then the initial data are continuous and piecewise differentiable. The evolution by system (3.2) preserves continuity and piecewise differentiability, as long as shocks do not form. The number of nondifferentiable points can change in time, and possibly in nature, due to their collisions. While a full analytical representation of the solution can be provided, for the present scope it is enough to simply derive equations for the motion of these nondifferentiable points, by describing them as shocks for the spatial derivatives of the dependent variables, using distributional calculus. We make the ansatz that the global behavior of the solution before a singularity time $t=t_{s}$ is

[TABLE]

where $W(x,t)$ and $V(x,t)$ are, respectively, local surface and velocity solutions of the Airy system (3.2), and $a=a(t)$, $b=b(t)$ are suitable time-dependent end points of the intervals ${\cal A}=(-a,a)$, ${\cal B}=(-b,b)$, ${\cal I}_{+}=(a,b)$, ${\cal I}_{-}=(-b,-a)$ which keep track of the nondifferentiable points of the global solution. Because spatial symmetry of solutions is preserved in time under Airy evolution, it is sufficient to follow positive half-line evolution only. Differentiating the Airy system with respect to $x$ yields

[TABLE]

When the solution ansatz (3.47) is substituted in these equations, three different kinds of distributions (Heaviside $H$, Dirac$-\delta$ and Dirac$-\delta^{\prime}$) appear and need to be separately balanced. The terms involving Heaviside functions cancel out because (3.47) is a piecewise solution of the Airy system, and the same happens with the terms involving Dirac$-\delta^{\prime}$ because of the continuity of (3.47). Therefore the only terms that give the evolution of the points $a(t)$ and $b(t)$ come from the Dirac$-\delta$’s. Using the relation

[TABLE]

yields the following set of ordinary differential equations for the evolution of the end points:

[TABLE]

for the terms supported at $x=a(t)$ and

[TABLE]

for the terms supported at $b(t)$. The (explicit) evolution of the singular points $a$ and $b$ follows from the ratio of the $W$ and $V$-field derivative jumps at the singular points,

[TABLE]

(The negative sign in the first equation is necessary since, for all times, $\upsilon>0$ and $\xi<0$.) Hence, the evolution of the location of the singular points $a$ and $b$ is simply

[TABLE]

The second equation is easily solved,

[TABLE]

The first equation requires the solution of the self-similar parabola parameters ($\nu,\mu,\gamma$). Since the auxiliary variable $\tau=(\gamma/\gamma_{0})^{1/3}$ is increasing in time (as shown in (3.42)), the $\tau$ evolution

[TABLE]

allows the parametric solution of this equation, via

[TABLE]

The evolution equation for $a$ becomes

[TABLE]

whose solution is

[TABLE]

When the point $x=a(t)$ reaches the origin it collides with its symmetric twin at $x=-a(t)$. The value $\tau_{c}$ of $\tau$ at the collision time $t_{c}$ can be obtained from (3.58) by imposing $a(\tau_{c})=0$, thus yielding

[TABLE]

From (3.42) one can explicitly obtain the collision time in terms of the initial data,

[TABLE]

In the limit $\mu_{0}\to 0$ the collision time coincides with the shock time, defined by (3.44), when the parabola collapses (closes) onto its axis,

[TABLE]

This equality is confirmed by the fact that in the same limit $\tau$ diverges and so does the curvature $\gamma$ of the parabolic part of the interface.

Two interesting quantities are the values of the interface and the velocity at $x=a$. By continuity, these are given by

[TABLE]

At the collision time, for strictly positive $\mu_{0}$, such quantities are

[TABLE]

As shown in figure 5(a) this implies that at the collision time the bottom $\mu$ has reached the background level $M$ and the parabola shrinks to a point and disappears, with the point expanding to flat region at level $M$. When $\mu_{0}=0$ the fields value at $t_{c}=t_{s}$ are

[TABLE]

As typical in these models, we also see that spatial and time limits do not commute,

[TABLE]

With $\mu_{0}>0$ the parabola disappears before a shock forms, with only two “(anti-)bumps” remaining. These bumps can be shown to be simple wave solutions moving leftward and rightward, respectively for the right and left bump. These simple waves eventually develop a standard shock, whose time of formation can be computed explicitly. When there is a dry point, i.e., $\mu_{0}=0$, the parabola persists, becoming increasingly narrow as shown in figure 5(b), to finally collapse to a vertical segment

[TABLE]

reminiscent of the schematic depicted by figure 1(a). Having lost regularity, this offers a mechanism for the “detachment” of the interface from the bottom (see figure 1(a)), since the interface at the origin has reduced to the segment (3.66), which can be eliminated to redefine the fluid domain with a new interface at the elevation $\eta(0,t_{s})=M/4$ away from the bottom. Of course, how the dynamics of the Airy system can be continued past this point is not immediately clear, as at the collision time a shock in the velocities also develops on a whole segment.

3.2.3 Compact-support self-similar solutions connecting to dry state

A variant of the parabolic self-similar solutions of the Airy’s system (3.2) which corresponds to a compactly supported fluid region is of particular interest as it illustrates, outside of the special class of simple waves and within a closed form expression, the evolution of the interface and its wet domain boundary. For this purpose, it is useful to change notation slightly to focus on the support along the floor, and rewrite the corresponding similarity equations in the new variables. Thus, we now seek solutions in the form

[TABLE]

where $\chi_{\cal B}(x)$ is the characteristic function of the interval ${\cal B}=(-\beta,\beta)$, and the coefficients $\nu$, $\beta$ and $\gamma$ depend only on time. If $\eta$ describes a parabola with negative concavity, i.e.,

[TABLE]

then the free surface $\eta(x,t)$ connects continuously to $\eta(x,t)=0$ for $x\notin{\cal B}$, while the velocity $u(x,t)$ connects through shocks at $x=\pm\beta(t)$ to the rest state $u=0$ for $x\notin{\cal B}$, for as long as $\nu(t)\neq 0$. In order for the fields (3.67) to be solutions of the Airy system (3.2), the coefficients $(\nu,\beta,\gamma)$ must evolve according to the ODE system

[TABLE]

with the three parameter family of initial conditions

[TABLE]

When $\gamma_{0}>0$ and $\beta_{0}\neq 0$ ($\beta_{0}>0$, say), the fields $\eta(x,t)$ and $u(x,t)$ always asymptote to zero at long times. This can be seen by direct quadrature of the ODE system (3.69), which again follows from the Riccati-like structure of the $\nu$ equation, and the conservation law implied by the second and third equation. It is easy to see that

[TABLE]

and so

[TABLE]

for some constant $K_{2}>2K_{1}/\beta_{0}$. The first quadrature (3.71) shows that $\gamma$ retains its initial sign for all times if $\beta$ does not go through zero. The functions $\nu(t)$ and $\gamma(t)$ are therefore determined by $\beta(t)$ through

[TABLE]

while $\beta(t)$ itself is determined implicitly by the last quadrature (3.72). This corresponds to the one-dimensional dynamics along the $\beta$-axis of a charged particle in the repulsive field generated by a same-signed charge at the origin. Thus, if $\dot{\beta}(0)<0$, which is implied by $\nu<0$ initially, $\beta(t)$ is a decreasing function of time down to the minimum

[TABLE]

attained at $t=t_{0}$, corresponding to the time when $\dot{\beta}=0$; for times $t>t_{0}$, the function $\beta(t)$ is monotonically increasing. Conversely, if $\nu_{0}>0$, then $\dot{\beta}>0$ initially and $\beta(t)$ is always monotonically growing for all times $t>0$. These properties imply $\eta(x,t)\to 0$ and $u(x,t)\to 0$ as $t\to\infty$, since $\gamma(t)$ and $\nu(t)$ are both proportional to inverse powers of $\beta(t)$. The time $t=t_{0}$ gives the maximum amplitude reached by $\eta(x,t)$, as

[TABLE]

Note that depending on the initial values $\nu_{0}$ and $\gamma_{0}$ the solution $\nu(t)$ may have a stationary point $\dot{\nu}=0$ at some intermediate time. For instance, if $\dot{\nu}(0)<0$ (i.e., $\nu_{0}^{2}-2\gamma_{0}>0$) and $\nu_{0}<0$, then $\nu^{2}$ increases until its time derivative changes sign, since in terms of the quadrature constants $\dot{\nu}=2(K_{2}-3K_{1}/\beta)/\beta^{2}$, which crosses zero when $\beta(t)=3K_{1}/K_{2}>\beta_{\rm{min}}$. This can be given a simple physical interpretation: the initial fluid velocity gradient parametrized by $\nu_{0}<0$ can be chosen in such a way that the parabolic “lump” shrinks, and if its initial peak amplitude $\gamma_{0}$ is sufficiently large the fluid will possess sufficient inertia to further steepen the velocity gradient as the lump’s amplitude grows. However, regardless of the initial sign of $\dot{\nu}$, even for the initial condition leading to initially growing lump’s amplitudes $\nu_{0}<0$, gravity eventually always dominates, causing the relaxation of the fluid velocity $u$ to zero as the peak amplitude $\eta=K_{1}/\beta_{\rm{min}}$ is attained at $x=0$ and $t=t_{0}$. From then on $\nu(t)>0$ and the fluid slumps, spreading indefinitely over ever larger intervals $x\in(-\beta,\beta)$. (Of course, this is also the only possibility for the case $\nu_{0}>0$ — i.e., $u$ with the same sign of $x$ — which leads to monotonic decrease of $\nu$ and $\gamma$ as $\beta$ grows to infinity.)

It is instructive to interpret the dynamics of this class of solutions from the perspective of the conservation laws of the Airy system (3.2). By connecting at $x=\pm\beta$ to a state with $\eta=0$, it is clear that all the balance (transport) laws

[TABLE]

whose fluxes $G[\eta,u]$ depend on $\eta$ and its derivatives effectively satisfy insulating boundary conditions at $x=\pm\beta$, which lead to conservation of the corresponding quantity $\int F[\eta,u]dx$. Thus, for instance, the mass density $F[\eta,u]=\eta$, with flux $G[\eta,u]=\eta u$, is conserved; the total mass of the fluid is given in terms of the similarity coefficients by

[TABLE]

Similarly, both the momentum

[TABLE]

(which by antisymmetry is always identically null) and the energy

[TABLE]

are conserved, as can be easily checked directly from system (3.69). It is worth remarking that energy conservation survives despite the presence of shocks in the form of the solution (3.67), unlike the generic setup for weak solutions of the Airy system, whereby conservation of mass and momentum can be enforced, but energy and other conservation laws (with higher order polynomial densities) are not conserved and decay according to the appropriate Rankine-Hugoniot conditions.

The quadrature (3.72) provides an exact functional form for $\beta(t)$ (expressed implicitly through its inverse $t(\beta)$). The detailed expression does not seem to be important, but the quadrature can be further analyzed to extract the asymptotic behavior of the solution at long times, to yield at leading order

[TABLE]

As expected, the asymptotic behavior of $\nu$ is independent of the initial data and the interface $\eta$ flattens asymptotically as $t^{-1}$, at any fixed location $x$, independently of the initial velocity.

This class of solutions can serve as a starting point of a perturbation analysis for a two-fluid system, which will be our focus from now, in the limit of vanishing density of the upper fluid, in particular for modeling the splashing process when a top plate is positioned at the maximum height (3.75).

4 Two-fluid system in a horizontal channel

We are now ready to move on to the stratified case of two fluids divided by an interface, constrained to move in a horizontal channel of height $h$ (see figure 2). As reviewed in § 2, a dispersionless model for such a physical setup is system (2.7). Such a system can be given a Hamiltonian structure, with Hamiltonian density (see, e.g., [8, 7])

[TABLE]

and Hamiltonian operator

[TABLE]

which leads to the equations of motion in the form

[TABLE]

where the functions $A(\eta,\sigma)$, $B(\eta,\sigma)$, and $C(\eta,\sigma)$ are

[TABLE]

System (4.3) becomes elliptic for large values of $|\sigma|$ (see, e.g., [15, 29] for the fluid-dynamics relevance of this feature).

The structure of the two-fluid system (4.3) is best analyzed in terms of non-dimensional variables, by introducing the parameter $r=(\rho_{2}-\rho_{1})/{\rho_{2}}$ and the rescaled variables

[TABLE]

Once again abusing notation by not distinguishing between dimensional and scaled variables, as long as this does not generate confusion, the Hamiltonian density (4.1) takes the form

[TABLE]

and equations (4.3) retain the same form with rescaled matrix elements

[TABLE]

Explicitly, in conservation form, system (4.3) (the nondimensional version of system (2.7)) reads

[TABLE]

The limit $r\to 1^{-}$ when the upper fluid’s density becomes negligible is related to the Airy system as long as the interface does not get into contact with the upper lid (see Appendix B ). The opposite limiting case $r\to 0^{+}$ gives the so-called Boussinesq approximation of negligible density effects on the fluids’ inertia. It is remarkable that in this limit the equations of motion can be mapped into the Airy’s model of a single layer fluid [16]. In fact, the map

[TABLE]

sends solutions of the $(\eta,\sigma)$-system (4.8) into solutions of the classical Airy system (3.2), with dependent variables now adorned by tildes,

[TABLE]

However, care must be paid when interpreting solutions $(\tilde{\eta},\tilde{u})$ of the Airy’s system (4.10) with respect to the original variables $(\eta,\sigma)$ since the map (4.9) is not one-to-one [20]. Thus, under the map (4.9), the same initial data $\tilde{\eta}(x,0)$ and $\tilde{u}(x,0)$ are the image, in general, of four different initial conditions in the original variables $(\eta,\sigma)$. Some consequences of this multivaluedness in the context of the boundary confinement effects will be briefly explored below in §4.3.1.

The results proved for the Airy system can be extended to the generic form (4.3) which includes all the Hamiltonian systems with the “canonical structure”

[TABLE]

and whose Hamiltonian density ${\cal H}(\eta,\sigma)$ does not depend on $x$-derivatives. In this case, the matrix entries in (4.3) have the general form

[TABLE]

The characteristic velocities of (4.3) are

[TABLE]

so that the sonic curve (i.e., the region in the $(\eta,\sigma)$-plane where $\lambda_{+}=\lambda_{-}$) has two components,

[TABLE]

Given the symmetry between these two cases, it suffices to analyze one of them, with the results easily extended to the other case. In the following we will be mainly interested in the first component of system (4.3), hence we focus on the first choice, $B(\eta,\sigma)=0$. Note that in the Airy case this amounts to $\eta=0$. We assume that $(\eta(x,t),\sigma(x,t))$ is a solution of system (4.3) such that $B(\eta(x_{0},0),\sigma(x_{0},0))=0$ for some $x_{0}$. Under suitable hypotheses, we want to show that there exists a moving point $x_{\kappa}(t)$ such that

[TABLE]

for all $t$ before shocks emerge. This can be done by considering the solution $x_{\kappa}(t)$ of the Cauchy problem

[TABLE]

Indeed, define $\beta_{\kappa}(t)=B(\eta(x_{\kappa}(t),t),\sigma(x_{\kappa}(t),t))$ and compute

[TABLE]

where \alpha_{1}(t)=-\left(B_{\eta}\sigma_{x}\right)\big{|}_{x=x_{\kappa}(t)}, \alpha_{2}(t)=-\left(B_{\sigma}C\right)\big{|}_{x=x_{\kappa}(t)}, and $\gamma_{\kappa}(t)\equiv\eta_{x}(x_{\kappa}(t),t)$. If the hypothesis $B_{\sigma}=0$ is satisfied (i.e., if $B$ depends only on the variable $\eta$), then the previous equation is a (linear) homogeneous equation for $\beta_{\kappa}(t)$. Since $\beta_{\kappa}(0)=0$, we obtain our thesis.

Now let us remove the hypothesis $B_{\sigma}=0$ and let us suppose that $\eta_{x}(x_{0},0)=0$. From the definition of $\gamma_{\kappa}(t)$, it is easily found that

[TABLE]

where $\alpha_{3}(t)=-\sigma_{xx}(x_{\kappa}(t),t)$, \alpha_{4}(t)=-\left((A_{\sigma}+B_{\eta})\sigma_{x}\right)\big{|}_{x=x_{\kappa}(t)}, and \alpha_{5}(t)=-A_{\eta}\big{|}_{x=x_{\kappa}(t)}. By uniqueness of solutions of system (4.17-4.18) we can conclude that

[TABLE]

As a consequence, one can easily check that $\frac{{\mathrm{d}}}{{\mathrm{d}}t}A(\eta(x_{\kappa}(t),t),\sigma(x_{\kappa}(t),t))=0$, so that the “contact point” $x_{\kappa}(t)$ moves at the constant speed

[TABLE]

In particular, if $A(\eta(x_{0},0),\sigma(x_{0},0))=0$ then $x_{\kappa}(t)=x_{0}$, so that $B(\eta(x_{0},t),\sigma(x_{0},t))=0$ and $\eta_{x}(x_{0},t)=0$ for all $t$.

Analogous conclusions can be drawn for the case $C=0$, under either of the conditions $C_{\eta}=0$ or $\sigma_{x}(x_{0},0)=0$.

We can apply these general results to the two-fluid system (4.8), where $B$ does not depend on $\sigma$. Hence, if $(\eta(x,t),\sigma(x,t))$ is a solution such that $B=0$ at $x=x_{0}$ and $t=0$, then, for as long as this solution maintains regularity, $B(\eta(x_{\kappa}(t),t),\sigma(x_{\kappa}(t),t))=0$, where $x=x_{\kappa}(t)$ is the solution of the Cauchy problem

[TABLE]

Since $B=0$ if and only if $\eta=0$ or $\eta=1$, we have the following conclusions, for as long as the solution $(\eta(x,t),\sigma(x,t))$ exists and is regular: (i) if the upper (lower) fluid initially touches the bottom (top) plate at a point $x_{0}$, then it will stay in contact with the bottom (top) at the moving point $x_{\kappa}(t)$; (ii) $\eta_{x}(x_{0},0)=0$ implies that $\eta_{x}(x_{\kappa}(t),t)=0$ and that $\dot{x}_{\kappa}(t)=A(\eta(x_{\kappa}(t),t),\sigma(x_{\kappa}(t),t))$ is constant; (iii) if the velocity field of the fluid is zero initially, then $\sigma(x_{0},0)=0$ and $\dot{x}_{\kappa}(t)=\dot{x}_{\kappa}(0)=0$, so that $x_{\kappa}(t)=x_{0}$, and the remarkable behavior noticed in points (i) and (ii) after equation (3.5) above can be observed in this setting as well (however, note that, as we have seen in the previous section, the contact point can certainly move if the interface has a discontinuous derivative at the contact point).

4.1 Behavior of simple waves near hyperbolic-elliptic transition curves

We recall that an alternative, more general definition of a simple wave for system (4.3) is a particular solution of the form

[TABLE]

According to this definition, we will refer to the map $\theta\mapsto(N(\theta),S(\theta))$ in the $(\eta,\sigma)$-plane as the simple wave curve. It can be easily checked that (4.22) is a solution of (4.3) if

[TABLE]

that is, the simple wave curves are solutions of the equation

[TABLE]

In order to give a geometrical interpretation of the condition $B_{\sigma}=0$ encountered in the previous section, we consider the curve $B(\eta,\sigma)=0$. We observe that equation (4.24) implies that any simple wave curve has a horizontal tangent at a point $(\eta_{0},\sigma_{0})$ where $B$ vanishes. In the generic case $B_{\sigma}\neq 0$, the behavior of simple waves near to the sonic line is [25]

[TABLE]

where the subscript 0 stands for evaluation at $(\eta_{0},\sigma_{0})$. If additionally this curve has vanishing derivative $B_{\sigma}(\eta_{0},\sigma_{0})=0$, then the tangent to the curve $B(\eta,\sigma)=0$ at $(\eta_{0},\sigma_{0})$ is also horizontal, so that the simple wave curve (4.22) and the curve defined implicitly by $B=0$ are tangent at $(\eta_{0},\sigma_{0})$.

The behavior of a simple wave near a generic point of intersection with the curve $B(\eta,\sigma)=0$, under the assumption $B_{\sigma}=0$ at this point, can now be analyzed. The asymptotic expansion for $\sigma\to\sigma_{0}$

[TABLE]

in equation (4.24) yields at the leading order

[TABLE]

Since by assumption $B_{0}=\left(B_{\sigma}\right)_{0}=0$, this asymptotic relation is satisfied for $\alpha=2$ and $\Omega=\left(B_{\eta}\right)_{0}/(4C_{0})$, so that the generic behavior of simple waves near this class of intersection points is parabolic, independently of the explicit system (4.3):

[TABLE]

The value $\sigma_{0}$ is by construction approximately constant along the simple wave curve (4.28). We now recall that the Riemann invariants $R_{\pm}$ are the transported quantities along the solutions of (4.3) with characteristic velocities $\lambda_{\pm}$, i.e., $\left(R_{\pm}\right)_{t}+\lambda_{\pm}\left(R_{\pm}\right)_{x}=0$. In the $(\eta,\sigma)$-plane, this means that the level curves of $R_{+}$ or $R_{-}$ are simple wave curves. This implies that $\sigma_{0}$ gives the local behavior of the Riemann invariants.

These results should be compared with the Airy’s case, analyzed in section 3, where the previous computations are exact. Indeed, we have $(\eta_{0},\sigma_{0})=(0,u_{0})$, $B(\eta,u)=\eta$ and $C(\eta,u)=1$, so that (4.28) takes the form

[TABLE]

Hence, the Airy case gives the local generic behavior of the more general system (4.3).

We remark that for simple wave curves getting into contact with the sonic lines, generically the matrix in system (4.3) will not be diagonalizable (being of Jordan form with an off-diagonal nonzero entry and corresponding eigenvalue of geometric multiplicity one). The special case $B=C=0$, i.e, when the matrix is proportional to the identity, is of interest. Simple wave curves behave linearly in a neighborhood of these points, or $\eta\propto\sigma$ (for generic analytic functions $B$ and $C$). For our systems, this occurs at the corners of the hyperbolicity domain (discussed in Appendix A, figure 21), which coincide with the locations of the physical domain boundary. A general study, when these points occur away from confining boundaries, is possibly more relevant, and a substantial step in this direction can be found in [17]. We leave this to future studies, and move on to examining the shock structure at the fluid domain boundaries next.

4.2 Two-fluid system’s shocks for simple-waves

We want to extend the study of §3.1 for the Airy system of solution behaviors near the bottom of the channel at the shock time, by adapting it to the results of §4.1.

As seen in that section, simple waves are found by solving the differential equation

[TABLE]

whose general solutions are not explicitly known in the two-fluid case. However, it is sufficient for the present purposes to consider approximations close to the channel boundary $\eta=0$, the other case $\eta=1$ being totally analogous. We look for a solution in a neighborhood of the point $(0,\sigma_{0})$: from (4.28) we have that

[TABLE]

see figure 6 for a comparison with the Airy case.

Moreover, from $\left(R_{+}\right)_{t}+\lambda_{+}\left(R_{+}\right)_{x}=0$ it follows that $\eta(\sigma(x,t))$ is a trasported quantity,

[TABLE]

The characteristic velocity $\lambda_{+}(\eta,\sigma(\eta))$ can be approximated by (4.31),

[TABLE]

so that an approximation of (4.32) is given by

[TABLE]

After rescaling $\eta\to\left(1-(1-r){\sigma_{0}}^{2}\right)\eta$, this reduces to the same equation as in (3.11). Therefore, the results of section 3.1 apply to the two-fluid system as well. In particular, the degree of regularity of the initial interface in contact with the bottom affects the shock position, and the critical exponent is still $\alpha=2$ (both for the initial conditions (3.17) and (3.22)).

Higher order terms in the approximation of the characteristic velocity $\lambda_{+}$ could correct this behavior. However, at the next order in $\eta^{1/2}$ (see (4.31)), the dependence of $\sigma$ on $\eta$ is given by

[TABLE]

and the characteristic velocity $\lambda_{+}$ in (4.33) becomes

[TABLE]

If we consider the initial condition (3.17), the study of the shock curve associated with the hodograph solution

[TABLE]

leads to the same qualitative behavior as in the single layer Airy’s case.

4.3 Self-similar solutions for two-layer fluids: the high and low density contrasts

The study of the self-similar solutions of the Airy system in section 3.2 is directly relevant to the case of two-fluid flows in two (opposite) limits mentioned in section 4: when $r\to 0$, i.e., the densities of the two fluids are close, the Boussinesq approximation can be invoked, so that under mapping (4.9) the two-layer system reduces to the Airy’s model; in contrast, when $r\to 1^{-}$, i.e., the upper layer fluid density becomes negligible, the two-fluid system becomes effectively single layer and the Airy’s model applies. However, asymptotic corrections to these limits introduce subtle issues, which we will briefly examine next.

4.3.1 The Boussinesq limit $\rho_{1}\sim\rho_{2}$

As mentioned at the beginning of this section, the Airy’s system can play a significant role for its two-fluid counterpart in the Boussinesq limit, as solutions of (4.10) under the inverse of map (4.9) yield solutions of system (4.8) in this limit.

Explicitly, the (multivalued, one-to-four) inverse of (4.9) is

[TABLE]

where

[TABLE]

In particular, under this map the class of self-similar solutions (3.31) yields the local behavior of solutions of the two-fluid system (4.8) near the bottom or top plate. For instance, when the interface $\eta(x,t)$ is in contact with the top plate, or $\eta(0,t)=1$ for some time interval $t\in[0,t_{s})$, the map is, explicitly,

[TABLE]

where the coefficients $\gamma$ and $\nu$ evolve in time according to (3.32), with $\mu_{0}=\mu(t)=0$ until the global shock time $t_{s}$ (3.61). The evolution of the local interface $\eta$ (for this class of initial data, this coincides at $t=0$ with half an ellipse with vertical and horizontal axes) and shear $\sigma$ given by this solution is depicted in figure 7. These local solutions can be connected to a background constant state $\eta=M\geq 1/2$ as in § 3.2.2. The evolution of the Boussinesq approximate solutions thus obtained can be compared with a direct numerical integration of the two-fluid system (4.8), when the densities are sufficiently close, with initial conditions in the class (4.40). This will be taken on again below in section 5.

4.3.2 “Air-water” limit $\rho_{1}\to 0$

Solutions of the Airy’s system in the form (3.67), studied in § 3.2.3, can be useful to investigate the case of a two-fluid system in a channel in the limit of vanishing density of the upper fluid. For this, it is convenient to go back to the original formulation in dimensional form of the two-fluid system in long-wave asymptotics (2.3), sketched in section 2, which retains the interfacial pressure $P(x,t)$. Thus, $\eta_{1}(x,t)$ and $u_{1}(x,t)$ (respectively, $\eta_{2}(x,t)$ and $u_{2}(x,t)$) denote the thickness and the vertically averaged horizontal velocity of the upper (respectively, lower) fluid at $x$, with $\rho_{1}$, $\rho_{2}$ their constant densities. In the hydrostatic, nondispersive approximation, the equations of motion are obtained by dropping the $D_{j}$ terms

[TABLE]

If $h$ is the height of the channel and the velocities vanish at infinity, we have the constraints

[TABLE]

the second one expressing volume conservation from incompressibility. These constraints link the upper fluid thickness $\eta_{1}$ and velocity $u_{1}$ to the ones ($\eta_{2}$ and $u_{2}$) of the lower fluid by

[TABLE]

Use of the density weighted shear $\sigma(x,t)$ is suggested by the elimination of the interfacial pressure gradient $P_{x}$, and in fact this can be seen as the counterpart, for two-layer stratification, of the well-known step of eliminating the pressure gradient by taking the curl of the momentum equation in (2.1). Just as in that case, if a solution $(\eta,\sigma)$ of system (2.7) is known, the interfacial pressure $P$ can then be found, in this specific case through relations (2.6), (4.41) and (4.42). In particular, the layer mean-velocities can be expressed solely in terms of $\eta$ and $\sigma$, e.g.,

[TABLE]

so that, in turn, the interfacial pressure gradient can be expressed in terms of $\eta$ and $\sigma$ through either one of the momentum equations (4.41).

The above relations allow to set up an asymptotic expansion of the solutions to system (4.41), or, equivalently, system (2.7), in the ‘air-water’ limit, that is, in the small parameter $\rho_{r}\equiv 1-r=\rho_{1}/\rho_{2}\ll 1$,

[TABLE]

starting from solutions of the Airy’s system (2.8). In fact, at leading order the upper layer momentum equation yields $P_{x}^{(0)}=0$ (i.e., the interfacial pressure is constant, as in the usual free-surface water wave problem), so that at this order the system for governing the evolution of $\eta^{(0)}(x,t)$ and $u_{2}^{(0)}(x,t)$ coincides with Airy’s. At the next order, we have

[TABLE]

where the leading order upper layer mean horizontal velocity $u_{1}^{(0)}$ is related to the leading order solution of the Airy’s equation $u_{2}^{(0)}$ through (4.43)

[TABLE]

Substituting this into (4.46) and using the Airy system for the time derivatives yields an expression for the interfacial pressure gradient in terms of $\eta^{(0)}$, $u_{2}^{(0)}$ and their $x$-derivatives,

[TABLE]

The next order correction to the leading Airy evolution is then obtained by solving the linearized (inhomogeneous) Airy system

[TABLE]

with zero initial data. An analytic expression for the solution of this equation can in principle be computed by characteristics (which coincide with those of the leading order Airy system’s solution), in particular for those emanating from the class of initial conditions (3.46). We will not attempt this here, however some information can be gleaned from the knowledge of the interfacial pressure gradient $P_{x}^{(1)}$ when such leading order solutions are used. For instance, when the initial data are such that the interface approaches the upper plate, the denominators in (4.48) vanish, signaling the possibility of a divergence of the interfacial pressure gradient (which can lead to the breakdown of the asymptotic expansion for $\rho_{r}\to 0$, by driving the first order corrections $\eta^{(1)}$ and $u_{2}^{(1)}$ to become of magnitude comparable to $O(1/{\rho_{r}})$). In contrast, no such mechanism for divergence exists when the interface contacts the lower plate.

More specifically, when the class of self-similar solutions (3.31) is used for the leading order fields, that is $\eta^{(0)}=\mu(t)+\gamma(t)x^{2}$ and $u_{2}^{(0)}=\nu(t)x$, the asymptotic approximation of the upper fluid velocity as $\rho_{r}\to 0$ evaluates to, from (4.43) and (4.47),

[TABLE]

where $\tau$ is the function of time defined implicitly by (3.42). The slope of the upper fluid velocity $u_{1}(x,t)$ at $x=0$ is therefore, at leading order,

[TABLE]

This expression has an overall minimum at some time $t>0$, $t_{m}$ say, corresponding to

[TABLE]

since $\tau(t)$ is a monotonic (decreasing) function of time (with $\tau(0)=1$ and $\tau(\infty)=0$). Thus, as $\mu_{0}\to h$ the amplitude of this minimum grows unbounded, while the corresponding time $t_{m}(\mu_{0})$ tends to zero, see figure 8. Of course, as noted above, the asymptotics may become invalid as $\eta^{(0)}(0,t)\to h$, and it is necessary to check whether the $O(\rho_{r})$ terms in (4.50) are negligible. However, the coefficient of $\rho_{r}$ in this expression multiplies the next order field $\eta^{(1)}$ and $u_{2}^{(1)}$, which are zero initially, thereby providing a way to preserve the asymptotic ordering in (4.50) initially in time. This is notably different from the case of the parabola touching the bottom plate, where, e.g., the slope of the velocity at the origin, $\partial_{x}u_{2}^{(0)}(0,t)=\nu(t)$ diverges only after a finite time $t=t_{s}$ has elapsed.

4.3.3 Critical depth-ratio and short time evolution

We now return to system (4.3), and consider a parabolic initial interface (the same as that studied for the Airy equation in §3.2) with vanishing initial velocities. More precisely, for this section we go back to dimensional variables and assume that

[TABLE]

We focus on the evolution of the concavity at $x=0$, with the aim of isolating, in the time dependence generated by such initial data, those critical phenomena that can be used for comparisons with numerical simulations of both the dispersionless system and the stratified incompressible Euler equations. We stress that the choice of initial conditions (4.53) is purely for convenience of comparison with the limiting cases we studied, where the evolution can be found analytically. The lack of explicit expressions for characteristics and Riemann invariants prevents the analog of such calculations to be carried out for system (4.3). We also stress that what follows is generic for any initial data with sufficient regularity and symmetry around the location $x=0$.

In fact, because of the symmetry of the initial interface, we have that $\eta_{x}(0,t)=0$ for all $t$. Using the equation of motion (4.3-4.7), yields $\eta_{t}(0,0)=0$ and

[TABLE]

We first remark that $\eta_{tt}(0,0)=0$ when $\mu_{0}=h$ (i.e., if the interface initially is in contact with the top of the channel). This is in agreement with the general results obtained in the previous sections, that is, $\eta(0,t)=h$ for all times $t$ as long as the solution maintains regularity. Moreover, and again in agreement with the previous findings, if $r=1$, that is, when $\rho_{1}=0$, then $\eta_{tt}(0,0)=-2gh\mu_{0}{\gamma_{0}}<0$, indicating that the interface detaches from the upper lid at $t=0^{+}$.

Continuing with the evolution of the concavity, it is easy to show that $\eta_{xxt}(0,0)=0$ and that

[TABLE]

This is negative when the initial height $\mu_{0}$ is such that

[TABLE]

Hence, when $r<1$ is fixed, there is a critical height $\eta(0,0)=\eta_{\rm cr}\equiv{\mu_{0}}_{\rm cr}$ given by the right-hand-side of (4.56) beyond which the curvature of the interface is initially growing at the maximum point $x=0$, even though the height $\eta(0,t)$ might be decreasing. Conversely, when $r=1$ the interface flattens for all values of $\mu_{0}$, in agreement with the results of §3.2. It is interesting to study how the dispersion effects neglected by system (4.3) affect these conclusions on the initial evolution of this class of data. While the calculation can be laborious, for the class of initial data under consideration the details can be worked out by going back to the dispersive system (2.3), and it is not too difficult to show that the threshold initial height separating the two curvature behaviors is determined by the modified curvature acceleration

[TABLE]

where the dispersion contribution $D$ is

[TABLE]

where the coefficients $D_{n}$, $n=1,\dots,4$, are

[TABLE]

The threshold initial height $\mu_{0}$ is obtained by setting the acceleration (4.57) to zero. To be asymptotically consistent with the remainder of the long-wave expansion for model (2.3), which corresponds to the limit $\gamma_{0}\to 0$, the zeros of the resulting polynomial expression must be computed in this limit, which yields

[TABLE]

While not asymptotically consistent, it is nonetheless remarkable that by following the root ${\mu_{0}}_{\rm cr,D}(\gamma_{0})$ of the full right-hand-side of (4.57) (that is, with only a subset of the higher order terms in the small parameter $\gamma_{0}h$ accounted for in the dependence of ${\mu_{0}}_{\rm cr,D}$ on $\gamma_{0}$), the qualitative dependence of the resulting critical height on the curvature $\gamma_{0}$, and in particular the nearly constant shift upward for the “large” curvature parameter range $0.2<\gamma_{0}h<1$, is reflected by the dynamics of parent Euler system’s evolution, as we shall see next.

As a final remark before leaving this section, notice that sending $r\to 1^{-}$ and then $\mu_{0}/h\to 1^{-}$ in (4.55) yields (recall $|\gamma_{0}h|\ll 1$)

[TABLE]

whereas exchanging the two limits gives

[TABLE]

This phenomenon is yet another instance, as already seen in § 4.3.2, of the fact that when the interface gets closer and closer to the channel boundaries the air-water limit $r\to 1^{-}$ is not in agreement with the single Airy system’s behavior.

5 Numerical simulations

The extent to which the two-fluid non-dispersive equations we have been studying so far, together with their density contrast limiting cases, manage to capture relevant elements of the dynamics of a real, stratified fluid is of course a question that asymptotics at a formal level cannot answer directly. To some extent, this question can be addressed by accurate numerical simulations of the model and parent Euler equations, and we now turn our attention to how predictions from asymptotic results are verified in this context. We choose to work in the context of specific examples where the initial velocity is taken to be zero, which eliminates the approximation errors associated with the reconstruction of the 2D velocity field from the layer-averaged horizontal component. The density field for the Euler simulations will be smooth, albeit with a sharp transition between two essentially homogeneous densities. We focus on critical phenomena, whereby the choice of initial data will lead to different evolutions that can be predicted by the models and discerned in the Euler direct numerical simulations.

5.1 Two-layer dispersionless model

Continuous initial data with zero velocity, i.e. $\sigma=0$, and $M<\eta(x,0)\leq 1$ (restoring nondimensional variables), with the constant $M>0$ chosen to be sufficiently large so as to avoid hyperbolic-elliptic transition for $t>0$ in the evolution governed by system (4.3), can be efficiently integrated numerically with the method of lines until shock forms. We use the version of this algorithm as implemented by the “NDSolve” command environment of Mathematica 8 software, with convergence tested against closed form solutions (when available) as well as by monitoring errors with the “MaxStepSize” option ranging from $0.01$ (low resolution) to $0.0005$ (the maximum resolution for reasonable computation times – in the order of tens of minutes – on a MPB Air laptop). We choose to work with the (nondimensional) class of initial data (4.53), i.e., piecewise initial conditions formed by arcs of parabolae centered at $x=0$ intersecting the constant elevation $\eta=M$ away from the origin. As the evolution out of these initial conditions bears strong similarities with their counterpart for Airy equations (3.2), which can be computed analytically, it can serve as term of comparison to highlight the essential differences induced by the two layer model with generic density contrast.

5.1.1 Positive interfacial curvature initial data: upper fluid extending into the lower one

We first look at the parabolae with positive curvature, i.e., with the upper layer fluid closer to bottom, in the class

[TABLE]

As the motion equation in divergent form (2.7) shows, the evolution for this class of data is invariant under the scaling $x\to x/\sqrt{\gamma_{0}}$, $t\to t/\sqrt{\gamma_{0}}$, so that different initial curvatures can simply be related to one single choice of $\gamma_{0}$. The densities will be taken to be $\rho_{1}=0.111$, $\rho_{2}=0.999$, corresponding to $r\simeq 0.888$, and the background constant $M=0.5$ throughout. As in the Airy’s case in §3.2, generically when $\mu_{0}>0$ the elevation of the interface $\eta(0,t)$ grows till it recovers the background level $\eta(0,t)=M$ at some time $t=t_{M}$, after which the interface correspond to a flat, growing interval around $x=0$ where $\eta=M$, sandwiched between two symmetric simple waves (or “wings”). Figure 9 shows the initial conditions $\eta(x,0)$ and snapshots of their evolution $\eta(x,t)$ and $\sigma(x,t)$ at later times for the case $\mu_{0}=0.1$ and $\mu_{0}=0$, where the interface is close or contacts the bottom at $x=0$, respectively. Unlike their Airy’s counterpart, for the former case the concavity of the rising interface at $x=0$ changes at some intermediate time before the flat interface region develops, while for the latter case the interface remains attached to the bottom with local positive curvature, similarly to the interface behaviour for the Airy system counterpart. In both cases, shocks develop at the “back” of the simple waves joining the local interface solution near $x=0$ at points where the derivative jumps.

5.1.2 Negative interfacial curvature initial data: lower fluid extending into the upper one

We now consider cases with negative initial curvature, where the lower layer, denser fluid approaches the top lid. We first look at an example where the Boussinesq approximation holds, and use this as a test for comparison with the exact solutions obtained by the map (4.38). Figure 10 shows snapshots of the evolution from initial data obtained through this map from $\tilde{\eta}(x,0)$ and $\tilde{u}(x,0)$ in the class (3.46) (corresponding to the choice depicted in figure 5 (b)). Here, to further test the limitations of the Boussinesq approximation, the density of the upper fluid is chosen to be $80\%$ that of the lower fluid, or $r=0.2$. Note that at time $t=3.0$ the map of the exact solution (3.31) with respect to the variables ($\tilde{\eta},\tilde{u}$) provides a close approximation to the “core” (parabola’s support in $(\tilde{\eta},\tilde{u})$ variables) region of the two-layer model, even with this relatively large departure from the Boussinesq limit $r=0$.

Even though the connection with the simple-wave wings may involve shocks, i.e., vertical (infinite) slopes at the connection points, the numerical solution tracks the analytical prediction until fairly close to the time $t_{s}$ of collapse of the parabola (not shown here), and this time predicts the upper bound for the time of singularity formation detected numerically.

Next, we focus on examples from the air-water limit counterpart. Unlike the Boussinesq limit, where there is up and down symmetry and the dynamics ensuing from the interface approach to the top lid mirrors the one to the bottom, the negative curvature cases can be tested for the critical phenomena discussed in §4.3.3. For the density ratio parameter $r=0.89$, the critical elevation predicted by the dispersionless model is ${\mu_{0}}_{\rm cr}=0.75$. We illustrate this by bracketing this value with different choices of $\mu_{0}$ and tracking the time history of curvature in figure 11.

Simulations of the dispersionless system (4.3) further indicate that this critical phenomenon appears to manifest itself in a different form at later times in the evolution out of initial data in the class (4.53). In fact, when the apex of the parabola is sufficiently close to the top lid, the curvature at $x=0$ changes sign before the bump completely flattens, which makes a “dimple” appear at the top of the bump during the slumping phase, while for initial elevations below the critical height ${\mu_{0}}_{\rm cr}$ the curvature maintains its (negative) sign throughout the slumping at $x=0$. This is illustrated by figure 12, where two snapshots at the same nondimensional time of the evolution from initial conditions $\mu_{0}=0.7$ and $\mu_{0}=0.95$ are shown side-by-side. We will see that such critical phenomena are capturing similar features in the evolution of the full Euler problem with continuous, nearly two-layer stratification.

The last computation before leaving this section focuses on the extreme case of initial conditions (4.53) with negative curvature, whereby the interface is in contact with the upper plate at the parabola apex $\eta(0,0)=1$. Figure 13 depicts the initial condition and the snapshots of the interface $\eta(x,t)$ and shear $\sigma(x,t)$ shortly before a shock develops in the latter field. The contact with the upper plate is not lost until this time, but, more importantly, the way the solution loses regularity is markedly different than the case when the contact point is at the bottom plate. The shock in $\sigma$ develops while the interface region over the lower fluid is still very broad with bounded slope; at the contact point there appears to be a change of concavity, with inflection points close to $x=0$ and a jump in derivative developing, in general agreement with the perturbative analysis of §4.3.2.

5.2 Full Euler evolution from zero-velocity initial data

We now present some numerical simulations of the full stratified, incompressible Euler equations (2.1) performed using the numerical algorithm VarDen [1], with the aim of checking whether the critical phenomena isolated above can actually be detected in the evolution of the full Euler system. We typically use a square grid with $1024$ initial points along the vertical, although we have run cases with higher (up to 2048 points) and lower (down to 128 points) resolution to assess convergence. We also used 6 levels for the adjustable mesh refinement limit. Lateral boundary conditions will be assumed to be periodic, while impermeable top and bottom plates (zero-vertical velocity) will be assumed. The code is optimized for simulations of stratified fluids and is second-order in both space and time using adaptive mesh refinement (see [1] for further details about the algorithm). All the numerical integrations were performed on the University of North Carolina’s KillDevil computing cluster.

For the simulations reported here, we again choose to work with a fairly high density contrast matching the one used for the dispersionless model, i.e., $\rho_{\text{min}}=.111$ and $\rho_{\text{max}}=.999$. Likewise, the magnitude of the acceleration of gravity is scaled to be unity. The initial conditions we consider are a close approximation to the class (4.53) for the two-layer limit, that is, the velocity is taken to be initially zero and the stratification is determined through the function

[TABLE]

(where the choice $\max$ or $\min$ is for negative or positive $\gamma_{0}$, respectively) which tracks the inflection point of a continuous density distribution

[TABLE]

To approximate the two layer limit continuously, we typically work with the scale $\lambda=1/100$. The evolution of the second derivative of the mean-density isoline at the center of the domain is computed by fitting a parabola centered at $x=0$, thus tracking the value of the best-fit second derivative over time. In all the figures reporting snapshots of direct numerical simulations, we limit ourselves to portraits of instantaneous density field, color coded so that red corresponds to $\rho_{2}=0.999$ transitioning to blue at $\rho_{1}=0.111$. The computation domain is kept to rectangular box with aspect ratio $32\times 1$, with only the intermediate window $x\in[-8,8]$ visualized in the figures. Axis and color map are reported in the last figure, 20, where the overlay of the Euler density field with the interface position from the corresponding two-layer model at comparable times in the evolution is shown.

5.2.1 Critical height and curvature evolution

As a first test of the Euler dynamics we look at the phenomenon identified by the two-layer model depicted in figure 12 for the curvature evolution. Figures 14 and 15 depict the initial density configuration, together with a snapshot of the density evolution at time $t=2.1$. The cases considered correspond to the initial height $\mu_{0}=0.7$, figure 14 and $\mu_{0}=0.95$, figure 15, respectively. The initial curvature is kept constant at $\gamma_{0}=1/16$ for both cases. The formation of the “dimple” for the higher initial elevation initial condition can be seen in the mean position of the density stratification, while the lower elevation case maintains the negative curvature sign around the origin. The general shape of the stratification is qualitatively similar to that of the dispersion model evolution at the same time depicted in figure 12, with the simple wave “wings” being represented at this time by the full Euler dynamics, although of course some amount of smoothing by the dispersion, plus the effects of shears developing with the smooth stratification alter the fine scale structure of the solution with respect to the hyperbolic model. It is worth remarking that the full Euler system admits an upper bound on internal wave speeds for the stratification considered herein; once this upper bound is estimated, it is easy to check that no internal waves generated, for instance, by the discontinuous derivative of the density isolines at the parabola shoulders can reach the region around $x=0$ at the time of these snapshots, which further compels to attribute this behavior to the dispersionless response of the system.

Next, we test how the Euler dynamics reflects the critical height of § 4.3.3 for the curvature evolution from initial conditions in the class (5.3). Figure 16 shows the time history of curvature for two different choices of the curvature parameter $\gamma_{0}$, $|\gamma_{0}|=1/256$ and $|\gamma_{0}|=1/16$ and three values of the initial heights $\mu_{0}$ which bracket the threshold elevation detected by monitoring how the initial curvature evolution switches from monotonic decreasing to increasing in time. Note the similarity, even at a quantitative level, with the hyperbolic two-layer model predictions of figure 11.

The data on the critical height with respect to height and curvature parameters $\mu_{0}$ and $\gamma_{0}$ from the simulations of the full Euler system are summarized in figure 17. We consider an array of initial conditions, with heights $\mu_{0}$ and curvature parameters $\gamma_{0}$ bracketing the theoretical values identified in section 4.3.3. As can be seen from figure 17, even when applied outside of its strictly asymptotic validity, the critical curve determined by equation (4.57) (lying outside and roughly paralleling the convex hull of the circle symbols which mark the curvature-increase parametric region) can provide a fair estimate of the critical height even for large initial curvature parameter $|\gamma_{0}|\simeq 1$. Of course, a more accurate quantitative agreement is beyond reach of this simple approximation: effects such as higher order dispersion in the model asymptotics as well as those originating from the full Euler system with smooth stratification must play a role; with this in mind, it is remarkable that the prediction is verified qualitatively at all. We also stress that the relative steepness of the linear, asymptotically consistent estimate of the relation between $\mu_{\rm cr,D}$ and $\gamma_{0}$ as $\gamma_{0}\to 0$ makes for a relatively sensitive dependence of the critical height on the curvature parameter and long-wave asymptotics, which can be seen with the log-scale of the inset of figure 17. This shows that dispersive effects can become relatively important at short times, even before gradient catastrophes have begun to form, in this case in the smooth quadratic region near the parabolic apex. This early-time effects of dispersion for smooth initial data are reminiscent of those quantified by pressure imbalances in stratified fluids investigated in [4, 5]. We remark that each grid point in the two dimensional parametric search summarized by figure 17 represents a relatively expensive numerical computation Euler equations, so that further refinement of the parametric search would require a more intensive investment of computation resources that is best left to a separate study.

5.2.2 Singular behavior at the bottom and top plates

In the same vein as the previous numerical simulations, it is interesting to see how the full Euler system evolves initial conditions for which the initial mean density isoline contacts the bottom or top plate, and observe the distinct behavior in the evolution of initial configurations brought forth by the two-layer model. In fact, as seen in figure 9(b), the bottom contact case evolves similarly to the Airy’s solution, with the interface remaining in touch with the plate, and the curvature remaining finite until a ‘global’ shock develops in the interface profile at the contact point. In contrast, for the case of the interface touching the top plate, while the persistence of contact until loss of regularity is also observed, this occurs with the formation of a corner (jump in derivative) in the interface profile at the contact point (accompanied by a standard shock – derivative going to infinity at an inflection point – at the origin for the corresponding velocity shear field $\sigma$), see figure 13.

These trends can in fact be detected by the full Euler simulations. Figure 18 shows snapshots of the evolution of the density field for the bottom contact case. As can be seen, density isolines remain attached to the bottom plate, with curvature staying finite for the isolines near the contact region, for a finite time of the evolution, with the interfacial region bearing strong similarities with the two-layer counterpart. The Euler simulation is not limited by the shock formation of its two-layer counterpart, and offers a glimpse of how the evolution continues past the time of loss of regularity with the formation of an upward moving jet of higher density fluid.

Snapshots of the density field evolution for the case of top contact are depicted in figure 19. The loss of regularity occurs much earlier than the bottom touching case. Density isolines remain attached to the top lid for times past that of the last panel in the figure, $t=1.7$, but eventually lack of resolution introduces artificial effects that can not truly represent the PDE solution of the Euler system. Shear instabilities may also develop at later times in the course of the evolution, however the “stem” region attaching the lower layer fluid to the top plate seemingly persists up to those times.

6 Discussion and conclusions

In this work we investigated the dynamic implications of interactions of stratified fluids with boundaries. These can be either of confining, rigid-wall type or of internal nature such as interfacial density jumps. We first looked at simple kinematic information that can be extracted from the full Euler equations, both for homogeneous fluids where a free surface contacts a bottom in “dry spot” configuration, and for a two-layer stratified fluid whose interface contacts either the bottom or top plate. In particular, we clarified in this context the interpretation of the typical characterization of boundary points as remaining on the boundary for all times. Next, working with limiting models, such as those that derive from two-fluid setups by nondispersive, long-wave asymptotics, including their limits to single homogeneous fluid with a free surface, we showed that detachment of an interface from smooth contacts with a rigid boundary can occur via a singularity development in the form of shocks. In particular, the limit of vanishing upper fluid density reveals substantial differences between the evolution of top and bottom contact points, whereby at top contact points shocks may develop immediately in this limit. Moreover, our study identified several critical phenomena with zero velocity initial data, and symmetric interfacial profiles, as well as the non-negligible role played even at short times by dispersive wave effects in some cases. One notable element that emerges from this analysis is a kind of upper lid “suction,” or tendency, above a critical threshold determined by the density ratio between the two fluids, for the peak region of a slumping bulge of heavier fluid to be delayed in its downward motion with respect to nearby lower fluid. This causes the curvature at the peak to increase temporarily, before reversing sign in its downward acceleration.

The critical cases, besides their inherent interest, can also serve as a testbed for numerical simulations, including direct ones for the parent Euler equations with smooth stratification (albeit with a sharp enough transition between two nearly uniform densities to approximate a two-fluid system). Performing these computations, with extensive testing of their numerical accuracy (see also [2, 10]), demonstrated a remarkable robustness of the nondispersive model predictions, at least for the class of initial conditions and time scales of evolution we have focussed on, at both qualitative and, with few exceptions, quantitative levels. This is borne out of a comparison between interface profiles obtained through the numerical solutions of the two-layer model (4.3) and the direct numerical simulation of their Euler counterpart by overlaying the profile on the late time snapshots from figures 14, 15, 18, 19, as reported in figure 20. A notable example that emerges from this comparison is the prediction of “wing” structures (simple waves for the hyperbolic models) and their locations. These wings connect the core symmetric state centered at the origin with the constant background state: while only qualitatively similar, all the simulations with the full Euler system displayed the emergence of these features, whose existence does not seem to be immediately predictable from Euler dynamics alone. On the other hand, the evolution of connection points, figure 19, where the pycnoclines have discontinuous first derivatives show incipient shear instability which cannot be associated with hyperbolic-elliptic transition in the two-fluid model.

Among the results above summarized, our work has identified possible scenarios in the evolution of initial data where internal and rigid confining boundary come into contact. In particular, filling a dry hole in a body of water may occur with a singularity of global kind, whereby an infinite derivative develops over a whole segment simultaneously, as opposed to the standard case of pointwise derivative divergence at first. Conversely, with a non-negligible upper fluid density, the same initial data develop a standard shock with an interface that closes the gap in a zipper-like fashion (though it remains to be seen whether a left over filament of upper fluid, while becoming increasingly thinner, nonetheless persists at all times, as possibly suggested by the numerical simulations). It is worth remarking, at this point, that so long as viscosity is neglected as in the parent Euler equations, the symmetry of the initial conditions we have considered can be viewed as enforcing the boundary condition of a vertical wall at $x=0$, since antisymmetry of the horizontal velocity is equivalent to an impermeable boundary condition at the origin. Thus, we see that, according to the model, the impinging mass of water smashing into the wall produces a jet whose apex remains lower than the background state if the initial condition has a dry point at the wall, which is recovered only at large times after collision.

Our study opens future avenues of investigation which we are currently pursuing. In a separate paper, we will report on closed form expressions for the simple-wave wings, and also study how the evolution of the hyperbolic models can be continued after shocks develop. Further, the effects we have studied can be expected to be altered by the orientation of the confining boundaries, for instance when the bottom is a sloping “beach.” The classical analysis in [11] provides a hint of what could happen in this case, for a single-layer fluid, and would have to be extended to treat more general configurations and stratifications. These are challenging questions which we plan to address in future work.

Acknowledgments

We thank R. Colombo and M. Garavello for discussions and useful comments on the theory of quasilinear PDE’s. R.C. and G.O. thank S.L. Gavrilyuk for bringing to their attention, after a first draft of this work had been submitted for publication, reference [32] during the Summer school “Dispersive hydrodynamics and oceanography: from experiments to theory” 27 August - 1 September 2017, Les Houches (France). Support by NSF grants DMS-0908423, DMS-1009750, DMS-1517879, RTG DMS-0943851, CMG ARC-1025523, ONR grants N00014-18-1-2490, DURIP N00014-12-1-0749, ERC grant H2020-MSCA-RISE-2017 PROJECT No. 778010 IPaDEGAN, and the auspices of the GNFM Section of INdAM are all gratefully acknowledged. R.C. and M.P. thank the Dipartimento di Matematica e Applicazioni of Università Milano-Bicocca for its hospitality. In addition, the Istituto Nazionale di Alta Matematica (INdAM) and the program “Singularities and Waves in Incompressible Fluids” in the Spring of 2017 at the Institute for Computational and Experimental Research in Mathematics (ICERM), supported by the National Science Foundation under Grant DMS-1439786, are gratefully acknowledged for hosting the visits of R.C. (INdAM, Summer 16) and R.C. & C.T. (ICERM, Spring 17) while some of this work was carried out. Last, but not least, we would like to thank the anonymous referees whose attentive reading of the manuscript greatly helped improve it.

Appendix A: The limiting case of air-water system

We have seen in section 4 that the Hamiltonian density (4.6)

[TABLE]

gives rise for $r=0$ to the Boussinesq two-layer model smoothly in $r$ (i.e., the Hamiltonian is continuous in $r$ at $r=0$). The opposite limit $r\to 1$ (that is, the air-water system) is more subtle. In fact, as long as the interface does not get into contact with the top boundary $\eta=1$ of the channel, expression (4.6) coincides with the classical Airy Hamiltonian density

[TABLE]

However, the limit $r\to 1$ preserves the memory of the upper lid when the interface is close to the top boundary. Indeed, expanding (4.6) with respect to small $\rho_{r}=1-r$ yields

[TABLE]

For every fixed value $\rho_{r}\ll 1$, if the interface gets sufficiently close to the upper lid, the correction term of the Hamiltonian density can no longer be assumed to be smaller than the leading order term. In other words, an asymptotic expansion in the small parameter $\rho_{r}=1-r$ fails in a neighborhood of $\eta=1$.

Let us now consider the hyperbolicity region of the Hamiltonian system generated by (4.6), i.e., the domain in the $(\sigma,\eta)$-plane where the characteristic velocities $\lambda_{\pm}$ are real and distinct (see figure 21). It is easy to show that its boundary (the so-called transition line) is given by $H_{\sigma\sigma}=0$ and $H_{\eta\eta}=0$. The first equation gives the physical boundaries of the system: $\eta=0$ and $\eta=1$. The second equation is

[TABLE]

It follows that

[TABLE]

When $r\to 1$, the hyperbolic domain fills the whole configuration space $\mathbb{R}\times[0,1]$. However, the transition line $H_{\eta\eta}=0$ does not merge smoothly with the horizontal line $\eta=1$, since the last derivative in (A.5) tends to $\mp{2}/{3}$. As a consequence, if the distance of some point of the interface is smaller than the infinitesimal $1-r$, the system will experience a hyperbolic-elliptic transition if the weighted vorticity $\sigma$ is larger than an infinitesimal of the same order.

Appendix B: An example of simple-wave behavior for a smooth interface

In this appendix we provide an explicit example of evolution from an analytic initial interface in the simple-wave setting for the Airy system.

We consider the following family of interfaces at $t=0$,

[TABLE]

where $K$ is a nonnegative parameter, $n$ is a positive integer and $H(x)$ is the Heaviside function. For a simple wave, by choosing the minus sign in (3.10), the velocity field is

[TABLE]

The Riemann constant value $c=-2\sqrt{1+K}$ is fixed by the asymptotic zero-value of the velocity at $x\to-\infty$. The simple wave equation is

[TABLE]

so that its hodograph solution is

[TABLE]

In the region $x>(3\sqrt{K}-2\sqrt{1+K})t$ the solution remains constant,

[TABLE]

while if $x<(3\sqrt{K}-2\sqrt{1+K})t$ the solution is implicitly given by

[TABLE]

The velocity at the interface minimum is $u_{0}(0)=2\sqrt{K}-2\sqrt{1+K}<0$ independently of $n$, therefore a shock will always form there. The shock curve (3.14) is given by

[TABLE]

As recalled in section 3.1 the shock value $\eta_{c}$ is a solution of ${{\mathrm{d}}t}/{{\mathrm{d}}\eta}=0$, so that

[TABLE]

The solution of this equation obviously depends on the values of $K$ and $n$. For example, when $n=1$ we have

[TABLE]

and

[TABLE]

We remark that even if

[TABLE]

both $x_{c}$ and $t_{c}$ diverge when $K\to\infty$.

If $K=0$, the shock value is

[TABLE]

and therefore it is at the bottom of the channel only if $n=1$, see below. The other critical values are

[TABLE]

Consider now the case when $K=0$ and $n=1$: the initial velocity

[TABLE]

is not differentiable at $x=0$, and the critical values of the field $\eta$ is $\eta_{c}=0$. Hence the shock forms at the boundary of the domain, and one can check that $x_{c}=-2/3$ and $t_{c}=1/3$. We can conclude that for this example the smoothing does not introduce qualitatively new phenomena on the position of the shock formation with respect to the non-smooth cases considered in section 3.1 for simple-waves solutions of the Airy equations.