Singularity formation as a wetting mechanism in a dispersionless water wave model

TL;DR

This paper analyzes the formation of singularities in a dispersionless water wave model, using analytical and numerical methods to understand the evolution of solutions with initial discontinuities and dry spots.

Contribution

It provides a closed-form description of solution evolution up to collapse and develops spectrally accurate numerical methods to study post-collapse dynamics.

Findings

Solutions develop finite-time velocity and gradient catastrophes.

Dry spots are terminated by collapse of the parabolic core.

Numerical methods accurately track evolution beyond collapse.

Abstract

The behavior of a class of solutions of the shallow water Airy system originating from initial data with discontinuous derivatives is considered. Initial data are obtained by splicing together self-similar parabolae with a constant background state. These solutions are shown to develop velocity and surface gradient catastrophes in finite time and the inherent persistence of dry spots is shown to be terminated by the collapse of the parabolic core. All details of the evolution can be obtained in closed form until the collapse time, thanks to formation of simple waves that sandwich the evolving self-similar core. The continuation of solutions asymptotically for short times beyond the collapse is then investigated analytically, in its weak form, with an approach using stretched coordinates inspired by singular perturbation theory. This approach allows to follow the evolution after collapse by implementing a spectrally accurate numerical code, which is developed alongside a classical shock-capturing scheme for accuracy comparison. The codes are validated on special classes of initial data, in increasing order of complexity, to illustrate the evolution of the dry spot initial conditions on longer time scales past collapse.

Figures

Equations

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

$$\eta =\gamma (t)x^2+\mu (t),u=\nu (t)x.$$

$$\dot{\nu }+{\nu }^2+2\gamma =0,\dot{\gamma }+3\nu \gamma =0,\dot{\mu }+\nu \mu =0.$$

$$\nu (0)={\nu }_0,\gamma (0)={\gamma }_0,\mu (0)={\mu }_0,$$

$$\nu (\sigma )=\pm \sqrt{4{\gamma }_0{\sigma }^3+({\nu }_0^2-4{\gamma }_0){\sigma }^2},\gamma (\sigma )={\gamma }_0{\sigma }^3,\mu ={\mu }_0\sigma ,$$

$${\dot{\sigma }}^2=({\nu }_0^2-4{\gamma }_0){\sigma }^4+4{\gamma }_0{\sigma }^5.$$

$$\gamma ={\gamma }_0\exp \left(-3{\int }_0^t\nu (t^{\prime })\mathrm{d}{t}^{\prime }\right),$$

$$\dot{\nu }+{\nu }^2+2\gamma =0,\dot{\gamma }+3\nu \gamma =0.$$

$${\dot{\sigma }}^2=4{\gamma }_0({\sigma }^5-{\sigma }^4),$$

$$\dot{\nu }(t)<0⟹\nu (t)<0⟹\dot{\gamma }(t)>0.$$

$$\dot{\sigma }=2\sqrt{{\gamma }_0}{\sigma }^2\sqrt{\sigma -1}.$$

$$t(\sigma)=\frac{\displaystyle\sqrt{\sigma-1}+\sigma\,\arctan\big{(}\sqrt{\sigma-1}\hskip 1.0pt\big{)}}{2\displaystyle\sqrt{\gamma_{0}}\,\sigma}\,,\qquad\gamma=\gamma_{0}\sigma^{3}\,,\quad\nu=-{\dot{\sigma}\over\sigma}=-2\sqrt{\gamma_{0}}\,\sigma\sqrt{\sigma-1}\,,$$

$$t_c=\underset{\sigma \to +\infty }{\lim }t(\sigma )=\frac{\pi }{4\sqrt{{\gamma }_0}},$$

$$R_{\pm }=u\pm 2\sqrt{\eta },{\dot{x}}_{\pm }\equiv {\lambda }_{\pm }=u\pm \sqrt{\eta }$$

$$\big{(}\nu(t)\pm 2\sqrt{\gamma(t)}\,\big{)}x=\pm 2\sqrt{\gamma_{0}}\,\,x_{0}\,.$$

$$x_{\pm }(\sigma ;x_0)=x_0\frac{\sqrt{\sigma }\pm \sqrt{\sigma -1}}{\sigma }.$$

$$\eta(x,0)=\left\{\begin{array}[]{lc}Q,&x<-a_{0}\\ \gamma_{0}x^{2},&-a_{0}\leq x\leq a_{0}\\ Q,&x>a_{0}\end{array}\right.\,,\qquad u(x,0)=0.$$

$$a(t(\sigma ))=\frac{\sqrt{Q\sigma }-\sqrt{Q(\sigma -1)}}{\sigma \sqrt{{\gamma }_0}}.$$

$$t_s=\frac{\pi }{4\sqrt{{\gamma }_0}},$$

$$\begin{split}&\eta\big{(}a(t(\sigma)\big{)},t(\sigma))=\gamma(t(\sigma))\,\big{[}a(t(\sigma))\big{]}^{2}=\sqrt{Q}\,\sigma\left(\sqrt{\sigma}-\sqrt{\sigma-1}\,\right)^{2}\,,\\ &u\big{(}a(t(\sigma)),t(\sigma)\big{)}=\nu(t(\sigma))\,a(t(\sigma))=-2\sqrt{Q}\left(\sqrt{\sigma^{2}-\sigma}-(\sigma-1)\right)\,.\end{split}$$

$$\eta (a(t_c),t_c)=\frac{Q}{4},u(a(t_c),t_c)=-\sqrt{Q}.$$

$$x_{\pm }(\sigma ;x_0)=\frac{x_0}{\sqrt{\sigma }}\pm \frac{\sqrt{(\sigma -1)\left(x_0^2+{\mu }_0/{\gamma }_0\right)}}{\sigma },$$

$$t_{s}=\frac{\sqrt{{Q}/{\mu_{0}}-1}+({Q}/\mu_{0})\arctan\big{(}\sqrt{{Q}/{\mu_{0}}-1}\hskip 1.0pt\big{)}}{2\sqrt{\gamma_{0}}\,\,{Q}/{\mu_{0}}}\,,$$

$$\begin{array}{rl}\eta (x,t)=& \gamma x^2{\chi }_{(-a,a)}+N(x,t){\chi }_{(a,b)}+N(-x,t){\chi }_{(-b,-a)}+Q{\chi }_{(-b,b{)}^c}\\ u(x,t)=& \nu x{\chi }_{(-a,a)}+V(x,t){\chi }_{(a,b)}-V(-x,t){\chi }_{(-b,-a)}.\end{array}$$

$$V=2\sqrt{N}-2\sqrt{Q}.$$

$$N_{t}+\lambda_{+}\big{(}R_{+},R_{-}\big{)}\,N_{x}=N_{t}+\frac{1}{4}\big{(}3R_{+}+R_{-}\big{)}\,N_{x}=N_{t}+\big{(}3\displaystyle\sqrt{N}-2\displaystyle\sqrt{Q}\big{)}N_{x}=0\,,$$

$$3\sqrt{N}-2\sqrt{Q}=3\sqrt{\gamma (t_0)}a(t_0)-2\sqrt{Q}.$$

$$x=\left(3\sqrt{\gamma (t_0)}a(t_0)-2\sqrt{Q}\right)(t-t_0)+a(t_0).$$

$$\begin{array}[]{c}N(x,t)=\gamma(\sigma_{0})a(\sigma_{0})^{2}=\sigma_{0}Q\left({\sqrt{\sigma_{0}}-\sqrt{\sigma_{0}-1}}\right)^{2}\,,\vspace*{0.3cm}\\ 3\displaystyle\sqrt{N(x,t)}-2\displaystyle\sqrt{Q}=\Lambda(\sigma_{0})\equiv 3\,\displaystyle\sqrt{Q}\,\sigma_{0}\left(1-\sqrt{1-{1\over\sigma_{0}}}-{2\over 3\sigma_{0}}\right)\,,\end{array}$$

$$x=\Lambda ({\sigma }_0)\left(t-\frac{\sqrt{{\sigma }_0-1}+{\sigma }_0\arctan \left(\sqrt{{\sigma }_0-1}\right)}{2{\sigma }_0\sqrt{{\gamma }_0}}\right)+\frac{\sqrt{Q{\sigma }_0}-\sqrt{Q({\sigma }_0-1)}}{{\sigma }_0\sqrt{{\gamma }_0}}.$$

Full text

1 Introduction

The interaction of fluid layers with boundaries are arguably among the most challenging problems in the mathematical formulation of fluid dynamics. The Euler system of perfect fluids, while highly idealized, offers a physically relevant simplification of this dynamics in many setups, however it often presents several obstacles from an analytical standpoint. While a direct investigation of Euler equations could always in principle be carried out numerically, much insight can be gained from the study of reduced models that can extract the essential elements of confining-boundary/interface interactions, by incorporating 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 layer long-wave problem which is the focus of this work. Here, the usual free-surface boundary conditions of the Euler equations are enforced ab initio through the process of averaging. The simplest set-up to investigate the interaction of a free surface with a bottom boundary is that of a two-dimensional channel. As can be easily shown (see, e.g., [4]), contact points between the free surface and the bottom boundary can never be wiped out as long as the solution of either the Euler or the Airy systems remains smooth.

Hyperbolic quasilinear systems can develop shocks, whose study is a classic topic of the literature (see, e.g., [30] in fluid settings, or [31] for a more general esposition). The general subject of hyperbolic-parabolic transitions, to which the problem of wetting of dry points/intervals (“vacuum states”) belongs in the present context, has also been the subject of many investigations, both in the physical and in the mathematical literature (see, e.g., [1, 8, 23, 21]). In particular, the authors of [23] remark how shocks cannot propagate into vacuum states, and how a theory of how vacuum states can be filled by rarefaction waves may still lack a general, rigorous foundation from the mathematical perspective. To this end, various methods of regularization in the study of such vacuum states have been proposed, including entropy methods [23] and dissipative regularizations [22]. In this work we confine ourselves to the purely quasi-linear setting. By using specific classes of initial data whose evolution is amenable to analytical methods, we demonstrate that the formation of a singularity can circumvent the seemingly unphysical conclusion of the persistence of the dry point. The class of initial data we consider are exact global solutions of the Airy system up to the time of singularity, and in certain sectors provide local solutions even after this time. These solutions evolve from initial conditions obtained by splicing together constant and time-dependent solutions in the form of local parabolic profiles for the free surface elevation, which can be viewed as self-similar in the proper time-like variables. The singularity for these data manifests itself as a “global” shock, that is, a gradient catastrophe occurs over a whole interval in the range of the dependent variables, as opposed to the generic case of an infinite derivative at single point in this range.

Our main goal is to study the shock problem emanating from the singularity for a suitable class of initial conditions. In this, our approach differs substantially from previous works on the subject (see, e.g, [16, 24, 13, 18, 11, 29, 2, 20]), and, in particular, from the one more recently set forth by Dubrovin and collaborators (see, e.g., [9, 10, 12]) where the main focus is the conjectured universality behaviour of hyperbolic (and elliptic) shocks, universality being intended both with respect to general dispersive Hamiltonian regularizations and to initial conditions. In this context, the Airy system, viewed as the nondispersive version of the defocussing Nonlinear Schrödinger equation (NLS), and the non-generic nature of gradient catastrophe in the presence of a vacuum point, was investigated using specific initial data in [25]; the evolution through singularity was illustrated by numerical simulations the full dispersive NLS, which resolved possible shocks into dispersive wave trains. Of course, in our case the evolution of weak (shock) solutions, and more generally the presence of jump discontinuities in the derivatives of the dependent variables, conflicts with the derivation of Airy’s model from Euler equations by long wave asymptotics. A full discussion of this point would fall beyond aims of this study, and we leave this to the fairly extensive literature (see, e.g, [30, 31]). However, it is worth mentioning here that in our recent work [4], while testing the validity of non-dispersive models for one and two layer fluids we found that, at least for the class of initial data we study, hyperbolic shallow water models do provide a reasonably good prediction (both qualitatively and quantitatively) of the dynamics of the parent equations, on short time scales before other physical effects neglected by the model, such as dispersion, become important.

The details of the class of initial data and the parabolic solutions we consider are reported in section 2, where we focus on the case of an initial dry contact point with the bottom boundary. We then briefly contrast, in §2.1, the evolution of the “dry” case and its singularity formation with that of a nearby class of initial data which maintain a fluid film over the bottom. The shock formation with parabola collapse is examined more closely in §2.2, with the aim of bringing forth the analytical obstacles to weak solution continuations. Section 3 introduces new classes of initial data designed to illustrate the main mechanisms of after-shock continuation by eliminating technical obstacles of the original case. These are overcome in §3.2 by designing “unfolding” coordinates that magnify space-time region near the singularity at the collapse time, turning the initial value problem on the whole real line into a boundary value problem on a strip whose end points are pinned to the propagating shocks. The new formalism allows for asymptotic solutions to be computed explicitly. Finally, in section 4, we illustrate the results with numerical algorithms for both unfolding and physical coordinate systems, and validate the numerical simulations with the analytical results, with an eye at isolating possible accuracy shortcomings of the algorithms.

2 Parabolic solutions for Airy model

In a previous paper [4] we introduced and studied, in various instances, piecewise smooth initial data for $1+1$-dimensional models approximating the Euler equations, where the density-related variable $\eta$ is obtained by splicing together constant heights with (self-similar) parabolic profiles. In the specific case of the classical Airy model for shallow water dynamics, in suitable non-dimensional variables,

[TABLE]

local solutions can be constructed, for $(x,t)$ belonging to certain sectors in space-time, using the simple spatial functions

[TABLE]

The evolution described by these functions can be viewed as (generalized) self-similar solutions of (2.1) since (i) they can be cast in the form $u(x,t)=k_{u}(t)f_{1}({x}/{l(t)})$, $\eta(x,t)=k_{\eta}(t)f_{2}({x}/{l(t)})$, for appropriate functions $k_{u}$, $k_{\eta}$, $f_{1}$, $f_{2}$ and $l$ of one variable, and (ii) inserting this prescription in the Airy equations leads to the closed exact system of ordinary differential equations (ODEs) for the coefficients $(\nu,\gamma,\mu)$

[TABLE]

With initial conditions

[TABLE]

the solutions of system (2.3) can be expressed implicitly by

[TABLE]

where the auxiliary variable $\sigma$ solves the first order nonlinear ODE

[TABLE]

The parabola curvature sign is invariant in time because the quadrature of system (2.3) implies

[TABLE]

so that $\gamma$ remains positive if it is initially so.

The study of self-similar solutions (2.2), which were independently noted in [26], has first appeared in [4]. Here we focus on the case of solutions with positive curvature $\gamma(t)>0$ starting initially with a dry point, a point where the surface touches the bottom ($\mu_{0}=0$). For the shallow water model, $\eta=0$ is a line that plays a double role: it is a degenerate characteristic in the hodograph plane, and it is the “sonic line,” i.e. the locus of points where the system loses hyperbolicity because the characteristic velocities coincide. Such peculiar feature has some interesting consequences. In particular, from (2.5) we have $\mu(t)=0$, and the minimum of the parabola (where $\eta=0$) remains dry for as long as the solution maintains its local regularity.

When initial velocities are sufficiently large, $\nu_{0}\geq 2\sqrt{\gamma_{0}}$, it can be shown that the curvature $\gamma(t)$ tends monotonically to zero and therefore no loss of regularity in the form of shocks can occur. For small velocities the local parabola can collapse to a vertical segment, a “global” shock, in finite times. This shock effectively reduces a portion of the free surface $\eta$ support containing the dry point to a set of measure zero, which can then be removed from the subsequent dynamics, thus providing a mechanism for the disconnection of the free surface from the bottom.

Next, we derive closed form solutions by quadratures of the first pair of equations in system (2.3),

[TABLE]

(Since the third equation of (2.3) is slaved to these two, its solution follows by quadratures; for the “dry spot” case $\mu_{0}(t)=0$ for all times, and no additional quadrature is required.) For sake of simplicity we restrict our study to zero initial velocities, i.e., $\nu_{0}=0$. In this case equation (2.6) becomes

[TABLE]

and, since $\gamma(t)>0$ and $\nu_{0}=0$, for all $t>0$,

[TABLE]

Therefore ${\gamma}(t)$ and $\sigma(t)$ are strictly increasing functions of time, ${\gamma}(t)\geq\gamma_{0}$ and $\sigma(t)\geq 1$, and equation (2.9) can be rewritten as

[TABLE]

Its implicit solution

[TABLE]

shows that $\sigma$ is a time reparameterization. As the variable $\sigma\to+\infty$ the curvature of the parabola $\gamma$ blows up at the finite time

[TABLE]

while the coefficient $\nu$ diverges, $\nu(t)\to-\infty$ as $t\to t_{c}^{-}$, see (2.5).

Next, we look at characteristics corresponding to parabolic initial data. As well known, see e.g. [31], the Riemann invariants and characteristic equations for system (2.1) are

[TABLE]

which yield, for the characteristic curves emanating from the initial position $x(0)=x_{0}$ in the support region of the parabolic initial data,

[TABLE]

Hence, using solutions (2.5) for $\nu(t)$ and $\gamma(t)$ in terms of the auxiliary time-like variable $\sigma$, the characteristic solutions can be expressed as

[TABLE]

The geometry of these characteristic curves is depicted in figure 1.

The initial configurations considered in [4] can be called truncated parabolae, inasmuch as the profile of the height $\eta$ is obtained by splicing a parabola with a constant height $\eta=Q$ for $|x|$ larger than a certain distance, $a_{0}>0$ say, from the origin, or $x_{0}\in[-a_{0},a_{0}]$. We can represent these initial conditions as

[TABLE]

Continuity requires $Q>0$ and $a_{0}=\sqrt{Q/\gamma_{0}}$. Time evolution initially preserves continuity and piecewise differentiability, although the original parabolic section splits into three sections: a shrinking (positively curved) parabola centred at $x=0$, and two convex curves symmetrically placed with respect to the origin, joining the parabola with the constant level $\eta=Q$. The behavior of the solution for $\mu_{0}=0$ is shown in figure 2. In summary, the extrema $\pm a(t)$ of the parabolic section evolve in time from their initial values $\pm a_{0}$, and as functions of the auxiliary variable $\sigma$ they can be found at

[TABLE]

By comparison with the characteristics (2.16), as expected by the general theory of singularities (see e.g. [31]), the point $a(t)$, being a point of continuity with discontinuous first derivative of the dependent variables, moves along a characteristic. When the curve $x=a(t)$ crosses $x=0$ it meets its symmetric counterpart $x=-a(t)$. Let $t=t_{s}$ denote such crossing time. The value $\sigma_{s}$ of $\sigma$ at time $t_{s}$ can be obtained from (2.18) by setting $a(\sigma_{s})=0$, so that $\sigma_{s}=+\infty$ and

[TABLE]

which coincides with the parabola collapse time $t_{c}$. This is special to the dry-point case $\mu_{0}=0$; in general the two times $t_{s}$ and $t_{c}$ are different, with $t_{s}<t_{c}$ if $\mu_{0}>0$, as we shall see below. Unless otherwise noted, we shall focus from now on onto the dry case, for which $t_{s}=t_{c}$, and use the latter notation to denote the time at which the singularity occurs. The height of the interface and the velocity at $x=a(t)$ are given by

[TABLE]

At time $t=t_{c}$ such quantities (for the solution with $x>0$ in the limit $t\to t_{c}^{-}$, say) are

[TABLE]

The support of the parabolic section disappears at $t=t_{s}$, and the solution develops a jump of amplitude $2\displaystyle\sqrt{Q}$ in the velocity. The segment $\eta\in[0,Q/4]$ to which the parabola has reduced can also be viewed as double (from left and right of $x=0$) shock for the surface elevation $\eta$. The removal of this segment, i.e., the value $\eta(0,t_{c})=0$ and the replacement with $\eta$’s limiting value $\eta(0,t_{c})=\lim_{|x|\to 0}\eta(x,t_{c})=Q/4$ makes the surface elevation continuous. Hence the gradient catastrophe may be viewed as the mechanism responsible for the detachment of the interface from the bottom (see figure 2). Notice that the gradient catastrophe in this exact solution differs from the more common case whereby the derivative diverges initially at a single point. This can be viewed as the consequence of characteristics of the same family originating from a finite length segment, such as $x_{0}\in(0,a_{0})$, all crossing at the same point in the $(x,t)$-plane (see figure 4). For this reason we will often refer to $t_{c}$ as the “global shock” time.

2.1 The “wet” case $\mu_{0}>0$

When the initial parabola’s minimum is above the bottom plate, missing it by a non-zero amount $\mu_{0}>0$, however small, the ($\eta$,$u$) evolution is substantially different from the contact case above, though the actual analytical calculations are similar to those in the previous section. We omit all the details (for a full account, see [4]) and report here the final results, focussing on the differences between the “near miss” of the surface contact and its contact counterpart $\mu_{0}=0$. The characteristics for $\mu_{0}>0$ are

[TABLE]

and those for the extrema of the parabola $x_{0}=\pm a_{0}$ cross at time $t_{s}<t_{c}$,

[TABLE]

at which the elevation $\eta(x,t)$ recovers its background value, $\eta(0,t_{s})=Q$.

The difference between the two evolutions is illustrated by figure 3. Unlike the case $\mu_{0}=0$, the surface $\eta(x,t)$ returns to the background elevation $Q$ in a neighborhood of the origin, and both the surface elevation $\eta$ and the fluid velocity $u$ maintain continuity beyond the time $t_{c}$. This is in contrast with the dry-point case at time $t=t_{c}$, when the boundary-touching parabola collapses to a segment, forming a global shock $0\leq\eta(0,t_{c})\leq Q/4$, with corresponding fluid velocity shock $-\sqrt{Q}\leq u(0,t_{c})\leq\sqrt{Q}$.

2.2 The connecting “shoulder” states

In this section we discuss the evolution of the “shoulders” that join the parabolic sectors with the constant far-field states $\eta=Q$, by using an approach akin to that of the classical piston problem (see e.g., [30]). The behaviour of the solution corresponding to the initial data (2.17) before the time $t_{c}$ (global shock or collapse time) is obtained by joining three regions of different analytical behavior. Thus, along with a far-field region, $|x|>b(t)$ say, where both $\eta$ and $u$ are constant, and with a region centered around the origin $|x|<a(t)$, where the solution maintains the parabolic form described in the previous section, we have connecting regions $a(t)<|x|<b(t)$ where the solution is determined by simple waves, which we denote by the pair $(N(x,t),V(x,t))$. With the standard notation for characteristic functions $\chi_{I}$ of the support interval $I\subseteq\hbox{{I\hskip-1.50696ptR}}$, the whole solution can be written explicitly as

[TABLE]

Here the functions $a(t),b(t)$ describe the evolution of the boundaries of the above-mentioned regions where the solution of system (2.1) is continuous but might not be differentiable. With reference to figure 4, these parts of the solution live in the $S$ region, where the pair $(N,V)$ constitutes a simple wave (thanks to the symmetry $x\mapsto-x$, we can restrict to $x\geq 0$).

The $P$ region in figure 4 depicts characteristics (2.22) and corresponds to the parabola part of the solution. In the $C$ regions both Riemann invariants are constant and both families of characteristics are straight lines. In the right $S$ region, the Riemann invariant $R_{-}=u-2\sqrt{\eta}\equiv V-2\displaystyle\sqrt{N}$ is constant (independent of $x_{0}$), since its value is $R_{-}(x_{0},0)=-2\displaystyle\sqrt{Q}$ at all points $(x_{0},0)$ with $x_{0}\geq a(0)=\sqrt{Q/\gamma_{0}}$, and this value is conserved along characteristics $x_{-}(t;x_{0})$ (dashed curves in figure 4). Hence the pair $(N,V)$ is functionally related, i.e., it is a simple wave with

[TABLE]

Thus, in this region $S$, the evolution is governed by

[TABLE]

so that the non-dashed characteristic curves in the right $S$ region are straight lines starting from $(a(t_{0}),t_{0})$ with (inverse of the) slope given by

[TABLE]

Note that the special value $N=4\,Q$/9 makes the characteristic line vertical, which in turn makes this value constant for all times at the corresponding $x$, x\equiv x_{d}=\big{(}\displaystyle\sqrt{3Q/\gamma_{0}}\big{)}/4, or $N(x_{d},t)=4\,Q/9$ for $t>t_{0}$. Since the value of $N$ is invariant along on straight lines defined by (2.27), we have that $N(x,t)=\gamma(t_{0})a(t_{0})^{2}$, where $t_{0}$ is given implicitly by

[TABLE]

It is convenient to use the parameter $\sigma_{0}$ corresponding to $t_{0}$ defined by relation (2.12). Using (2.18) and the second equation in (2.5) yields

[TABLE]

where $\sigma_{0}=\sigma_{0}(x,t)$ is implicitly given by

[TABLE]

From the shoulder solutions’ viewpoint, an alternative definition of the global shock time $t_{c}$ is that of the catastrophe condition, obtained by finding the first time at which the partial derivative $N_{x}$ becomes infinite. From (2.29) and (2.30) we see that this happens when

[TABLE]

This shock time for the shoulders coincides with the shock time for the parabola’s collapse to a segment (and with the crossing time $t=t_{s}$ at which $a(t)=0$),

[TABLE]

Similarly, the shock position of the shoulders coincides with the shock position by the collapsed parabola, i.e., $x(t_{c})=0$. Once again, the vertical segment $\eta\in[0,Q/4]$ at $x=0$, $t=t_{c}$ to which the parabola has collapsed, can be suppressed by setting $\eta(0,t_{c})=Q/4$ and $\eta(x,t_{c})=N(|x|,t_{c})$, thus restoring continuity at $x=0$ for the surface elevation $\eta$.

Past the global shock time $t=t_{c}$, the evolution of solutions of system (2.1) continues in weak form. As we shall see below, weak solutions of (2.1) at short time $t>t_{c}$ are determined by the local behaviour of the simple-wave shoulders around the origin, i.e., before the shocks that emanate from $x=0$ have travelled sufficiently far away from a neighbourhood of the origin. While of course simple wave solutions are globally known implicitly by the characteristics for (2.26), an explicit, albeit asymptotic, form is useful for deriving analytic expressions of $\eta$ and $u$ near $x=0$. This requires knowledge of a new initial condition, $N(x,t_{c})\equiv N_{0}(x)$ say, in the shifted time $\tilde{t}\equiv t-t_{c}$, for the Hopf-like equation (2.26). The form of $N_{0}$ is given by the joined left and right shoulders at $\tilde{t}=0$, continued through $x=0$ by taking for $N_{0}(0)$ its limiting value $Q/4$. Putting $t=t_{c}$ in (2.30) and seeking the asymptotic expansion $\sigma_{0}\to\infty$ of the right-hand side yields

[TABLE]

or, inverting,

[TABLE]

This asymptotic solution for the implicit relation (2.30) gives the leading order behaviour for $N_{0}(x)$,

[TABLE]

The corresponding asymptotics for the velocity component can of course be obtained from this expression and the simple wave relation (2.25). Thus, around the origin the spatial behavior of the left and right shoulders, connected at $x=0$, $\eta=Q/4$, is that of a branch point (cusp singularity) $O(x^{2/3})$. As we shall see below, this singular behaviour has non-trivial effects on the global evolution of interface $\eta$ and velocity $u$ after the shock time.

A few remarks are now in order. First, while system (2.8) and its solution (2.12) have been introduced for the case of parabolic-shape layer thickness at the dry-point, with zero velocity initial data, most of the corresponding results actually apply to any (symmetric) initial condition for $\eta$ that admits a convergent Taylor series expression in the neighbourhood of the contact point. As shown in the Appendix, in the presence of dry point $\eta=0$, the evolution of the leading order coefficients of the Taylor series for $\eta$ and $u$ decouples from all their higher order counterparts. That is, if by $\eta_{n+1}(t)$ and $u_{n}(t)$ we denote the $n$-th order pair of Taylor coefficients (cf. (A1)), the system of ODE’s governing the evolution of $\eta_{1}(t)$ and $u_{0}(t)$, equation (A2), coincides with system (2.8). Second, all the higher order coefficients in the Taylor series are determined recursively from the preceding ones through a well defined hierarchy of ODE’s. At each stage beyond the leading order pair $(\eta_{1},u_{0})$, these ODE’s are linear in the new unknowns $\eta_{n+1}(t)$ and $u_{n}(t)$, and hence each $n$-th pair inherits the singularity determined by that of the first pair equations, i.e., system (2.8), which as shown above admits a closed form solution. Hence, the catastrophe time at $x=0$ is determined by the local curvature blowup, or $t=t_{c}$, through relation (2.13) for any smooth initial data $\eta(x,0)$ and $u=0$ that have the assumed property of nonzero curvature at the contact point. This generalizes a result reported in [25] obtained through a different approach for the case (in our notation) $\eta(x,0)=\tanh^{2}(x)$, $u(x,0)=0$, for the Airy system viewed as the zero dispersion limit of the defocusing Nonlinear Schrödinger equation. Further, the singularity time expression (2.13) signals that as the curvature of $\eta$ at the dry-point vanishes, or $\gamma_{0}\to 0$, the singularity might actually never occur in finite time, and that when the contact point at $x=0$ is a zero of higher order for $\eta(x,0)$, a different result might be expected. In fact, it is easy to show from the Taylor series approach, and the ensuing hierarchy of time dependent coefficients, that the first nonzero Taylor coefficient of $\eta$ is time independent, that is, the local shape of $\eta$ in a neighbourhood of the origin is essentially invariant in time. Lastly, note that the limiting form of the simple wave (2.35) suggests a property that was missed in previous investigations with smooth initial data: the spatial dependence of $\eta$ and $u$ at the singularity is a non-analytic cusp scaling like $x^{2/3}$. Thus, the smoothness loss occurs with a derivative $\eta_{x}$ that diverges as $x=0$, unlike the case of the divergence of curvature to a corner whereby the derivative would jump between finite values to the left and right of $x=0$.

From the initial condition $N_{0}(x)$ at time $\tilde{t}=0$ the evolution of the shoulders generates multivaluedness, as the characteristics $x_{-}(t)$ with $x_{0}>a_{0}$ from right-region $C$ cross their counterpart with $-a_{0}<x_{0}<0$ from the left-region $S$ (and vice-versa, by symmetry, for $x_{+}(t)$ characteristics). As well known for hyperbolic systems, multivaluedness can be remedied by selecting a shock location, $|x|=x_{s}(t)$ say, and allowing the evolution of solutions governed by system (2.1) to continue in the weak sense. With respect to the shifted time $\tilde{t}=t-t_{c}$ (and immediately dropping the tilde as long as this does not generate confusion), the initial condition for the shock position is clearly $x_{s}(0)=0$; for later purposes it is informative to eliminate, asymptotically for $t\to 0^{+}$ (or $t\to t_{c}^{+}$ in the non-shifted time), the implicitness of the functional dependence $N(x,t)$ at the shock location. This can be achieved directly from (2.26) taking $N_{0}(x)$ as a new initial condition. In fact, along the characteristic originating from $x=x_{0}>0$, which for clarity we simply denote by $x=X(t)$ since we are treating this as a new initial value problem, we have

[TABLE]

Hence the characteristics through the point $x=x_{s}(t)$ at time $t$ are defined implicitly by $X(t)=x_{s}(t)$, i.e., by an initial value $x_{0}$ such that

[TABLE]

For $t\to 0^{+}$ the asymptotic form of $N_{0}$ (2.35) can be used in this relation, which results in a cubic equation for $x_{0}^{1/3}$

[TABLE]

where we have defined the shorthand

[TABLE]

Since both $t$ and $x_{s}(t)$ are small (and expected to be of the same order) as $t\to 0^{+}$, the dominant balance in equation (2.36) is

[TABLE]

so that the explicit asymptotic form of $N(x,t)$ at the shock location for $t\to 0^{+}$ is

[TABLE]

3 Beyond the collapse time: analytical approach

The analysis developed in Section 2 describes the evolution of our class of initial conditions until the catastrophe time $t_{c}$. In particular, we have seen how the singularity of the velocity field $u(x,t)$ occurring at $t=t_{c}$ for dry-spot initial data leads to a reconnection of the left and right water masses into a single mass with a cuspy “well” of depth $Q/4$ at $x=0$ (see figure 2a, last time/panel). Our next task is to discuss the behavior of the system after this collapse and ensuing global (vertical segment) shock. We shall tackle this problem from an analytical perspective first, with a numerical study reported in the following section §4.

3.1 The three main examples

We begin by briefly recalling some fundamental properties of shock fitting in the context of our model equation (2.1). We shall focus here on the salient features of the original problem (2.17), after the parabola collapse and the zero-measure stem connecting the bottom to the elevation $Q/4$ is removed. As discussed above, the general structure giving rise to shocks can hence be taken to be that of simple waves “shoulders” which have merged at the collapse time. Thus, the solution (2.29) and (2.30) for the Airy model (2.1) is comprised of two simple waves corresponding to the Riemann invariants (2.14), $R_{\pm}=\hbox{const.}$ respectively for $x>0$ and $x<0$. For both signs the related characteristic velocities $u\pm\sqrt{\eta}$ increases when $u$ increases, and the Lax conditions [21] are satisfied and lead to shock jump relations enforcing conservation of mass and momentum, with densities $\eta$ $\eta u$ respectively. The shock speed $\dot{x}_{s}(t)$ for the shock position $x=x_{s}(t)$ is then given (see, e.g., [31]) by

[TABLE]

and the consistency between the two expressions can be manipulated to

[TABLE]

where, as customary, $[f]=f_{+}-f_{-}$ denotes the jump across the shock from right ($+$ subscript), to left ($-$ subscript). Next, we consider three realizations of the structure described above in order of increasing complexity, with the parabola case (2.17) at collapse being the most involved due to the loss of smoothness at $t=t_{c}$ as evidenced by (2.35).

3.1.1 The “double-Riemann” (toy) model

The simplest “cartoon” of evolution after the collapse of the dry-spot parabola of §2 is arguably that of the Riemann problem obtained by taking an initial condition determined solely by the asymptotic value $\eta(0,t_{c})=Q/4$ and a jump at $x=0$ of the velocity $u(0,t_{c})$ of amplitude $2\displaystyle\sqrt{Q}$,

[TABLE]

Such an initial condition for the Airy system (2.1) can be viewed as mimicking the local behavior around $x=0$ at the global shock time $t=t_{c}$ created by the merging of the shoulder simple waves. This Riemann problem can itself be considered as that emerging at the catastrophe time from the evolution of a piecewise continuous initial condition obtained by reverse evolution of the classical expansion wave emanating from a step-like initial condition, as well known (see, e.g., [17]). Thus, prior to the catastrophe, we consider a solution of the form, which we dub “double-Riemann,”

[TABLE]

with

[TABLE]

When $t=t_{c}$ a global shock forms, with the elevation becoming constant, $\eta=Q/4$, except at $x=0$ where $\eta=0$, and with the velocity jumping at $x=0$, $u=-\displaystyle\sqrt{Q}\,{\rm sgn}(x)$. The evolution is not affected by the value of the new initial data at one point, so we can assume that $\eta={Q/4}$ everywhere. Thus, in our case for $x>0$, $f_{+}$ corresponds to the constant “background” solution $\eta=Q/4$ and $u=-\displaystyle\sqrt{Q}$, while $f_{-}$ is determined by the solution to the cubic equation for $\eta_{-}$ which follows from (3.2) when $u_{-}=0$. We obtain

[TABLE]

where

[TABLE]

Here $Q^{*}$ is given by the unique solution of the cubic equation for $\eta_{-}$ which follows from the compatibility condition (3.2)

[TABLE]

subject to $\eta_{-}>Q/4$. The decimal value of this root is approximately

[TABLE]

which leads to the (constant) shock speed

[TABLE]

The evolution is depicted in figure 5.

3.1.2 The “double-Stoker” model

In order to capture the time evolution and spatial dependence of the initial conditions that emerge from the parabola collapse, the toy double-Riemann problem model described above is not adequate and can be modified as follows. We consider the evolution from an initial condition obtained by “cutting” and “splicing” together two Stoker (parabolic) rarefaction wavescrossing at $x=0$, $\eta(0,0)=Q/4$ (for relevant definitions, see e.g., [30], §10.8),

[TABLE]

where

[TABLE]

This choice of initial data mimics the state reached at the collapse time for the “full” case $N(x,t_{c}),V(x,t_{c})$, when the global shock at the origin can be removed from the simple wave shoulders in contact at $x=0$, $\eta(0^{\pm},t_{c})=Q/4$. Just as in that case, the evolution from initial data (3.11) maintains the elevation $\eta(x_{d},t)=4Q/9$ and $u(x_{d},t)=\mp 2\displaystyle\sqrt{Q}/3$ at all times (for as long as shocks do not reach the “hinge” position $x=\pm x_{d}$). The location $x_{Q}$ corresponds to the moving point where the rarefaction waves reaches the background elevation $Q$ and the fluid velocity vanishes. The analogous point for the full case moves at the same characteristic speed $\sqrt{Q}$, but at the initial time $t=t_{c}$ putting $g_{0}=\gamma_{0}$ yields a different initial position for this point. We dub this initial condition the “double-Stoker” case.

The evolution out of these initial data forms, at time $t=0^{+}$, two shocks which move away from the origin, just as in the double-Riemann case of section 3.1.1. However, unlike that case, a time-dependent shock speed and a non-trivial spatial and temporal dependence for $\eta$ and $u$ in the region between shocks can now be expected. For $t>t_{c}$, $|x|>x_{s}(t)$, the functions $N_{S}$, $V_{S}$ continue to be solutions of Airy’s system outside the region $|x|<x_{s}(t)$ between shocks.

3.1.3 The “full” case

This is the case of initial data (in the shifted time $t-t_{c}$)

[TABLE]

analyzed at the end of §2.2. The surface elevation $\eta$ is continuous with a cusp point $O(x^{2/3})$ at the origin, while the corresponding velocity $u(x,0)$ has a jump at the origin of amplitude $2\displaystyle\sqrt{Q}$. The local behavior at early times $t\to 0^{+}$ in a neighbourhood of the origin is summarized by the asymptotic behaviour (2.39). The schematic of the initial condition and subsequent short time evolution is qualitatively indistinguishable from that sketched in figure 6 for the double-Stoker case (though for the velocity $u(x,t)$ there are substantial differences between the two cases, as we shall see below).

3.2 Shock-fitted coordinates: unfolding the collapse singularity

Motivated by the tools of singular perturbation theory (see, e.g., [19]), in order to zoom onto the region $|x|<x_{s}(t)$ it is convenient to define “unfolding” coordinates $(\xi,\tau)$ to set the shock positions at fixed locations in time, e.g., $\xi=\pm 1$, which can be achieved by the new independent variables

[TABLE]

so that

[TABLE]

With this mapping, the Airy’s system (2.1) assumes the form (with a little abuse of notation by maintaining the same symbols for dependent variables)

[TABLE]

Here $\dot{x}_{s}$ is a placeholder for the expression that couples the evolution equation for the physical time $t$ to the new evolution variable $\tau$. The shock position evolves according to the equation that defines $\dot{x}_{s}$ in terms of the jump amplitudes $[\eta]$ and $[\eta u]$:

[TABLE]

so that, in the new variables,

[TABLE]

Thus, the system governing the evolution of the “inner” solution between shocks consists of this ordinary differential equation together with the partial differential equations (3.16), to be solved within the strip $\xi\in[0,1]$ (by symmetry only half the $\xi$ domain $[-1,1]$ can be used) subject to the boundary conditions

[TABLE]

The first equality is a consequence of the antisymmetry of the velocity, $u(\xi,\tau)=-u(-\xi,\tau)$. The second relation expresses the compatibility condition (3.2) in terms of the new independent variables. Here and in what follows, we suppress subscript labels that differentiate between the different cases of simple wave solutions which bracket the region between shocks, so long as it does not generate confusion; further, recall that for simple waves $V$ and $N$ are functionally related so that for all our cases $V$ can be thought of as a placeholder for the quadratic relation $V=2\displaystyle\sqrt{N}-2\displaystyle\sqrt{Q}$.

With the shorthand notations

[TABLE]

and

[TABLE]

where $s_{0}$ is the initial shock speed defined by (3.10), the full system of evolution equations for the inner shock region $(\xi,\tau)\in[0,1]\times\hbox{{I\hskip-1.50696ptR}}$ can be written compactly as

[TABLE]

with boundary conditions

[TABLE]

As $\tau\to-\infty$ “initial data” for $u$ and $\eta$ can be assigned so that they are compatible with the boundary conditions (3.19). This can be done by taking $u(\xi,\tau)\to 0$ and $\eta(\xi,\tau)\to Q^{*}\simeq 3.49396\,N(0,0)$ as $\tau\to-\infty$, because of the definition (3.9) of $Q^{*}$ as the solution of the equation

[TABLE]

constrained by the condition $Q^{*}>N(0,0)$ (with $N(0,0)=Q/4$ for all cases considered here).

We remark that system (3.22) shares the same Riemann invariants as those of the original model, although the corresponding characteristic eigenvalues now depend explicitly on the unfolding coordinates $\xi$ and $\tau$ as well. The Riemann invariant form of system (3.22) is

[TABLE]

We do not attempt here to find solutions in closed form of the initial boundary value problem for system (3.22); this is most efficiently approached numerically. Analytical progress can be made however by focussing on the short time behavior after shock formation, i.e., as $\tau\to-\infty$. In this limit, the initial evolution is governed by the perturbation of system (3.22) around the initial data $\eta=Q^{*}$, $u=0$ and $t=0$. Denoting by $\tilde{\eta}$ the difference $\tilde{\eta}=\eta-Q^{*}$, and abusing notation a little by immediately dropping the tilde while maintaining the same symbols for the dependent variables, unless necessary to avoid confusion, the governing equations for short times after shock formation are

[TABLE]

with boundary conditions

[TABLE]

Here we have defined

[TABLE]

and the gradient components of $\phi$ and $\psi$ with respect to the fields $\eta$ and $u$ are, respectively,

[TABLE]

Note that, as $\tau\to-\infty$, from the definitions (3.9) and (3.10) of the elevation $Q^{*}$ and shock speed $s_{0}$ at $t=0^{+}$, respectively, the boundary “forcing” terms $\phi_{0}(\tau)$ and $\psi_{0}(\tau)$ both vanish, with asymptotic rate dictated by the analytic feature of simple-wave shoulder solutions $N(e^{\tau},t(\tau))$ in this limit. Thus, for instance, for the double-Stoker case $N$ (and $V=2\displaystyle\sqrt{N}-2\displaystyle\sqrt{Q}$) are analytic functions of $e^{\tau}$ and $t(\tau)$, so that both $\phi_{0}(\tau)$ and $\psi_{0}(\tau)$ are of order $O(e^{\tau})$ as $\tau\to-\infty$, while for the full case the branch-point singularity evidenced by (2.35) leads to order O\big{(}e^{{2\over 3}\tau}\big{)} for these quantities in the same limit. Conversely, in general all gradient components $\phi_{\eta}$, $\phi_{u}$ and $\psi_{\eta}$ in (3.2) have finite non-zero limits as $\tau\to-\infty$.

The general solution of the perturbation system (3.26) (which is essentially the classical wave-equation in a space- and time-varying medium) satisfying the boundary condition at $\xi=0$ is

[TABLE]

for any function $F(\cdot)$ of sufficient regularity, with the function $\Phi$ defined by

[TABLE]

The asymptotic conditions on the system’s solution $(\eta,u)$ as $\tau\to-\infty$ requires $F(0)=0$, and substitution of these expressions evaluated at $\xi=1$ into the boundary condition (3.27) leads to a functional equation for $F$, coupled to the evolution equation for $t(\tau)$ in system (3.26).

From this general formulation of the perturbed unfolding problem, it is useful to derive identities that can be used for validation purposes. The initial time evolution of the curvature of $\eta$ and the slope of $u$ in a neighbourhood of the origin can be explicitly computed directly from the perturbed evolution equation (3.26). In fact, at $\xi=0$, the first two equations of system (3.26) yield

[TABLE]

and

[TABLE]

Further progress towards closed form expressions depends on the particular asymptotic behaviour of the simple-wave shoulder solution $N(e^{\tau},t(\tau))$. We next examine each of our three representative cases above, the “double-Riemann,” “double-Stoker” and “full” cases, highlighting their respective differences in the unfolding coordinate formulation.

3.2.1 Shock unfolding for double-Riemann initial data

The weak solution and shock fitting after the collapse of the initial parabola for this case has already been achieved in §3.1.1 following classical results for the Riemann problem. In fact, it is easy to see that in the unfolding formulation above the inner problem between shocks is time independent and consistent with the initial conditions as $\tau\to-\infty$, since the fields $N$ and $V$ are constant, $N=Q/4$ and $V=-\displaystyle\sqrt{Q}$, and compatible with these initial data. Thus, the second boundary condition in (3.23) is constant, and the evolution equation (3.22) admits the trivial solution $u(\xi,\tau)=0$ and $\eta(\xi,\tau)=Q^{*}\simeq 0.87349\,Q$ for all $\tau\in\hbox{{I\hskip-1.50696ptR}}$.

3.2.2 Shock unfolding for double-Stoker initial data

Next, and in order of increasing complexity, we look at the double-Stoker case (3.11), where the $(x,t)$-dependence of the boundary forcing terms $N$ and $V$ is known explicitly, and hence so is the functional form of $N(e^{\tau},t(\tau))$ and $V(e^{\tau},t(\tau))$.

We use asymptotics for $\tau\to-\infty$ in equations (3.28)-(3.31). Expanding the function $F$ in Taylor series as its argument goes to zero, or $e^{\tau}-\Phi(\tau)=O(e^{\tau})$ as $\tau\to-\infty$, yields at leading order

[TABLE]

where the constants $A$ and $B$ follow from the asymptotic limits as $\tau\to-\infty$

[TABLE]

Here we have used the asymptotic expression

[TABLE]

which follows from the leading order of the third equation in system (3.26). The same leading order evaluation of this equation yields the asymptotics for $\Phi$

[TABLE]

so that equation (3.34) provides a closed form relation to the parameters of the double-Stoker data for the derivative

[TABLE]

With this formula at hand, the asymptotic leading order behaviour of the boundary values $\eta(1,\tau)$ and $u(1,\tau)$ can be readily found,

[TABLE]

Note that the spatial dependence of $u(x,t)$ in the inner region between shocks x\in\big{(}-x_{s}(t),x_{s}(t)\big{)} at leading order as time $t\to 0^{+}$ can already be read off from the second equation in (3.37), since

[TABLE]

that is, the velocity field between shocks jumps from a vertical step of amplitude $2\displaystyle\sqrt{Q}$ at $t=0$ to a linear behaviour with finite (negative) slope

[TABLE]

The first equation in (3.37) shows that from the shock value $Q^{*}$ the free surface $\eta(0,t)$ initially evolves in the inner region (restoring the tilde notation to return to physical variables) with an $x$-independent base-point that increases linearly in time according to

[TABLE]

the decimal value of the slope $\mu^{(1)}$ for this linear increase in time is $\mu^{(1)}=-\nu^{(0)}Q^{*}\simeq 0.20762\,\sqrt{g_{0}}\,Q\,$. These scaling behaviours with respect to time and parameters are consistent with the general relations (3.32) and (3.33) that follow directly from the perturbed evolution equations (3.26). The choice of notation $\mu^{(1)}$ and $\nu^{(0)}$ in the above formulae is deliberate, since these parameters play a similar role as their counterparts for the self-similar solutions (2.2); however, we stress that in this case (3.38) and (3.40) provide asymptotic approximations, rather than exact solutions, for the evolution at short times $t>0$.

The above expansions can be continued to order $O(e^{2\tau})$, so that the asymptotic functional dependence on $\xi$ for $\eta$ and on $\tau$ for $u$, respectively, can be computed. A glance at the general solution (3.30) shows that the respective structure of these higher order corrections is that of quadratic and linear dependence on $\xi$. With the quadratic corrections the asymptotic form of $\tilde{\eta}$ and $u$ as $\tau\to-\infty$ becomes

[TABLE]

where the Taylor coefficient $F^{\prime\prime}(0)$ and that for the asymptotic expansion of $\Phi(\tau)$ (extending (3.35), where $\Phi_{0}=\displaystyle\sqrt{Q^{*}}/s_{0}$),

[TABLE]

at second order are

[TABLE]

The second order correction also enters the asympotics for the shock location $x_{s}$; from system (3.26), it can be shown that

[TABLE]

as $\tau\to-\infty$ (and hence $t\to 0^{+}$). We omit the somewhat lengthy details of the calculation for these expressions, but use these results to validate the numerical approaches in section 4 below. An alternative construction of the shock-continuation solution for this case is also sketched in the Appendix. In the original $\eta(x,t)$, $u(x,t)$ variables, the asymptotic expressions (3.41),(3.42) lead to

[TABLE]

We remark that, similarly to what happens at leading $O(e^{\tau})$-order for the slope $u_{x}(0,t)$ for this choice of initial data, the curvature $\eta_{xx}$ for the solution between shocks has a finite initial value at $t=0^{+}$, in this case set by the $O(e^{2\tau})$-order. That both the $u$-slope and the $\eta$-curvature turn out to be time-independent at their respective orders is a peculiarity of this double-Stoker problem, and their time evolution can be expected to be defined by higher order terms and later time evolution.

The construction of the solution outlined above is schematically illustrated in figure 6. As we shall see next, most of the qualitative features of this solution’s continuation past the shock time persist for the full case of initial data (2.17).

3.2.3 Shock unfolding for the full case

The result (2.39) of section 2.2 for the full case shows that the asymptotic behaviour as $t\to t_{c}^{+}$, which yields the corresponding asymptotics as $\tau\to-\infty$ of $N(e^{\tau},t(\tau))$, is

[TABLE]

where we have used the leading order for $t(\tau)$ from system (3.26), i.e., $t\sim x_{s}(t)/s_{0}$, in (2.39). The asymptotics for $\tau\to-\infty$ of equations (3.28)-(3.31) imply that the dependence of the function $F$ on its argument as this limits to zero must be through a functional composition with a $2/3$-power. Hence, the counterpart of (3.34) in the full case becomes the asymptotic relation

[TABLE]

where the coefficient $A$ for the full case, $A_{f}$ say, is now

[TABLE]

The leading order asymptotic for $\Phi$ (3.35) applies to the full case as well, hence the counterpart of equation (3.36) becomes

[TABLE]

The leading order asymptotics as $\tau\to-\infty$ determines the functional form of the perturbing fields $\tilde{\eta}(\xi,\tau)$ and $u(\xi,\tau)$,

[TABLE]

or, in the original variables and in the asymptotic limit $t\to t_{c}$, the solution between shocks is explicitly

[TABLE]

Here, we have restored the original (unshifted) time $t$ and elevation $\eta$ to further highlight the dynamics that develops shortly after the collapse time $t=t_{c}$. In particular, the velocity’s slope $u_{x}(0,t)$ and elevation’s curvature $\eta_{xx}(0,t)$ relax in time from their infinite values at $t=t_{c}$ ($u$ jumps from $u(0^{-},t_{c})=\displaystyle\sqrt{Q}$ to $u(0^{+},t_{c})=-\displaystyle\sqrt{Q}$, and $\eta(x,t_{c})$ has cusp singularity at $x=0$) with a monotonic increase from $-\infty$ which initially scales according to, respectively,

[TABLE]

Of course, these relations are consistent with those computed directly from the perturbed evolution equations (3.26) through (3.32) and (3.33), as can be easily verified.

The higher order correction to the shock position, which at leading order is $x_{s}\sim s_{0}\,t$, can now be determined. From the time advancing equation (3.26), the O\big{(}e^{{5\over 3}\tau}\big{)} terms are $t(\tau)\sim x_{s}/s_{0}+\hbox{const.}\,\,x_{s}^{5/3}$ so that

[TABLE]

i.e., the shock accelerates, as a function of the parameters $\gamma_{0}$ and $Q$, from the initial Riemann shock speed $s_{0}$ at a rate $O(t^{-1/3})$ as $t\to 0^{+}$. Note that unlike the double-Stoker case, the slope $u_{x}$ does not jump from (negative) infinity to a finite value at $t=t_{c}^{+}$, as the slope increase in the full case diverges with scaling $O(t-t_{c})^{-1/3}$.

4 Numerics

The theoretical results described in sections 2 and 3 can be illustrated by simulating numerically the evolution governed by the Airy model (2.1) for the class of initial data we have considered. In turn, the closed form solutions together with the asymptotic approximations we have derived can be used to monitor the accuracy of numerical algorithms after shocks develop.

We first consider the numerical approximation of the Airy’s system in the shock-fitted, unfolding coordinates (3.16). The resulting nonlinear wave equation is expected to always have regular solutions, thus avoiding the need for shock-capturing schemes. However, the numerical task for devising an algorithm is by no means trivial given the highly nonlinear nature of the coupled equations for the time advancement (3.18) and for the boundary condition (3.19) at $\xi=1$, which both drive and are determined by the solution of the wave-system (3.22). This structure is reminiscent of the classical water wave problem (see, e.g., [31]), in that the solution of a PDE in an interior domain is required to determine its boundary conditions, whose evolution in turn depends on the PDE solution.

Conversely, over the whole physical domain the original system of PDEs (2.1) with the initial data considered in Section 3 is expected to develop shocks even for continuous initial data. Implementing a numerical algorithm for a hyperbolic system of conservation laws up and beyond the shock formation while maintaining accuracy can also be a challenging task, especially when high precision is required, and it is worth describing our approach for this problem in some detail.

The first two parts of this section are devoted to the description of two approximation schemes, respectively for the shock-fitted variables $(\xi,\tau)$ and for the original system in (physical) variables $(x,t)$; the third part consists of a panoramic view of the numerical results, including a discussion on the validation of the numerical computations by the theoretical predictions.

4.1 Shock-fitted scheme

The Airy system in shock-fitted coordinates (3.16) consists of a couple of nonlinear wave equations set on the interval $[0,1]$. For a proper numerical treatment of the boundary conditions, it is convenient to rewrite the system in terms of the Riemann variables (3.25), which leads to the coupled equations

[TABLE]

The original variables can be expressed in terms of the Riemann invariants,

[TABLE]

and with this change of variables the boundary conditions (3.19) can be rewritten as

[TABLE]

where for simplicity we introduced the notation (and new variable) $R_{d}=R_{+}-R_{-}$. To ensure the well-posedness of problem (4.1), the boundary conditions must be imposed in the form of inflow conditions for the appropriate Riemann invariant. Since the variable $R_{+}$ has a positive velocity, the condition (4.3) at $\xi=0$ must be imposed on $R_{+}$; conversely, since $R_{-}$ has a negative propagation velocity, the condition (4.4) for $\xi=1$ must be imposed on $R_{-}$.

For the numerical approximation of system (4.1), we choose the Chebyshev Collocation Method (CCM) (briefly sketched below, see e.g., [5, 6] for a comprehensive exposition). The CCM is a special instance of a spectral method, meaning that the solution to the original PDE is expanded in a suitable set of basis functions (with an approximate solution then obtained by truncating the expansion). The truncated equations are the approximate analogue of the original PDE, so that the approximation error decreases exponentially fast with the number $M$ of basis functions retained, at least as long as the solution is sufficiently regular.

To set up the CCM, we introduce a set of collocation points $\{\xi_{i}\}_{i=1}^{M}\subset[0,1]$ that will be used both for the definition of the basis functions and in the determination of the expansion coefficients,

[TABLE]

The points $\{\xi_{i}\}$ defined above are often termed the Gauss–Lobatto–Chebyshev (GLC) points, to emphasize that the endpoints of the interval $[0,1]$ are included in the set of collocation points. It is now possible to introduce the set of Lagrange interpolant polynomials on the GLC points, as

[TABLE]

The CCM is based on writing the approximate solutions $R_{+}^{M}$, $R_{-}^{M}$ as a linear combination of the polynomials $T_{j}$,

[TABLE]

where the unknown coefficients $r_{j},s_{j}$ are determined by imposing that the original PDE (4.1) holds pointwise at the set of collocation points $\{\xi_{i}\}_{i=1}^{M}$

[TABLE]

for all $i=1,\dots,M$. Here primes denote derivatives with respect to either of the independent variables $\xi$ or $\tau$, as appropriate for each functional argument. Note that by virtue of (4.6), $T_{j}(\xi_{i})=\delta_{ji}$, and the first derivatives $T_{j}^{\prime}(\xi_{i})$ are known in closed form, see e.g. [5, p. 89]. Furthermore, due to the interpolant nature of the polynomials $T_{j}(\xi)$, the expansion coefficients $r_{j},s_{j}$ can be interpreted as the pointwise values $R_{+}(\xi_{j})$ and $R_{-}(\xi_{j})$, respectively. Note also that the nonlinearity of these equations is not confined to the quadratic terms, but enters through the coupling to the boundary conditions in the $\tau$ dependent factor $t^{\prime}(\tau)\,e^{-\tau}$, which depends on the Riemann invariants evaluated at $\xi=1$.

Equations (4.8) and (4.9) are a set of first order, nonlinear, coupled ordinary differential equations that govern the time evolution of the coefficients $\{r_{j}\}$ and $\{s_{j}\}$. In this work, we approach equations (4.8) and (4.9) by means of explicit, either Euler or fourth-order Runge–Kutta (RK4), methods. Thus, for instance, the Euler scheme time-advancing leads to the difference equations

[TABLE]

where we have adopted the shorthands $r_{j}^{n}=r_{j}(\tau_{n})$ (and similarly for $s_{j}$), $t_{n}=t(\tau_{n})$ and $\phi^{n}=\phi(e^{\tau_{n}},t_{n},\eta^{n}_{M},u^{n}_{M})$ for the values of the numerical approximations on the equispaced, constant $h$-duration subintervals $[0,\tau_{1},\dots,\tau_{n},\tau_{n+1},\dots]$.

Since both the forward Euler and the RK4 methods are explicit, a Courant-Friedrichs-Levy (CFL) condition holds on the time step size $h$. It is known from theoretical considerations that the velocity will achieve its maximum (in modulus) at $\xi=1$, and the collocation grid size achieves its minimum precisely at the same point, hence it is sufficient to check that the CFL condition is satisfied at the last grid point, $\xi_{M}=1$.

After advancing in time with the RK4 analogue of equations (4.10) and (4.11), the boundary condition (4.3) at $\xi=0$ can be imposed as follows:

[TABLE]

thus making sure that the relative inflow boundary condition is enforced. For the boundary condition at $\xi=1$, we found by trial and error that the accuracy in the approximation of the nonlinear relation (4.4) has a strong influence on the stability of the approximation method. To address this, we set up a Newton iterative scheme to make sure that, after propagating $r_{M}^{n+1}$ according to the RK4 analogue of equation (4.10), $s_{M}^{n+1}$ is chosen so that the nonlinear relation (4.4) holds up to machine precision accuracy.

Lastly, we mention that the unavoidable aliasing errors are controlled by guaranteeing a sufficiently high resolution in space, see e.g. [5, §3.10]. This approach towards de-aliasing, usually deemed impractical in two- or three-dimensional computations, is not unreachable in one-dimensional settings.

4.2 Shock-capturing scheme

The shock-capturing approximation scheme adopted for this work is a high-order, finite difference method commonly referred to as Weighted Essentially Non-Oscillatory (WENO). For an overview of the WENO method we refer to [28], and to [27] for a discussion on the details of a practical implementation.

Suppose that the Airy system (2.1) is set on the interval $I_{a}=[-c,c]$, and consider the set of $M$ ordered, equispaced, grid points $\{x_{i}\}_{i=0}^{M-1}\subset I_{c}$, chosen so that $x_{0}=-c$, $x_{M-1}=c$. For convenience, the distance between two consecutive grid points is assumed to be the same for all the grid points, and say $\delta=x_{i+1}-x_{i}$. Here we are concerned with finding an approximation for the pointwise values of $\eta$ and $u$ at the grid points.

For this purpose, it is convenient to rewrite equation (2.1) in conservative form,

[TABLE]

where, respectively, the conserved quantities $m=\eta\,u$ and $\eta$ can be interpreted physically as mass and momentum densities, and

[TABLE]

are the corresponding fluxes. In the following, we denote with $x_{i+1/2}$ the arithmetical average between the consecutive grid points $x_{i}$ and $x_{i+1}$, and with $x_{x-1/2}$ the arithmetical average between the consecutive grid points $x_{i-1}$ and $x_{i}$. Similarly, we denote with $m_{i}$, $\eta_{i}$ the pointwise values of $m$ and $\eta$ at the grid point $x_{i}$, and with $G_{i\pm 1/2}$, $L_{i\pm 1/2}$ the pointwise values of the fluxes $G$ and $L$ at $x\pm 1/2$.

The finite difference WENO scheme provides a numerical approximation to the solutions of system (4.14) through the ODE system

[TABLE]

where the space derivative of the fluxes is expressed formally as the difference quotient of some approximations $\widehat{G}$, $\widehat{L}$ of $G$ and $L$, evaluated at the two half-coordinate grid points centered in $x_{i}$. The time derivatives in equation (4.16) will be discretised with a finite difference method based on the knowledge of the values of $m$ and $\eta$ at the grid point $x_{i}$ at the previous time steps. The time advancement is done by means of a third-order, strong stability preserving, Runge–Kutta method, for the details of which we refer to [15].

In this work, we consider a WENO scheme with fifth-order accuracy in space, meaning that the flux approximation $\widehat{G}$ fulfills

[TABLE]

and a similar property holds for $\widehat{L}$. WENO schemes construct the approximations $\widehat{G}$, $\widehat{L}$ from the pointwise values of $G$ and $L$ on the grid $\{x_{i}\}_{i=0}^{M-1}$. For the Airy model, the expression for the fluxes is particularly simple, since it only involves quadratic nonlinearities, and it does not depend on the derivatives of $u$ and $\eta$.

To avoid spurious oscillations for the approximate fluxes, $\widehat{G}_{i\pm 1/2}$ and $\widehat{L}_{i\pm 1/2}$ are computed by means of Lagrange interpolation formulas based on three sets of points, namely:

[TABLE]

Then, by means of a local regularity estimator for $u$ and $\eta$, the approximate fluxes can be reconstructed using all of the three stencils in (4.18) if the solution does not have a discontinuity in the interval $[x_{i-1},x_{i+3}]$, thus giving a fifth-order approximation; otherwise, the approximate fluxes are reconstructed by neglecting the stencils in which the solution presents a shock, leading to an approximation which is locally of third order. For stability issues, it is convenient to split the fluxes as

[TABLE]

with the scheme adopted here being the Lax–Friedrichs splitting,

[TABLE]

where $\alpha$ is the maximum modulus among the eigenvalues of the Jacobian of the system of conservation laws. For the Airy system, this value can be computed explicitly as

[TABLE]

4.3 Numerical results

Among various alternatives for illustration and validation of the numerical schemes’ performance, it is particularly interesting to focus on the full case and monitor quantities along the $x=0$ symmetry line, such as the short time evolution of $\eta$, $u_{x}$, $\eta_{xx}$, and consider pointwise differences between solutions by algorithms and their closed form analytical counterparts when available. We focus mainly on the “full” case both because of its interest in the theoretical question of the wetting of a dry spot, and because of its challenging non-analytic behaviour for the continuation past the collapse point. Nonetheless, the analogous results about the “double-Riemann” and the “double-Stoker” problems will be discussed briefly to illustrate how the numerics performs with respect to the increasing level of complexity at the collapse point. Unless otherwise specified, all spectral computations are performed with $M=64$ modes and timestep $h=10^{-4}$, while for the shock-capturing scheme the grid will vary, depending on a target accuracy, from $M=2^{14}$ to $M=2^{18}$ with timesteps from $10^{-8}$ to $10^{-6}$.

4.3.1 Double-Riemann simulations

As already remarked, the dynamics for this case in shock fitting coordinates is trivial and does not pose a challenge to the spectral code implementation. On the other hand, the WENO algorithm is designed to deal with problem like these, and so this is an excellent test for this scheme, allowing to check pointwise accuracy and sharpness in resolving discontinuities in both space and time. In figure 7, we show the results of a Riemann problem with constant initial density $\eta_{0}=Q/2$ and discontinuous initial momentum $m=-{\rm sgn}(x)\displaystyle\sqrt{Q}(Q/2)$. For the numerical test we take $Q=1/2$, and run the code on the interval $[-1,1]$ for several space and time resolutions. The inflow conditions are the same as the initial conditions, namely $\eta=Q/2$ and $m=\mp\displaystyle\sqrt{Q}(Q/2)$ are imposed respectively at $x=\pm 1$. The results show that the code is indeed high order at a sufficient distance from the discontinuity, as the height of $\eta$ after the collision matches the theoretical prediction to 5 digits even for a modest resolution of $2^{10}$ grid points and $h=10^{-4}$. However, by focusing on the plot for $\eta(x=0,t)$, the numerical solution after shock formation shows a significant overshoot followed by a region of damped oscillations. While the temporal support of these oscillations decreases with increasing resolution, their amplitude does not. We can safely interpret these oscillations as numerical artifacts, as their frequency is significantly smaller than the Nyquist frequency for all the timesteps considered. These oscillations also make estimating the shock position and velocity inherently imprecise, as the discontinuity’s support jumps erratically from one timestep to the next over a set of a few spatial grid points.

4.3.2 Double-Stoker simulations

The purpose of this section is to check the consistency of the theoretical results of section 3.1.2 with the numerical computations. More precisely, we consider the initial datum (3.11) with $Q=1/2$ and $g_{0}=1/16$, and aim at checking the asymptotics of equations (3.41) and (3.42).

To avoid the proliferation of figures, we consider only the shock-fitted spectral method, referring to Section 4.3.3 for a complete comparison between the shock-capturing and the shock-fitted schemes for the full case. The results that we report here are the values of $\tilde{\eta}$, $\tilde{\eta}_{\xi\xi}$ and ${u}_{\xi}$ evaluated at $\xi=0$ as a function of $\tau$ (figure 8, first row); the analogous values of $\tilde{\eta}$, $\tilde{\eta}_{xx}$ and ${u}_{x}$, evaluated at $x=0$ as a function of $t$ (figure 8, bottom row); the profiles of $\tilde{\eta}$ and $u$ as a function of $\xi$, for a few representative values of $\tau$.

It can be seen that in all the cases the agreement with the theoretical predictions is quite reasonable. The agreement between asymptotic results and numerical computations is to some extent less clear for the curvature of $\eta$ seen in the original variables $x,t$. We attribute this effect to the fact that for computing $\tilde{\eta}_{xx}$ starting from $\tilde{\eta}_{\xi\xi}$, it is necessary to compute a ratio between $\tilde{\eta}_{\xi\xi}$ and $e^{2\tau}$. Both these quantities are very small for small times, making the resulting ratio an ill-conditioned problem. Such cancellation errors fade away as $\tau$ and $\tilde{\eta}_{\xi\xi}$ grow, and this is confirmed by the good fit for times $t$ between $0.01$ and $0.04$.

Finally, we check the spectral code with the analytic expression in equations (3.41) and (3.42) for the profiles of $\eta$ and $u$ in shock-fitted coordinates. The results are summarized by figure 9, which shows how the theoretical curves are indistinguishable from the numerical computed solutions within plotting resolution.

4.3.3 Full case simulations

For the shock-fitted spectral method, the initial datum corresponds to an asymptotic condition for $\tau\to-\infty$. We approximate the asymptotic condition by taking $x_{s}=10^{-8}$, and a corresponding initial value of $\tau=\log x_{s}\simeq-18.4207$. Initially, the Riemann variables are set to their asymptotic values of $R_{\pm}=\pm\displaystyle\sqrt{Q^{*}}$. There is clearly an inconsistency arising from the use of the asymptotic values for the Riemann variables while the initial value for $\tau$ is a finite number. A result of this inconsistency is the generation of spurious, small amplitude, waves in the initial instants of the numerical computations. These “acoustic” waves propagate in the computational domain and are reflected at the boundaries. Due to the absence of numerical diffusivity in Chebyshev spectral methods, the acoustic waves are dissipated only by interaction with the shock wave at $\xi=1$. A side-effect of these spurious waves is that the numerical results can be trusted only after a few units of $\tau$. However we verified that the dissipation of the acoustic waves is sufficiently fast that the small time asymptotics can still be easily detected. Being interested in the short time asymptotics, the boundary conditions at $\xi=1$ are further simplified by adopting the expansion (2.39) for $N(e^{\tau},t(\tau))$.

We first test the spectral code against the asymptotic expressions in equations (3.52) and (3.53) for the profile of $\eta$ and $u$ with respect to the shock-fitted coordinates $(\xi,\tau)$. Figure 10 shows time snapshots of the profiles, and the good agreement between theory and numerics can be clearly seen.

Next, the two numerical methods considered here are compared with the theoretical results in figure 11. The purpose of this plot is to check the predicted growth rates of $\eta$, $u_{x}$ and $\eta_{xx}$ along the centerline for small times (here we perform the comparison up to $t=0.04$). The plots generally show good agreement of both numerical schemes with the theoretical rates for $\eta$ and $u_{x}$, with relative errors always below 10%; for $\eta_{xx}$ the method based on shock-fitted variables still provides accurate results, while for the shock-capturing method the error is substantially higher, being on average about 80%. A further drawback of the shock-capturing algorithm shows up as oscillations appearing just after the shock formation (see figure 7). These numerical artifacts are to be contrasted with the small amplitude acoustic waves of the shock-fitted spectral method, which perturb the numerical results for a much smaller (in fact, exponentially smaller) period of time.

Lastly, we check the prediction of equation (3.51) for the scaling of $\eta$ and $u$ with the initial condition parameters $\gamma_{0}$ and $Q$. This comparison is more conveniently done with the numerical scheme based on shock-fitted coordinates and the corresponding spectral code. The value of $\eta(\xi=0,\tau=0)$ is computed by extrapolating to $\tau=0$ the value $\tilde{\eta}(\xi=0,\tau)$ obtained numerically, in this case with $\tau=-6$. This provides the following approximation for $\tilde{\eta}(0,0)$,

[TABLE]

Next, evaluating equation (3.52) at $\xi=0,\tau=0$ gives

[TABLE]

and replacing in equation (4.23) the value computed numerically for $\eta$, it is possible to obtain an estimate for $F^{\prime}(0)$. The computations described above are repeated for several values of $\gamma_{0}$ and $Q$ ranging between 0.05 and 1. For a clearer comparison, we show in figure 12 (left) the values of $\tilde{\eta}(0,0)\gamma_{0}^{-1/3}$, so that the resulting data should depend only on $Q$, and not on $\gamma_{0}$. For analogous reasons, on the right plot of figure 12 we show the values of $\tilde{\eta}(0,0)Q^{-2/3}$ as a function of $\gamma_{0}$. The numerical results agree with the theoretical predictions, confirming the numerical coefficient of $F^{\prime}(0)$ to within three decimal digits.

We conclude the presentation of the numerical computations by showing the results of the shock-capturing (WENO) scheme. The shock-capturing scheme is set on the interval $[-\sqrt{3/2},\sqrt{3/2}]$, with $M=2^{18}$ grid points, and a time step $h=2\cdot 10^{-7}$. The initial datum in this case is obtained by solving for the implicit expression (2.29) at the grid point coordinates using the Mathematica package. In figure 13 we show an example of the numerical evolution using this approach and also compare the WENO evolution with the shock-fitted spectral code, mapped into $(x,t)$ coordinates. The profiles of $\eta$ and $u$ by the algorithms coincide up to three decimal digits, and indeed the results obtained with the two methods can hardly be distinguished. However, it should be stressed that the WENO evolution is affected by numerical uncertainty, similar to that occurring for the Riemann problem shown in figure 7 and displayed here in figure 14, for the time of shock formation. There is also a slight difference in location of the endpoints of the inner solution between shocks in each snapshot, which follows from the relative lack of accuracy of the shock-capturing WENO code in detecting the shock position. This latter quantity is smeared over a few grid points in the shock-capturing method, making it less precise as far as the shock position and its velocity are concerned.

5 Discussion and conclusions

In this work we have investigated the evolution of a class of exact solutions of the Airy shallow water model. These correspond to a particular choice of initial data with zero initial velocity, and surface elevation in the form of a downward parabola sandwiched between, and joined continuously to, constant flat states $\eta=Q$. When the parabola contacts the bottom in a “dry point,” the evolution of the elevation $\eta$ and velocity $u$ is markedly different than their counterpart in the “wet” case, where the interface is at a finite, however small, distance from the bottom. In fact, the dry point persists until the time when strong shocks develop, corresponding to the collapse of the initial parabola to a segment of height $Q/4$. In contrast, for the wet case, the parabola’s minimum eventually lifts up to the background elevation $Q$, and a zero velocity, constant elevation state $\eta=Q$ develops in a neighbourhood of where the minimum was originally located. The dry case parabola collapse can be viewed as a mechanism for the fluid to “fill the hole” that goes all the way to bottom, thereby allowing the reconnection of the fluid domain.

It is remarkable that both the dynamics prior to reconnection, and that ensuing afterwards, can be followed entirely by analytical means, either in exact form, for pre-collapse times, or asymptotically, for short times after collapse. Until the collapse time, all the evolution details can be followed with closed form exact solutions. The self-similar form of the parabolic part can be followed by solving a set of ODEs for its coefficients (for the dry point case, these are simply the curvature $\eta_{xx}$ at the minimum $x=0$ and the slope $u_{x}$ at the same point), while the connection to the background rest state is achieved by simple waves whose analytic form is expressed implicitly through an appropriate characteristic parameter. The continuation past collapse and strong shock formation is more delicate and requires proper shock fitting. We have achieved this by introducing stretched coordinates to magnify the region where the shock first appears. This tool borrows ideas from singular perturbation theory, and allows the analysis to proceed in the form of asymptotics for $\eta$ and $u$ between shocks, thereby continuing the solution in the limit of short times beyond that of collapse.

In order to illustrate, and provide intuition on, how this can be implemented, we constructed two ad hoc classes of initial data that capture the essential mechanisms of shock formation and propagation, but avoid the technical difficulties posed by the original case’s branch singularities at the collapse location. The first of these two classes, which we dubbed “double-Riemann,” is a variant of the classic shock problem for piecewise constant data. The second case, dubbed “double-Stoker,” is obtained by splicing together two mirroring simple wave solutions of the dam-break type illustrated in Stoker’s book [30]. This increases the level of complexity with respect to the double-Riemann case, and gets qualitatively closer to the “full” case by adding the fundamental component of a time- and space-dependent medium through which the shocks must propagate.

In all cases, the elevation $\eta$ at the origin jumps from the level $Q/4$ (of course, this particular choice is dictated by the pre-collapse evolution for the full case, and it is not arbitrary as in the ad-hoc constructions) at $t=t_{c}^{-}$ to a new level $Q^{*}\simeq 0.87\,Q$ at $t=t_{c}^{+}$, in the form of a thin jet between shock locations at $x=0$. The shocks propagate outward, with nontrivial dynamics, unlike the double-Riemann problem where the shock speed is constant, thickening the jet while its maximum grows with developing curvature. These are signature features of the dry point case with respect to its wet counterpart, where the background state with a flat interface $\eta=Q$ and $u=0$ is always fully recovered, for a finite time dependent on the initial minimum thickness of the fluid layer, in a neighbourhood of the origin. In this respect, a question that deserves to be investigated is whether this full recovery is a consequence of discontinuities in the derivatives of our initial conditions, and whether the possible generation of more generic singularities emerging at earlier times, and travelling towards the vacuum regions, can significantly alter the qualitative features of the evolution we have identified. However, these details would presumably be heavily dependent on the nature of the initial conditions; for instance, the corners at $\eta=Q$ could in principle be smoothed out with $C^{\infty}$ compactly supported bump functions, so that the core parabolic solution and the splicing point dynamics $a(t)$ could be maintained with the resulting smooth initial data. A study along these lines is undergoing and will be reported in the future.

With the analytical results at hand, we have implemented algorithms for the numerical solution of the Airy model, including the continuation past shock formation, and gauged their performance with respect to the theory. Not unexpectedly, the WENO scheme is able to give a qualitatively accurate global perspective in the physical coordinates. However, our study shows that a quantitatively accurate solution is more problematic for this scheme, at least on short time scales after collapse. Thus, for instance, the scheme is plagued at short times by oscillations that cannot be eliminated even pushing grid refinement to levels where computational costs become unreasonable. As a result of this, the time-scaling of monitoring quantities such slope and curvature of $u$ and $\eta$ at the origin can not be pinned down accurately, making the validation against the short time theory impractical. Further, while these oscillations die out, presumably due to dissipation at the shocks, the longer time evolution seems to be affected by the initial errors accumulated by the scheme, though for the double-Riemann case, where the shock propagates in a constant medium, the correct (theoretical) values are eventually recovered. However, even in this simplest case the shock elevation $Q^{*}$ is reached after an initial spiky overshoot that persists in magnitude while its time-support decreases as the space grid is refined and the time step is reduced. We remark that these oscillations might be of different nature than the numerical ones the WENO scheme is designed to suppress. The implementation of the unfolding coordinate scheme, and its perturbation formulation, in fact reveals that left- and right-propagating waves in the interval $\xi\in[-1,1]$ can be caused by a (small) compatibility mismatch between the initial (as $\tau\to-\infty$) and boundary conditions. This bouncing back and forth modes in the interval $\xi\in[-1,1]$ can persist until dissipated at the shocks, and, at least initially, can appear as high frequency oscillations in the compressed interval $[-x_{s},x_{s}]$ of the physical space. It must be noted, however, that while the numerical algorithm we have developed for the unfolding coordinate is far more accurate than WENO’s, as attested by comparison with the theory (also a consequence of the local solution remaining regular in these coordinates, allowing our implementation to be spectrally accurate), numerical round-off errors are nonetheless magnified when converting back to physical coordinates, owing to the multiplication of small floating point numbers for the solution by large, $O(e^{-\tau})$ for negative $\tau$, values. A thorough investigation to identify optimal algorithms able to handle the kind of evolutions that we have studied analytically is interesting from a numerical perspective, but would need a dedicated study outside of the scope of this work.

The theoretical framework developed here for shocks emanating from the class of initial conditions we have focussed on may be extended to more general set-ups, and possibly be generalized to a broad class of hyperbolic conservation laws. We plan to pursue this in future work, where we will also examine the interplay between dissipation at the shock location vs. its realization via the boundary condition in the shock fitted domain. A relevant goal, which we touched upon above and will be reported more in depth elsewhere, is that of extending our results for the class of initial data we studied to more general, and in particular smooth, initial data. Furthermore, a thorough investigation of the differences between initial data that are limiting sequences to those with dry spots is being carried out and will be reported in the future. We remark here that such a study would bear some similarity with the investigation on the formation of “splash” type singularities in models of free surface flows [7], and it would be interesting to explore this direction further.

We also remark that the reflection symmetry of the solutions we have studied, together with the assumption of negligible viscosity for the fluid applications, allows us to interpret our results as equivalent to the case of a vertical impermeable boundary condition, a “wall,” at $x=0$. Thus, for the full case, the initial bottom contact of the interface at $x=0$ would imply that the upward jet resulting from the impact of the water layer with the vertical wall would be initially more than $10\%$ lower than it would be if a wetting liquid film were maintained over the bottom. Of course, for such a simplified model this estimate cannot be expected to be quantitatively accurate for a real physical situation, nonetheless it isolates a trend that could be significant enough not to be swamped by other real-fluid properties, such as viscosity and surface tension, and possibly be observable in experiments.

Finally, we stress that the solution features uncovered by our analysis and simulations cannot obviously be expected to represent in all details the physics of fluids modelled by the Airy’s system, based as this is on a long-wave asymptotics of the Euler equations, which would be clearly violated by shock formation. However, by enforcing the fundamental conservation laws of mass and momentum, this system is known to be quite robust in capturing general features of the dynamics of shallow water, with shocks representing hydraulic jumps. Initial conditions of the kind studied here, while hard to implement in actual experiments, can be studied in silico, and in future work we plan to present results for the parent Euler equations governing ideal fluids. Generalizations of our class of initial conditions and a thorough study of the energy dissipation due to shock formation are also interesting new developments, which will be reported in future work.

Acknowledgments

RC thanks the hospitality of the ICERM’s program “Singularities and Waves in Incompressible Fluids” in the Spring of 2017, and the support by the National Science Foundation under grants RTG DMS-0943851, CMG ARC-1025523, DMS-1009750, DMS-1517879, and by the Office of Naval Research under grants N00014-18-1-2490 and DURIP N00014-12-1-0749. All authors gratefully acknowledge the auspices of the GNFM Section of INdAM under which part of this work was carried out. R.C., M.P. and G.P. thank the Dipartimento di Matematica e Applicazioni of Università Milano-Bicocca for its hospitality, G.F., G.O., M.P. and G.P. thank the Carolina Center for Interdisciplinary Applied Mathematics at the University of North Carolina for hosting their visits in 2018. Support by grant H2020-MSCA-RISE-2017 Project No. 778010 IPaDEGAN is also gratefully acknowledged.

Appendix: Shock continuation via local Taylor series

We seek a local solution of system (2.1) as $|x|\to 0$ by an expansion in powers of $x$,

[TABLE]

where the choice of even and odd powers respectively for $\eta$ and $u$ enforces the symmetry $\eta(-x,\cdot)=\eta(x,\cdot)$, $u(-x,\cdot)=-u(x,\cdot)$. In particular, as well known, when the initial data are analytic functions of $x$, the Cauchy-Kovalevskaya (C-K) Theorem (see, e.g., [14]) assures that these power series have a finite radius of convergence, so that all coefficients are analytic in a neighbourhood of $t=0$. The exact solutions (2.2) are of course in this class, and are special in that they derive from a hierarchy that truncates at the leading orders, $n=0,1$, so that $\eta$ and $u$ are second and first degree polynomials, respectively. Of course, this property is special, and the series (A1) cannot in general be expected to truncate to a polynomial form; once substituted in (2.1) the series generate a recursive infinite hierarchy of ODE’s for the coefficients $\eta_{n}(t)$ and $u_{n}(t)$. The first equations in the hierarchy are (cf. (2.3))

[TABLE]

Note that, when applied to zero velocity smooth initial data with a dry point $\eta_{0}(0)=0$, then $\eta_{0}(t)=0$ for all times and these equations reduce to a closed system for $u_{0}$ and $\eta_{1}$ which coincides with (2.8). In fact, the presence of the dry point has important implications, since when $\eta_{0}(t)=0$ the equations in the hierarchy for the pairs $\eta_{n+1}$ and $u_{n}$ are linear , with coefficients and inhomogeneous terms that depend solely on the ($\eta_{k+1},u_{k}$) pairs, $k<n$. The recursive solution for all coefficients thus shows that time singularities are determined locally in an neighbourhood of $x=0$ solely by the nonlinear system for the first pair $(\eta_{1},u_{0})$, i.e., equations (A2) with $\eta_{0}=0$.

While the C-K Theorem refers to smooth initial data, we assume next that a local solution in a neighbourhood of the origin away from shocks can still be expressed in the form (A1), for suitable initial conditions such as the double-Stoker case §3.1.2, and that equations (A2) govern the leading order approximation to the solution between shocks for short times $t>0$. Accordingly, we seek the coefficients $u_{0}(t),\eta_{0}(t)$, and $\eta_{1}(t)$ as Taylor expansions in time,

[TABLE]

where, as mentioned in §3.2.2, we use this notation to remind the reader of the role played by the analogous coefficients in the exact counterpart from initial data (2.17) for system (2.1). With the above assumptions, the shock position can be similarly expanded as $t\to 0^{+}$,

[TABLE]

As in section 3, the shock speed and the consistency condition are given by (3.1) and (3.2), respectively. In this case $f_{+}$ corresponds to the rarefaction (simple wave) solution (3.11) while $f_{-}$ is determined by the local solution (A1).

Carrying out the above expansions through the recursive steps, we obtain for the first few coefficients the following decimal values

[TABLE]

where the ordering in this list reflects that of the recursion method for finding the asymptotic series coefficients, beginning with the solution of the cubic equation that sets the shock amplitude, $\mu^{(0)}$, at time $t=0^{+}$. These are identical to the corresponding values computed in §3.2.2, however we remark that higher order coefficients lead to calculations that can be even lengthier than their counterpart in the unfolding coordinate approach.