Derivation of a viscous Boussinesq system for surface water waves

Abstract

In this article, we derive a viscous Boussinesq system for surface water waves from Navier–Stokes equations for non-vanishing initial conditions. The relevance of such initial conditions is proved. We use neither the irrotationality assumption, nor the Zakharov–Craig–Sulem formulation. During the derivation, we find the bottom shear stress and also the decay rate for shallow water. In order to justify our derivation, we derive the viscous Korteweg–de Vries equation from our viscous Boussinesq system and compare it to the ones found in the bibliography. We also extend the system to the 3D geometry.

Keywords

water waves shallow water Boussinesq system viscosity KdV equation

1. Introduction

1.1. Motivation

The propagation of water waves over a fluid is a long run issue in mathematics, fluid mechanics, hydrogeology, coastal engineering, … In the case of an inviscid fluid, the topic stemmed many researches and even broadened with time. Various equations have been proposed to model this propagation of water waves. Since the full problem is very complex, the goal is to find reduced models on simplified domains with as little fields as possible, should they be valid only in an asymptotic regime.

This article is a step forward in the direction of a rigorous derivation of an asymptotic system for surface water waves in the so-called Boussinesq regime, taking into account the viscosity. While viscous effects can be neglected for most oceanic situations, they cannot be excluded for surface waves in relatively shallow channels.

In the inviscid potential case, the complete and rigorous justification of most asymptotic models for water waves has been thoroughly carried out and summarized in the book [16] and the bibliography therein. This book includes the proof of the consistency and stability of some models, the proof of the existence of solutions both of the water waves systems and of the asymptotic models on the relevant time scales and the proof of “optimal” error estimates between these two solutions. The curlfree assumption allows to use the Zakharov–Craig–Sulem formulation of the water waves system and facilitates the rigorous derivation of the models, through expansions of the Dirichlet to Neumann operator with respect to a suitable small parameter.

Things are more delicate when viscosity is taken into account and a complete justification of the asymptotic models is still lacking. The main difficulty, for not only a derivation but also for a rigorous proof, arises from the matching between the boundary layer solution coming from the bottom and the “Euler” solution in the upper part of the flow.

In this article, we derive an asymptotic system (Boussinesq system) for the viscous flow in a flat channel of water waves in the Boussinesq regime, that is to say in the long wave, small amplitude regime with an ad hoc balance between the two effects.

1.2. Literature

When deriving models of water waves in a channel and taking viscosity into account, numerous pitfalls must be avoided in order to be rigorous.

Since there are various dimensionless parameters, a linear study must be done so as to determine the most interesting regime between the parameters. One must also either assume linearized Navier–Stokes equation (NSE), or justify that the nonlinear terms can be dropped. This is not so obvious because numerous authors extend the inviscid theory by assuming the velocity to be the sum of an inviscid velocity and a viscous one. Then they force only one condition (for instance the vanishing velocity on the bottom) to be satisfied by the total velocity, once the inviscid velocity is assumed unchanged by viscosity. This assumption deserves to be justified.

At a certain level of the derivation, a heat-like equation arises. Most people solve it with a time Fourier transform while the only physical problem is a Cauchy one, so with an initial condition. The only possibility is to use either Laplace (in time) transform or a sine-transform (in the vertical dimension) with a complete treatment of the initial condition.

One must also derive the bottom shear stress because it is meaningful for the physicists who deal with sediment transport.

Last, the order up to which the expansion is done must be consistent throughout the article.

To the best of our knowledge, no article does all this. Yet various articles have been written on this topic. Let us review those that retained our attention and interest.

Boussinesq did take viscosity into account in 1895 [2]. Lamb [15] also derived the decay rate of the linear wave amplitude by a dissipation calculation (done in paragraph 348 of the sixth edition of [15]) and by a direct calculation based on the linearized NSE (paragraph 349 of the sixth edition of [15]). Both of them used linearized NSE on deep-water and computed the dispersion relation. We do not know who is the first. The imaginary part of the phase velocity gave the decay rate: $\frac{\partial A}{\partial t} = - 2 ν k^{2} A,$ (1) where A is the amplitude of the wave, ν the kinematic viscosity and k the wavenumber.

In [24], Ott and Sudan made a formal derivation (in nine lines) of a dissipative KdV equation (different from ours). They used the linear damping of shallow water waves already given by Landau–Lifshitz. This led them to an additional term to KdV, which looks like a half integral. They also found the damping in time of a solitary wave over a finite depth as ${(1 + T)}^{- 4}$ (already found by [13] and later by [11,12,22] (p. 374)).

Byatt-Smith studied the effect of a laminar viscosity (in the boundary layer where a laminar flow takes place) on the solution of an undular bore [3]. He found the (almost exact) Boussinesq system of evolution with a half differentiation but with no treatment of the initial condition. He did an error when providing the solution to the heat equation: his convolution in time is over $(0, + \infty)$ instead of $R$ (Fourier convolution) or $(0, t)$ (Laplace convolution).

In 1975, Kakutani and Matsuuchi [12] found a minor error in the computation of [24]. They started from the NSE and performed a clean boundary layer analysis. First, they made a linear analysis that gave the dispersion relation and, under some assumptions, the phase velocity as a function of both the wavenumber of the wave and the Reynold’s number Re. They distinguished various regimes of Re as a function of the classical small parameter of the Boussinesq regime. Then, they derived the corresponding viscous KdV equation. We want to stress that, at the level of the heat equation, they used a time Fourier transform. As a consequence, they may not use any initial condition. So, the problem they solve is not the Cauchy’s one.

In [21], one of the authors of the previous article [12] tried to validate the equation they had been led to. He showed that their “modified K-dV equation agrees with Zabushy–Galvin’s experiment with respect to the damping of solitary waves, while it produces disagreement in their phases” (see the conclusions). One might object that the numerical treatment seems light because the space step was between one and 10 percent, the numerical relaxation was not very efficient, the unbounded domain was replaced by a periodic one though there is “non-locality of the viscous effect” (p. 685), there was no numerical validation of the full algorithm, and the regime was not the Boussinesq one (the dispersion’s coefficient was about 0.002 and the viscous coefficient was 0.03). Moreover, the phase shift numerically measured was given with three digits while the space step was of the order of magnitude of some percents. The author, very fairly, added that “the phase shift obtained by the calculations is not confirmed by [the] experiments”.

In 1987, Khabakhpashev [14] extended the derivation of the viscous KdV evolution equation to the derivation of a viscous Boussinesq system. He studied the dispersion relation and predicted a reverse flow in the bottom, in case of the propagation of a soliton wave. He used a Laplace transform (instead of Fourier as [12] did) with vanishing initial conditions since he assumed starting from rest. Although he stressed this assumption, he acknowledged that “the time required for the boundary layer to develop over the entire thickness of the fluid [is] much greater than the characteristic time of the wave process”. The equations were not made dimensionless, so the right regime was not discussed and a very inefficient numerical method was used (Taylor series expansion is replaced in the convolution term).

In the book [11] (Part 5, pp. 356–391), Johnson found the same dispersion relation as [12], studied the attenuation of the solitary wave by a multiscale derivation, reached a heat equation, but solved it only with vanishing initial condition. He exhibited a convolution with a square root integrated on $(0, + \infty)$ (like Byatt-Smith [3]). Some numerical simulations (already partially done by [3]) enabled him to recover the mechanism of undular bore slightly damped.

Later, Liu and Orfila wrote a seminal article [18] (LO hereafter) in which they studied water waves in an infinite channel (so without meniscus). They derived a Boussinesq system with an additional half integration (seen as a convolution), and an initial condition assumed to be vanishing, but implicitely added to the system when numerical simulation must be done.

More precisely, the authors took a linearized Navier–Stokes fluid, used the Helmholtz–Leray decomposition and defined the parameters (index $LO$ denotes their parameters): $\begin{matrix} α_{LO}^{2} = ν / (l \sqrt{g h_{0}}), \\ ε_{LO} = A / h_{0}, \\ μ_{LO} = h_{0} / l, \end{matrix}$ where the following notations will be used throughout the present article: A is the characteristic amplitude of the wave, $h_{0}$ is the mean height of the channel, g is the gravitational acceleration and l is the characteristic wavelength of the wave. They made expansions up to order $α_{LO}$ , which square is a kind of a Reynold’s number inverse. They used the classical Boussinesq approximation: $ε_{LO} ≃ μ_{LO}^{2}$ , but they also set the link between the viscosity and $ε_{LO}$ by requiring $O (α_{LO}) ≃ O (ε_{LO}^{2}) ≃ O (μ_{LO}^{4})$ without further justification. Although “the boundary layer thickness is of $O (α_{LO})$ ”, they stretched the coordinates by a larger factor $α_{LO} / μ_{LO} ≃ μ_{LO}^{3}$ (see their (2.9)). More important, and maybe linked, they kept the $α_{LO} μ_{LO}$ terms (in their (2.8) or (2.21) for instance) and yet dropped $o (α_{LO})$ terms! This can explain why their final system (3.10)–(3.11) had a $α_{LO} / μ_{LO} = O (ε_{LO}^{3 / 2})$ term before the half integration, while we will justify an $O (ε_{LO})$ term for our system.

Let us stress that assuming $α_{LO}^{2} = ε_{LO}^{4}$ as did [18] amounts to $Re = ε_{LO}^{- 7 / 2}$ with our (further redefined) Reynold’s number: $Re = ν / (h_{0} \sqrt{g h_{0}})$ , while we prove below that the regime at which gravity and viscosity are both relevant is $Re = ε_{LO}^{- 5 / 2}$ . Our regime was also exhibited by [3,11,12]. So, [18] studied a regime different from ours.

Last, they claimed the shear stress at the bottom to be: $τ_{bottom} = - \frac{1}{2 \sqrt{π}} \int_{0}^{t} \frac{u (x, T)}{\sqrt{{(t - T)}^{3}}} d T,$ where $u (x, T)$ is the depth averaged horizontal velocity. Indeed this integral is infinite as they acknowledged in a later corrigendum where they claimed the right formula to be: $τ_{bottom} = \frac{1}{\sqrt{π}} \frac{u (x, 0)}{\sqrt{t}} + \frac{1}{\sqrt{π}} \int_{0}^{t} \frac{u_{T}}{\sqrt{t - T}} d T,$ where $u_{T}$ denotes the time derivative. However, they did not provide any justification. Moreover, their solution (2.15) to the heat equation, computed in [23] (pp. 153–159), assumed vanishing initial condition. So the treatment of the initial condition was not done. One of our goals in the present article is precisely to provide a better treatment of this initial condition. In this article [18], the authors also raised the question of the eligible boundary condition. Indeed, they considered to be well known that for a laminar boundary layer, the phase shift between the bottom shear stress and the free stream velocities is $π / 4$ . So it prohibits any bottom condition of the Navier type $τ_{x y} = - k u_{bottom}$ as is sometimes assumed (and not derived).

Although we presented some criticisms, we acknowledge the modeling, derivation and explanations of this article are insightful and, last but not least, very well written. Yet our criticisms apply to all subsequent articles of the same vein.

In [19], Liu et al. experimentally validated LO’s equations in the particular case of a solitary wave over a boundary layer. By Particle Image Velocimetry (PIV), they measured the horizontal velocity in the boundary layer over which the solitary wave run and confirmed the theory.

In [20], Liu et al. extended the derivation of the viscous Boussinesq system of [18] to the case of an uneven bottom. They compared the viscous damping and shoaling of a solitary wave, propagating in a wave tank from the experimental and numerical point of view. They provided a condition on the slope of the bottom and paid attention to the meniscus on the sidewalls of the rectangular cross section.

In [17], Liu and Chan used the same process to study the flow of an inviscid fluid over a mud bed modeled by a very viscous fluid. They also studied the damping rate of progressive linear waves and solitary waves. In [25], Park et al. validated this model with experiments. They also studied the influence of the ratio of the “mud bed thickness and the wave-induced boundary-layer thickness in the mud bed”.

In 2008, Dias et al. [7] took the (linearized) NSE of a deep water flow with a free boundary and used the Helmholtz–Leray decomposition. Both Bernoulli’s equation, through an irrotational pressure, and the kinematic boundary condition were modified. Then they made an ad hoc modeling for the nonlinear term. Starting from such a model, they provided the evolution equation for the envelope A of a Stokes wavetrain which, in case of an inviscid fluid, is the Non-Linear Schrödinger (NLS) equation. The provided equation happens to be a commonly used dissipative generalization of NLS.

Although it was published earlier (2007), [10] was a further development of [7] to a finite-depth flow. In this article, the authors still linearized NSE and generalized by including additional nonlinear terms.

In a later article [9], Dutykh linearized NSE and worked on dimensioned equations, considering the viscosity ν to be small (in absolute value). The author generalized by “including nonlinear terms” and reached a viscous Boussinesq system (his (11)–(12)). Making this system dimensionless triggered very odd terms and its zeroth order was no longer the wave equation. He further derived a KdV equation by making a change of variable in space (but not the associated change in time $τ = ε t$ ). He also made a study of the dispersion relation by assuming an exponential function ansatz of the type $e^{i (k x - ω t)}$ , but then he froze the half derivative term. Indeed it is well known that such an ansatz amounts to make a Fourier or Laplace transform. Here, the Fourier/Laplace transform of the half derivative of u is very simple: $| ξ |^{1 / 2} \hat{u}$ and could have been used instead of freezing this half derivative term.

In [4], Chen et al. investigated the well-posedness and decay rate of solutions to a viscous KdV equation which has a nonlocal term that is the same as Liu and Orfila’s [18] and [9], but not the same as [12] nor the same as ours. The theoretical proofs were made with no dispersive term ( $u_{x x x}$ ), but with a local dissipative term ( $u_{x x}$ ). The tools were either theoretical by finding the kernel and studying its decay rate, or numerical. In the numerical study, they took the dispersive term into account. As expected, they noticed that the “local dissipative term produces a bigger decay rate when compared with the nonlocal dissipative term”.

In [5], the authors proved the global existence of solution to the viscous KdV derived by [12] with the dispersive term and even, for sufficiently small initial conditions, without this dispersive term. In addition, they numerically investigated the decay rate for various norms.

In the present article, we first make a linear study of NSE in our domain (Section 2). We compute the dispersion relation and state various asymptotics that give different phase velocities, and so we give the decay rate in finite depth. In Section 3, we make the formal derivation of the viscous Boussinesq system by splitting the upper domain and the bottom one (the boundary layer). The explicit shear stress at the bottom is computed. In intermediate computations, there remains evaluations of the velocity at various heights which are expanded so as to replace these terms by the velocity at a generic height z without the assumption of irrotationality. This enables to state the viscous Boussinesq system. In Section 4, we state the 2D system, and cross-check with [12] that we get a similar viscous KdV equation. We also discuss the various viscous KdV equations proposed in the bibliography.

2. Linear theory

In order to make a linear theory, we need first to obtain dimensionless equations. This is done in the next subsection. Then we investigate two asymptotics in the following subsections.

Fig. 1.

The dimensionless domain.

2.1. Dimensionless equations

Let us denote $\tilde{u} = (\tilde{u}, \tilde{w})$ the velocity of a fluid in a 2D domain $\tilde{Ω} = {(\tilde{x}, \tilde{z}) / \tilde{x} \in R, \tilde{z} \in (- h, \tilde{η} (\tilde{x}, \tilde{t}))}$ . So we assume the bottom is flat and the free surface is characterized by $\tilde{z} = \tilde{η} (\tilde{x}, \tilde{t})$ with $\tilde{η} (\tilde{x}, \tilde{t}) > - h$ (the bottom does not get dry). Let $\tilde{p}$ denote the pressure and $\tilde{D} [\tilde{u}]$ the symmetric part of the velocity gradient. The dimensionless domain is drawn in Fig. 1. We also denote ρ the density of the fluid, ν the dynamic viscosity of the fluid, g the gravitational acceleration, k the unit vertical vector, n the outward unit normal to the upper frontier of $\tilde{Ω}$ , ${\tilde{p}}_{atm}$ the atmospheric pressure. The original system reads: $\{\begin{matrix} ρ (\frac{\partial \tilde{u}}{\partial \tilde{t}} + \tilde{u} \cdot \tilde{\nabla} \tilde{u}) - ν \tilde{Δ} \tilde{u} + \tilde{\nabla} \tilde{p} = - ρ g k & in \tilde{Ω}, \\ \tilde{div} \tilde{u} = 0 & in \tilde{Ω}, \\ (- \tilde{p} I + 2 ν \tilde{D} [\tilde{u}]) \cdot n = - {\tilde{p}}_{atm} n & on \tilde{z} = \tilde{η} (\tilde{x}, \tilde{t}), \\ {\tilde{η}}_{\tilde{t}} + \tilde{u} {\tilde{η}}_{\tilde{x}} - \tilde{w} = 0 & on \tilde{z} = \tilde{η} (\tilde{x}, \tilde{t}), \\ \tilde{u} = 0 & on \tilde{z} = - h, \end{matrix}$ (2) where we write the second order tensors and the vectors with bold letters. Of course, we need to add an initial condition and conditions at infinity.

So as to get dimensionless fields and variables, we need to choose a characteristic horizontal length l which is the wavelength (roughly the inverse of the wave vector), a characteristic vertical length h which is the water’s height, and the amplitude A of the propagating perturbation. Moreover, we denote U, W, P the characteristic horizontal velocity, vertical velocity and pressure, respectively. We may then define: $\begin{matrix} c_{0} = \sqrt{g h}, ε = \frac{A}{h}, β = \frac{h^{2}}{l^{2}}, U = ε c_{0}, \\ W = \frac{U l}{h} = \frac{c_{0} ε}{\sqrt{β}}, P = ρ g A, Re = \frac{ρ c_{0} h}{ν}, \end{matrix}$ (3) where $c_{0}$ is the phase velocity. As a consequence, one may make the fields dimensionless and unscaled: $\tilde{u} = U u, \tilde{w} = W w, \tilde{p} = {\tilde{p}}_{atm} - ρ g \tilde{z} + P p, \tilde{η} = A η,$ (4) and the variables: $\tilde{x} = l x, \tilde{z} = h (z - 1), \tilde{t} = t l / c_{0} .$ (5) With these definitions, the new system with the new fields and variables writes in the new domain $Ω_{t} = {(x, z), x \in R, z \in (0, 1 + ε η (x, t))}$ and the new outward unit normal still denoted n: $\{\begin{matrix} u_{t} + ε u u_{x} + \frac{ε}{β} w u_{z} - \frac{\sqrt{β}}{Re} u_{x x} - \frac{1}{Re \sqrt{β}} u_{z z} + p_{x} = 0 & in Ω_{t}, \\ w_{t} + ε u w_{x} + \frac{ε}{β} w w_{z} - \frac{\sqrt{β}}{Re} w_{x x} - \frac{1}{Re \sqrt{β}} w_{z z} + p_{z} = 0 & in Ω_{t}, \\ β u_{x} + w_{z} = 0 & in Ω_{t}, \\ (η - p) n + \frac{1}{Re} (\begin{matrix} 2 u_{x} \sqrt{β} & u_{z} + w_{x} \\ u_{z} + w_{x} & 2 / \sqrt{β} w_{z} \end{matrix}) \cdot n = 0 & on z = 1 + ε η, \\ η_{t} + ε u η_{x} - \frac{1}{β} w = 0 & on z = 1 + ε η, \\ u = 0 & on z = 0 . \end{matrix}$ (6) Like Kakutani and Matsuuchi [12], we could have eliminated $η - p$ in one of the two equations of stress continuity at the free boundary. After simplification by $1 / Re$ , this would have led us to the “simplified” system: $\{\begin{matrix} η - p + (- ε η_{x} (u_{z} + w_{x}) - 2 u_{x} \sqrt{β}) / Re = 0, \\ (1 - {(ε η_{x})}^{2}) (u_{z} + w_{x}) = 4 ε \sqrt{β} η_{x} u_{x} . \end{matrix}$

Notice that the number of dynamic conditions is linked to the Laplacian’s presence. If, in a subdomain, the flow is inviscid (Euler or $Re \to \infty$ ), then one must not keep the two above equations. Yet, once we have simplified the $1 / Re$ term above we might forget that the second equation must be taken off as if there remained a $1 / Re$ term before every term. So this “simplification” can be misleading.

Unlike us, the authors of [12] use the same characteristic length in the two space directions and so, for them, $h / l = 1$ . Our vertical velocity (scaled by W) is not the same as in [12]. Our choice of scale for W raises some $\sqrt{β}$ terms that [12] avoids. It suffices to set $β = 1$ in our equations to get those of [12]. Although the authors make their system dimensionless, they did not really unscale the fields nor the variables. Our fields are unscaled.

Our characteristic pressure is $P = ρ g A$ while [12] use $ρ g h$ . This explains why [12] has an ε more in front of the pressure p in their equations.

2.2. The dispersion relation

We are looking for small fields. So we linearize the system (6) and get: $\{\begin{matrix} u_{t} - \frac{\sqrt{β}}{Re} u_{x x} - \frac{1}{Re \sqrt{β}} u_{z z} + p_{x} = 0 & in R \times [0, 1], \\ w_{t} - \frac{\sqrt{β}}{Re} w_{x x} - \frac{1}{Re \sqrt{β}} w_{z z} + p_{z} = 0 & in R \times [0, 1], \\ β u_{x} + w_{z} = 0 & in R \times [0, 1], \\ η - p - \frac{2 \sqrt{β} u_{x}}{Re} = 0 & on z = 1, \\ u_{z} + w_{x} = 0 & on z = 1, \\ η_{t} - \frac{1}{β} w = 0 & on z = 1, \\ u = 0 & on z = 0 . \end{matrix}$ (7) First, we eliminate the pressure from (7)₁ and (7)₂: $u_{z t} - \frac{\sqrt{β}}{Re} u_{x x z} - \frac{1}{Re \sqrt{β}} u_{z z z} - w_{x t} + \frac{\sqrt{β}}{Re} w_{x x x} + \frac{1}{Re \sqrt{β}} w_{x z z} = 0 .$ Then we eliminate u from the previous equation thanks to (7)₃ by differentiating with respect to x and after some simplifications: $(\partial_{z}^{2} + β \partial_{x}^{2}) (\partial_{z}^{2} + β \partial_{x}^{2} - Re \sqrt{β} \partial_{t}) w = 0 .$ (8) If w is of the form $A (z) e^{i k (x - c t)}$ with a non-negative k and a (complex) phase velocity c, we can define a parameter with non-negative real part similar to the one used by [12]: $μ^{2} = β k^{2} - Re \sqrt{β} i k c .$ (9) Thanks to this notation, the solutions of (8) are such that $\begin{array}{rclr} A (z) & = & C_{1} cosh \sqrt{β} k (z - 1) + C_{2} sinh \sqrt{β} k (z - 1) \\ + C_{3} cosh μ (z - 1) + C_{4} sinh μ (z - 1) . & (10) \end{array}$

Up to now we have eliminated u and p only in the volumic equations. We still have to use the boundary conditions of (7) to find the conditions on the remaining field w.

The first equation of (7)₇ is $u (0) = 0$ . After a differentiation with respect to x and (7)₃, we get $w_{z} (0) = 0$ .

The second equation of (7)₇ is $w (0) = 0$ and needs no treatment.

Equation (7)₅ can be differentiated with respect to x and, thanks to (7)₃ leads to $w_{z z} - β w_{x x} = 0$ at height $z = 1$ .

The equation (7)₆ enables to compute/eliminate η.

The equation (7)₄ must be differentiated with respect to t for η to be replaced. Then we get $\frac{w}{β} - p_{t} + \frac{2}{Re \sqrt{β}} w_{z t} = 0 .$ We may differentiate the previous equation with respect to x so as to have a $p_{x}$ term which can be replaced thanks to (7)₁ to have new u terms. It suffices then to differentiate this equation and use the incompressibility (7)₃ to get the last condition. The full conditions on w are: $\begin{matrix} w_{z} (0) = 0, \\ w (0) = 0, \\ w_{z z} (1) - β w_{x x} (1) = 0, \\ w_{x x} (1) - w_{z t t} (1) + \frac{3 \sqrt{β}}{Re} w_{x x z t} (1) + \frac{1}{Re \sqrt{β}} w_{z z z t} (1) = 0 . \end{matrix}$ (11) The solutions (10) will satisfy a homogeneous linear system in the constants $C_{1}$ , $C_{2}$ , $C_{3}$ , $C_{4}$ . Its matrix is: $(\begin{matrix} \sqrt{β} k sinh (\sqrt{β} k) & - \sqrt{β} k cosh (\sqrt{β} k) & μ sinh μ & - μ cosh μ \\ cosh (\sqrt{β} k) & - sinh (\sqrt{β} k) & cosh μ & - sinh μ \\ 2 k^{2} β & 0 & μ^{2} + β k^{2} & 0 \\ - k^{2} & \sqrt{β} k^{3} c^{2} + \frac{2 i β k^{4} c}{Re} & - k^{2} & \frac{2 μ \sqrt{β} i k^{3} c}{Re} \end{matrix}) .$ (12) It suffices to compute its determinant to get the dispersion relation: $\begin{array}{rclr} 4 β k^{2} μ (β k^{2} + μ^{2}) + 4 μ k^{3} β^{3 / 2} (μ sinh (k \sqrt{β}) sinh μ - k \sqrt{β} cosh (k \sqrt{β}) cosh μ) \\ - {(β k^{2} + μ^{2})}^{2} (μ cosh (k \sqrt{β}) cosh μ - k \sqrt{β} sinh (k \sqrt{β}) sinh μ) \\ - k \sqrt{β} {Re}^{2} (μ sinh (k \sqrt{β}) cosh μ - k \sqrt{β} cosh (k \sqrt{β}) sinh μ) = 0 . & (13) \end{array}$ This relation is identical to the one of [12] except that our process of non-dimensionalizing makes a difference between x and z. So instead of k (for [12]), we have $k \sqrt{β}$ .

2.3. Asymptotic of the phase velocity (very large

Re

)

In this subsection, we prove the following proposition.

Proposition 1.
Under the assumptions $\begin{array}{rclr} k \sqrt{β} Re c \to + \infty, & (14) \\ k = O (1), & (15) \\ β \to 0, & (16) \\ Re \to + \infty, & (17) \\ c = O (1) (and c bounded away from 0), & (18) \end{array}$ if there exists a complex phase velocity c solution of ( 13 ), then it is such that: $c = \sqrt{\frac{tanh (k \sqrt{β})}{k \sqrt{β}}} - \frac{e^{i π / 4} {Re}^{- 1 / 2} {(k \sqrt{β})}^{1 / 4}}{2 {tanh}^{3 / 4} (k \sqrt{β})} + o (β^{- 1 / 4} {Re}^{- 1 / 2}) .$ (19) Moreover, the decay rate in our finite-depth geometry, which is viscous, is: $Im (k c) = \frac{- 1}{2 \sqrt{2}} \frac{k^{5 / 4} β^{1 / 8}}{\sqrt{Re} {tanh}^{3 / 4} (k \sqrt{β})} + o (β^{- 1 / 4} {Re}^{- 1 / 2}) .$ (20)

We denote $o (f)$ (respectively $O (f)$ ) a function which ratio with f tends to zero (respectively is bounded).

Our decay rate is not the same as Boussinesq’s or Lamb’s. The reason is that our geometry is not infinite. In the regime $Re = R ε^{- 5 / 2}$ and $β = b ε$ with constant R, b it gets: $Im (k c) = \frac{- \sqrt{k}}{2 \sqrt{2} \sqrt{Re \sqrt{β}}} + o (β^{- 1 / 4} {Re}^{- 1 / 2}) = \frac{- \sqrt{k} ε}{2 \sqrt{2} \sqrt{R \sqrt{b}}} + o (ε) .$ (21)

Our proposition is stated in [12] but not fully proved. One must also notice that the viscosity modifies both the real and the imaginary part of the phase velocity at the same order.
Proof of Proposition 1.
The definition of μ ( $ℜ (μ) ⩾ 0$ ) and assumptions (14), (15), (16), (18) enable to state that $μ^{2} \to \infty$ and the $k^{2} β$ term tends to zero. So we have: $μ = e^{- i π / 4} \sqrt{k \sqrt{β} Re c} + O (\frac{β^{3 / 4}}{Re}),$ (22) where the leading term tends to ∞ and its real part tends to $+ \infty$ , while the error term tends to zero. As a consequence, $tanh μ = 1 + O (e^{- μ})$ and $1 / cosh μ = O (e^{- μ})$ . Dividing (13) by $cosh μ$ and using a generic notation $P (β, μ)$ for an unspecified polynomial in β, μ, we have: $\begin{array}{rclr} O (P (β, μ) e^{- μ}) + 4 μ k^{4} β^{2} (μ \frac{sinh (k \sqrt{β})}{k \sqrt{β}} - cosh (k \sqrt{β})) \\ - {(k^{2} β + μ^{2})}^{2} (μ cosh (k \sqrt{β}) - k \sqrt{β} sinh (k \sqrt{β})) \\ - k^{2} β {Re}^{2} (μ \frac{sinh (k \sqrt{β})}{k \sqrt{β}} - cosh (k \sqrt{β})) = 0 . & (23) \end{array}$ The leading term of the second monomial is $4 k^{4} β^{2} μ^{2} sinh (k \sqrt{β}) / (k \sqrt{β})$ while the leading term of the fourth (last) is $- k^{2} β {Re}^{2} μ sinh (k \sqrt{β}) / (k \sqrt{β})$ . Seen the assumptions, their ratio is $4 k^{2} β μ {Re}^{- 2} = O (β^{5 / 4} {Re}^{- 3 / 2})$ . Under the assumptions (16), (17), this ratio tends to zero. So, in a first step, we can neglect the second monomial with respect to the fourth. If we look for a non-vanishing solution, we need to have a compensation of the only two remaining leading terms. One may then rewrite (23) as: $- (μ^{4} + hot) (μ cosh (k \sqrt{β}) + hot) - k^{2} β {Re}^{2} (μ \frac{sinh (k \sqrt{β})}{k \sqrt{β}} + hot) + hot = 0,$ where $hot$ denotes higher order terms. This reads after easy computations: $c^{2} = \frac{tanh (k \sqrt{β})}{k \sqrt{β}} + hot .$ (24) Such a relation is well known. It confirms the assumption (18). To pursue the expansion we come back to (23) and expand its various monomials starting with the second: $- 4 i k^{5} β^{5 / 2} Re c \frac{sinh (k \sqrt{β})}{k \sqrt{β}} + O (β^{3}) + O (β^{9 / 4} {Re}^{1 / 2}) .$ Indeed, even the leading term of this second monomial will be negligible in comparison with $O (β^{7 / 4} {Re}^{3 / 2})$ that we will have further. The third monomial of (23) is more complex and we must keep: $+ e^{- i π / 4} {(k \sqrt{β} Re c)}^{5 / 2} cosh (k \sqrt{β}) - k^{3} β^{3 / 2} {Re}^{2} c^{2} sinh (k \sqrt{β}) + O (β^{7 / 4} {Re}^{3 / 2}) .$ The fourth monomial of (23) is expanded: $- k^{2} β {Re}^{2} (e^{- i π / 4} \sqrt{k \sqrt{β} Re c} \frac{sinh (k \sqrt{β})}{k \sqrt{β}} - cosh (k \sqrt{β})) + O (β^{7 / 4} {Re}^{3 / 2}) .$ Using these expansions, Eq. (23) can be rewritten: $\begin{matrix} e^{- i π / 4} {(k \sqrt{β} Re)}^{5 / 2} \sqrt{c} cosh (k \sqrt{β}) [c^{2} - \frac{tanh (k \sqrt{β})}{k \sqrt{β}} + \frac{e^{i π / 4}}{\sqrt{k \sqrt{β} Re} \sqrt{c}}] \\ + O (P (β, μ) e^{- μ}) - k^{3} β^{3 / 2} {Re}^{2} c^{2} sinh (k \sqrt{β}) + O (β^{3}) + O (β^{7 / 4} {Re}^{3 / 2}) = 0 . \end{matrix}$ We would like to justify that the term between square brackets vanishes. For that purpose, we must check that the various other monomials are negligible in comparison with the third of those written between the square brackets which expands into: $O ({(\sqrt{β} Re)}^{5 / 2} [{(\sqrt{β} Re)}^{- 1 / 2}]) = O (β {Re}^{2})$ if we assume (18). Once it is checked (this is easy computation left to the reader), we can claim we proved that only the terms enclosed by square brackets remain: $c^{2} = \frac{tanh (k \sqrt{β})}{k \sqrt{β}} - \frac{e^{i π / 4}}{\sqrt{k \sqrt{β} Re c}} + o (β^{- 1 / 4} {Re}^{- 1 / 2}),$ (25) and the proof is complete by computing the square root of (25) and replacing the first order of c into (25) which leads to (19).

The computation of the decay rate is straightforward. □

We must stress that the complex phase velocity (19) contains two terms. The first is the classical gravitational term ( $\sqrt{tanh (k \sqrt{β}) / (k \sqrt{β})}$ ) which may be expanded when β tends to zero: $1 - k^{2} β / 6 + O (β^{2})$ . The second is purely viscous and can be expanded: $- \sqrt{2} (1 + i) {(4 \sqrt{k})}^{- 1} {(Re \sqrt{β})}^{- 1 / 2} + o {(Re \sqrt{β})}^{- 1 / 2}$ . So the dependences of c both on the gravitational and on the viscous effects are of the same order of magnitude when β and ${(Re \sqrt{β})}^{- 1 / 2}$ are of the same order. In this regime of very large Re, studied hereafter, the dependence of $Re$ on β is such that: $Re ≃ β^{- 5 / 2} .$ (26)
2.4. Second asymptotics of the phase velocity (moderate $Re$ )

The definition of $μ^{2}$ is $μ^{2} = k^{2} β - i k \sqrt{β} Re c$ and we still assume a long-wave asymptotics ( $β \to 0$ ). So one term or the other dominates in $μ^{2}$ . The extremes are either $μ^{2} \to \infty$ (see above) or $μ^{2} \to 0$ .

In the present subsection, we investigate the latter case and exhibit a more precise expansion than the one justified in [12]. Indeed, we prove the following proposition.

Proposition 2.
Under the assumptions $\begin{array}{rclr} k is bounded from zero and infinity, & (27) \\ μ \to 0 and so \sqrt{β} Re c \to 0, & (28) \\ β \to 0 (long waves), & (29) \\ Re \to + \infty, & (30) \end{array}$ if there exists a complex phase velocity c solution of ( 13 ), then it is such that: $c = - \frac{i k \sqrt{β} Re}{3} - \frac{19 i k^{3} β^{3 / 2} {Re}^{3}}{90} + o (β^{3 / 2} {Re}^{3}),$ (31) and necessarily ( 28 ) implies: $\sqrt{β} Re \to 0,$ (32) and so the phase velocity tends to zero.

Notice that if we assume $\sqrt{β} Re \to 0$ , the conclusion is the same and the proof much simpler.
Proof of Proposition 2.
Since $μ \to 0$ , we must expand all the functions in (13). In this expansion, we pay special attention to the fact that $Re \to + \infty$ and it may not be considered as a constant parameter of an expansion in β (hidden in $O (β^{2})$ as [12] did). After tedious expansions, there remains from (13): $\begin{array}{rclr} O (β {Re}^{2} c^{2} μ^{5}) + β {Re}^{2} μ^{7} (1 / 7! + o (μ)) + O (β^{3 / 2} Re c μ^{5}) + O (β^{2} μ^{5}) \\ + O (β^{3} {Re}^{2} c^{2} μ) + O (β^{7 / 2} Re c μ) + O (β^{4} {Re}^{2} μ) \\ + μ {Re}^{2} k^{2} β c [(c + \frac{i k \sqrt{β} Re}{3}) - \frac{i k \sqrt{β} Re c^{2}}{2} + \frac{4 k^{2} β {Re}^{2} c}{5} + 2 k^{2} β c \\ + \frac{8 i k^{3} β^{3 / 2} Re}{5} + \frac{i k^{3} β^{3 / 2} Re μ^{2}}{3 \times 5!}] = 0 . & (33) \end{array}$ Thanks to the assumptions (27)–(28) we can use that $\sqrt{β} Re c \to 0$ . Then, if we denote $T_{i}$ the ith term (among the eight) of this equation, and compare some of them, either to the first ( $μ {Re}^{2} k^{2} β c^{2}$ ) or to the second ( $- i μ {Re}^{3} k^{3} β^{3 / 2} c / 3$ ) of the terms inside the square brackets, we have: $\begin{array}{rcl} \frac{T_{1}}{μ {Re}^{2} k^{2} β c^{2}} = O (μ^{4}), \frac{T_{3}}{- i μ {Re}^{3} k^{3} β^{3 / 2} c / 3} = O (\frac{μ^{4}}{{Re}^{2}}), \\ \frac{T_{4}}{T_{2}} = O ({Re}^{- 2}), \frac{T_{5}}{μ {Re}^{2} k^{2} β c^{2}} = O (β^{2}), \\ \frac{T_{6}}{- i μ {Re}^{3} k^{3} β^{3 / 2} c / 3} = O (\frac{β^{2}}{{Re}^{2}}), \frac{T_{7}}{T_{2}} = O (\frac{β^{3}}{μ^{6}}) = O (1) . \end{array}$ As a consequence, the terms $T_{1}$ , $T_{3}$ , $T_{4}$ , $T_{5}$ , $T_{6}$ and $T_{7}$ can be taken off from (33). Then, if we simplify by $β {Re}^{2} μ$ and define D a constant, this equation writes: $\begin{array}{rclr} μ^{6} (D + O (μ)) + [c^{2} + \frac{i k (\sqrt{β} Re c)}{3} - c^{2} \frac{i k (\sqrt{β} Re c)}{2} + \frac{4 k^{2} {(\sqrt{β} Re c)}^{2}}{5} + 2 k^{2} β c^{2} \\ + \frac{8 i k^{3} β (\sqrt{β} Re c)}{5} + \frac{i k^{3} β μ^{2} (\sqrt{β} Re c)}{3 \times 5!}] = 0 . & (34) \end{array}$ Seen the assumptions (27)–(28), the highest order term in the square brackets is $c^{2}$ which must then vanish: $c \to 0$ . Moreover, one may see that $μ^{6} = O (β^{3}) + O {(\sqrt{β} Re c)}^{3}$ because of the complex definition of $μ^{2}$ . As a consequence, Eq. (34) simplifies first to $c = - \frac{i k \sqrt{β} Re}{3} + o (\sqrt{β} Re) .$ (35) Since we proved that $c \to 0$ , so does $\sqrt{β} Re$ as stated in (32). Moreover, by the definition of $μ^{2}$ and because of (35), one may claim $μ^{2} = - i k \sqrt{β} Re c (1 + O (1 / {Re}^{2})) \sim - k^{2} β {Re}^{2} / 3 .$ In a second step, before pursuing the expansion of c, one may notice that $μ^{6} = O (β^{3} {Re}^{6})$ which may then be neglected in (34). So, there remains only the terms in the square brackets (simplified by c): $\begin{array}{rclr} c + \frac{i k \sqrt{β} Re}{3} + \frac{i k^{3} β^{3 / 2} {Re}^{3}}{18} + o (β^{3 / 2} {Re}^{3}) - \frac{4 i k^{3} β^{3 / 2} {Re}^{3}}{15} + o (β^{3 / 2} {Re}^{3}) \\ - \frac{2 i k^{3} β^{3 / 2} Re}{3} + o (β^{3 / 2} Re) + \frac{8 i k^{3} β^{3 / 2} Re}{5} + \frac{i k^{3} β^{3 / 2} Re μ^{2}}{3 \times 5!} = O (β^{3}) + O (β^{3} {Re}^{6}) . & (36) \end{array}$ Among all the terms, one may justify that only the first to fourth must be kept: $c + \frac{i k \sqrt{β} Re}{3} + \frac{i k^{3} β^{3 / 2} {Re}^{3}}{18} - \frac{4 i k^{3} β^{3 / 2} {Re}^{3}}{15} = o (β^{3 / 2} {Re}^{3}),$ which gives the announced result (31) and completes the proof. □

Our phase velocity is different from Kakutani and Matsuuchi’s [12] because they assume constant Re (while it tends to infinity) and make expansions with the other parameter.
3. Formal derivation

We are going to consider the influence of viscosity on the solution of the Navier–Stokes equations in the domain $Ω_{t}$ . On the basis of the linear theory of the previous section, we assume a large Re and $Re ≃ β^{- 5 / 2}$ (37) as announced in (26). This is the case when viscous and gravitational effects balance in their influence on the variation of the phase velocity. We further make the Boussinesq assumption: $β \sim b ε,$ (38) where ε is already defined and b is a given positive number. Our main purpose here is to derive an asymptotic system of reduced size from the global Navier–Stokes equations in the whole moving domain. In the inviscid case, we would derive the classical Boussinesq system.

As will be justified below, these assumptions make that the NSE contain $u_{t} + ε u u_{x} + w u_{z} / b - ε^{2} u_{z z} / (R \sqrt{b}) + p_{x} = O (ε^{3}) .$ This enables to infer that in almost the whole domain, the second order derivative can be neglected but not in the boundary layer. Because of the square of ε and the second order derivative, this layer will precisely be of size ε vertically and close to the physical boundary.

In order to prove our main result, we proceed in the same way as [12] and distinguish two subdomains: the upper part ( $z > ε$ ) where viscosity can be neglected, and the lower part ( $0 < z < ε$ ) which is a boundary layer at the bottom and where viscosity must be taken into account. All the other geometrical characteristics have already been depicted. Our first main proposition is stated hereafter.

Proposition 3.

Let $η (x, t)$ be the free boundary’s height. Let $u^{b, 0} (x, γ)$ for $γ \in (0, + \infty)$ (resp. $u^{u, 0} (x, z)$ for $z \in (0, 1 + ε η (x, t))$ ) be the initial horizontal velocity in the boundary layer (resp. in the upper part of the domain). If $u^{b, 0} (x, γ)$ is uniformly continuous in γ and $u^{b, 0} (x, γ) - u^{u, 0} (x, z = 0)$ satisfies (58), then the solution of the Navier–Stokes equation with this given initial condition satisfies the Boussinesq system: $\begin{matrix} u_{t} + η_{x} + ε u u_{x} - β η_{x x x} \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ η_{t} + u_{x} (x, z, t) - \frac{β}{2} η_{x x t} (z^{2} - \frac{1}{3}) + ε {(u η)}_{x} + \frac{ε}{\sqrt{π R \sqrt{b}}} η_{t} * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \int_{γ^{'} = 0}^{\sqrt{(R \sqrt{b}) / (4 t)} γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}), \end{matrix}$ (39) where the Laplace convolution, denoted with ∗, is in time, the parameters ε, β, $Re$ have been defined and $z \in (0, 1 + ε η (x, t))$ .

If the initial velocity is a Euler flow, then $u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0) = 0$ (there is no viscous flow in the boundary layer) and the system writes: $\{\begin{matrix} u_{t} + η_{x} + ε u u_{x} - β η_{x x x} \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ η_{t} + u_{x} (x, z, t) - \frac{β}{2} η_{x x t} (z^{2} - \frac{1}{3}) + ε {(u η)}_{x} + \frac{ε}{\sqrt{π R \sqrt{b}}} η_{t} * \frac{1}{\sqrt{t}} = O (ε^{2}), \end{matrix}$ (40) where the Laplace convolution is still in time and defined by $\int_{0}^{t} η_{t} (x, t - t^{'}) / \sqrt{t^{'}} d t^{'}$ .

Before starting the proof, we must justify our non-obvious choice of method and a non-obvious term.

Remark 4.

Of course, the domain in the boundary layer $γ \in [0, + \infty [$ is not physical. Indeed, $u^{b}$ should be considered for γ between $γ = 0$ and γ large with respect to 1, but small with respect to $1 / ε$ so as to ensure $z = ε γ < 1$ . For instance, one may choose γ between 0 and $1 / \sqrt{ε}$ (equivalently $z = \sqrt{ε}$ ) or any value between $γ = 1$ and $γ = + \infty$ such that $z = ε γ ≪ 1$ .

The same applies in the upper part. Indeed, $u^{u} (x, z, t)$ should be considered for $z \in (ε, 1)$ and $u^{u} (x, z = 0, t)$ should be written $u^{u} (x, z = ε, t)$ .

The two domains overlap and we will write the matching conditions anywhere on the overlapping subdomain. For instance, one can then write the boundary condition at any height like $γ = 1 / \sqrt{ε}$ and force that the final result does not rely on this choice.

As is classical in boundary layer analysis, these more justified notations would give the same result as our choice. So we will use the simplest notation though less justified and consider $u^{u}$ for $z \in (0, 1 + ε η (x, t))$ and $u^{b}$ for $γ \in (0, + \infty)$ . Yet, the matching conditions will be for large γ (say $1 / \sqrt{ε}$ ).

Remark 5.

The double integral term in (39) is new and surprising because of its dependence on the initial condition. Assuming vanishing initial conditions in the boundary ( $u^{b, 0} = u^{u, 0}$ or equivalently that the initial flow is of Euler type), would greatly simplify the computations. But would it be physical? A physical question is then to know whether an initial (inviscid) flow in the boundary layer (where Navier–Stokes applies) establishes as a Navier–Stokes flow fast or not. The answer is negative for at least two reasons.

On the one hand, we claim that the characteristic time for the viscous effects to appear is roughly $T_{NSE} = ρ h_{0}^{2} / ν$ or $T_{NSE} = ρ l^{2} / ν$ . Then, its ratio with the characteristic time of the inviscid gravity flow ( $l / c_{0}$ ) is either $Re \sqrt{β} = ε^{- 2}$ or $Re / \sqrt{β} = ε^{- 3}$ , respectively. In any case, it is large and the flow in the boundary layer does not establish fast enough. It prohibits to claim that a Euler initial condition is physically compatible with Navier–Stokes equations for moderate times. Khabakhpashev [14] already discussed it although he started from rest!

On the second hand, let us assume the x differentiated discrepancy between the initial horizontal velocity in the boundary layer ( $u^{b, 0}$ ) and the trace at $z = 0$ of the initial horizontal velocity in the upper part ( $u^{u, 0}$ ) be an exponential. Namely $u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0) = G (x) e^{- α γ^{″}}$ . Such an expression is physical. In addition it leads to an explicit computation which provides for the whole term: $\frac{ε G (x) \sqrt{R \sqrt{b}}}{α^{2} \sqrt{π t}} + O (\frac{1}{t^{3 / 2}}) .$ As could have been guessed, the solution to the NSE tends to zero as $1 / \sqrt{t}$ . Such a slow convergence proves that one may not get rid of this term. Assuming the initial flow to be of Euler type while simulating NSE is incompatible.

In the first Section 3.1 we treat the upper part where convenient equations of (6) are kept. Then in Section 3.2, after a rescaling, we solve in the boundary layer the convenient equations extracted from (6). The solutions are forced to match through a continuity condition at the boundary ( $z = ε$ ) discussed in Section 3.3. At this stage, the system still has $u^{u} (1)$ and $\int_{0}^{1} u^{u}$ terms. So Section 3.4 is devoted to making explicit and simple the dependence on z so as to get rid of these extra terms.

3.1. Resolution in the upper part

The upper part is characterized by $ε < z < 1 + ε η (x, t)$ and $x, t \in R$ . We start from the system for the fields in the upper part and write u, w, p instead of $u^{u}$ , $w^{u}$ , $p^{u}$ for the sake of simplification. The height of the perturbation η is only defined in the upper part and so will always be denoted the same in the boundary layer. The system of PDE in the upper part is extracted from (6): $\{\begin{matrix} u_{t} + ε u u_{x} + \frac{ε}{β} w u_{z} - \frac{\sqrt{β}}{Re} u_{x x} - \frac{1}{Re \sqrt{β}} u_{z z} + p_{x} = 0, & ε < z < 1 + ε η, \\ w_{t} + ε u w_{x} + \frac{ε}{β} w w_{z} - \frac{\sqrt{β}}{Re} w_{x x} - \frac{1}{Re \sqrt{β}} w_{z z} + p_{z} = 0, & ε < z < 1 + ε η, \\ β u_{x} + w_{z} = 0, & ε < z < 1 + ε η, \\ - ε η_{x} (η - p) + \frac{1}{Re} (- 2 ε \sqrt{β} u_{x} η_{x} + (u_{z} + w_{x})) = 0 & on z = 1 + ε η, \\ η - p + \frac{1}{Re} (- ε η_{x} (u_{z} + w_{x}) - 2 \sqrt{β} u_{x}) = 0 & on z = 1 + ε η, \\ η_{t} + ε u η_{x} - \frac{1}{β} w = 0 & on z = 1 + ε η . \end{matrix}$ (41)

Since we assume $Re ≃ ε^{- 5 / 2}$ , the terms $\sqrt{β} / Re$ are of the order of $ε^{3}$ and the terms $1 / (Re \sqrt{β})$ of the order of $ε^{2}$ . This simplifies (41)₁ and (41)₂ and justifies to take off the Laplacian. As a consequence, we must not keep the two dynamic conditions (41)₄ and (41)₅ since they are associated to a Laplacian. We decide to drop (41)₄ and keep (41)₅.

Alternatively, one can stress that (41)₅ gives $η - p = O (ε^{3})$ and so the LHS of (41)₄ is $O (ε^{4}) + O (ε^{5 / 2})$ . Since we expand until the order two, one may claim the equation reduces to $0 = 0$ . But one could also simplify by $1 / Re$ ( $≃ ε^{5 / 2}$ ) and be driven to a new equation. This equation would provide one more condition to the two equations for two fields. It is not surprising to see that the final solution would then be $u = 0$ . The error is that we must drop one boundary condition. The above argument to get rid of (41)₄ is sufficient.

On this topic, the literature uses the same equations, but the argument for dropping one boundary condition is rarely explicited. In [12], Kakutani and Matsuuchi claim “the condition [(41)₄] is automatically satisfied” (p. 242 paragraph 3) which is either wrong (the equation disappears) or incomplete (what if they simplify by $Re = ε^{5 / 2}$ ?). In [10], Dutykh and Dias solve the same problem and write two equations (their (3) and (4)) among which they keep only one for the derivation without explaining this drop.

Let us come back to the resolution in the upper part. The equation (41)₃ gives w up to a constant that can be found in (41)₆: $w (x, z, t) = - β \int_{1 + ε η}^{z} u_{x} (x, z^{'}, t) d z^{'} + β (η_{t} + ε u (1 + ε η) η_{x}),$ (42) and we stress that this equation is exact. For the computations later, we need to expand this equation up to the third order: $w (x, z, t) = β (η_{t} + \int_{0}^{1} u_{x}) - β \int_{0}^{z} u_{x} + ε β {(u (1) η)}_{x} + O (ε^{3}) .$ (43) The second order of the previous equations suffices to determine p from (41)₂ up to a constant: $\begin{array}{rclr} p (x, z, t) & = & p (x, 1 + ε η, t) - β (η_{t t} + \int_{0}^{1} u_{x t}) (z - 1) \\ + β \int_{1}^{z} \int_{0}^{z^{'}} u_{x t} (x, z^{″}, t) d z^{″} d z^{'} + O (ε^{2}) . & (44) \end{array}$ Thanks to (41)₅ the constant may be found ( $p (1 + ε η) = η + O (ε^{3})$ ) and so: $\begin{array}{rclr} p (x, z, t) & = & η - β (η_{t t} + \int_{0}^{1} u_{x t}) (z - 1) + β \int_{1}^{z} \int_{0}^{z^{'}} u_{x t} (x, z^{″}, t) d z^{″} d z^{'} + O (ε^{2}) \\ = & η - β η_{t t} (z - 1) + β \int_{1}^{z} \int_{1}^{z^{'}} u_{x t} (x, z^{″}, t) d z^{″} d z^{'} + O (ε^{2}) . & (45) \end{array}$ Then the remaining field u satisfies (41)₁ at the first order: $u_{t} + η_{x} + ε u u_{x} + ε u_{z} (η_{t} + \int_{z}^{1} u_{x}) - β η_{x t t} (z - 1) - β η_{x x x} {(z - 1)}^{2} / 2 = O (ε^{2}),$ (46) where we have replaced the $u_{x x t}$ by $- η_{x x x}$ as usual.

We still have to solve the equations in the lower part.

3.2. Resolution in the boundary layer

We need first to recall some classical properties of Laplace transforms.

3.2.1. Some useful properties

Before solving the equations in the lower part, we list here some classical properties of the Laplace transform. We start from the definition $L (f) (p) = \hat{f} (p) = \int_{t \in R^{+}} f (t) e^{- p t} d t .$ (47) It is well known that the Laplace transform of the derivative is given by $L (f^{'}) (p) = - f (0) + p L (f) (p) .$ (48) If the two transforms $L (f) (p)$ and $L (g) (p)$ converge absolutely for $p = p_{0}$ , and if both f and g are absolutely integrable and bounded in every finite interval that does not include the origin such as $(p_{1}, p_{2})$ where $0 < p_{1} ⩽ p_{2}$ , then the Laplace transform of the convolution exists for p such that $ℜ (p) ⩾ ℜ (p_{0})$ ([8], Theorem 10.1), even converges absolutely, and satisfies: $L (f) (p) L (g) (p) = L (f * g) (p) .$ (49) Below, we use the following definition of the convolution, linked to the Laplace transform: $f_{1} * f_{2} (t) = \int_{0}^{t} f_{1} (u) f_{2} (t - u) d u .$ (50) These formulas will be useful in the next subsection.

3.2.2. The fields in the boundary layer

The lower part of the domain ( $0 < z < ε$ ) is a boundary layer. We start from the system for the bottom fields, written u, w, p instead of $u^{b}$ , $w^{b}$ , $p^{b}$ for the sake of simplification and extracted from (6): $\{\begin{matrix} u_{t} + ε u u_{x} + \frac{ε}{β} w u_{z} - \frac{\sqrt{β}}{Re} u_{x x} - \frac{1}{Re \sqrt{β}} u_{z z} + p_{x} = 0 & for 0 < z < ε, \\ w_{t} + ε u w_{x} + \frac{ε}{β} w w_{z} - \frac{\sqrt{β}}{Re} w_{x x} - \frac{1}{Re \sqrt{β}} w_{z z} + p_{z} = 0 & for 0 < z < ε, \\ β u_{x} + w_{z} = 0 & for 0 < z < ε, \\ u (z = 0) = 0 and w (z = 0) = 0 . \end{matrix}$ (51)

As is justified in Section 2.3, the viscous and gravitational effects balance when $Re ≃ β^{- 5 / 2}$ (same as (26)). So we remind the reader of the assumptions $Re = R ε^{- 5 / 2}$ and $β = b ε$ for constant R, b. We are naturally led to change the scale in z as in any boundary layer. Let us introduce a new vertical variable $γ = z / ε$ . The new fields should be denoted in another way. Nevertheless, we will not change the notation for the sake of simplification. The new system writes: $\{\begin{matrix} u_{t} + ε u u_{x} + \frac{1}{b ε} w u_{γ} - \frac{\sqrt{b}}{R} ε^{3} u_{x x} - \frac{u_{γ γ}}{R \sqrt{b}} + p_{x} = 0, \\ w_{t} + ε u w_{x} + \frac{1}{b ε} w w_{γ} - \frac{\sqrt{b}}{R} ε^{3} w_{x x} - \frac{w_{γ γ}}{R \sqrt{b}} + \frac{p_{γ}}{ε} = 0, \\ ε β u_{x} + w_{γ} = 0, \\ u (γ = 0) = 0 and w (γ = 0) = 0 . \end{matrix}$ (52) One must notice that the Laplacian lets some remaining terms of zeroth degree (in ε) in this system. So the viscosity is relevant inside the boundary layer.

We can find the vertical velocity from (52)₃ and (52)₄: $w (x, γ, t) = - ε β \int_{0}^{γ} u_{x} (x, γ^{'}, t) d γ^{'} .$ (53) Carrying backward the previous equation in (52)₂, one has $p_{γ} = O (ε^{3})$ . So as to determine p, we need to use the continuity relation for the pressure ( $p (x, γ = 1, t) = p^{u} (x, z = ε, t)$ ) unless we cannot go on. Since we know the pressure in the upper part $p^{u}$ from (45), we can write: $p (x, γ, t) = p (x, γ = 1, t) + O (ε^{3}) = p^{u} (x, ε, t) + O (ε^{3}) = η (x, t) + O (ε) .$ (54) Using this equation and (53) in (52)₁, we have at zeroth order: $u_{t} + η_{x} - \frac{u_{γ γ}}{R \sqrt{b}} = O (ε) .$ (55) This equation must be completed with initial condition $u (x, γ, t = 0) = u^{b, 0} (x, γ)$ (56) and boundary condition: $\{\begin{matrix} u (x, γ = 0, t) = 0, \\ u (x, γ \to + \infty, t) = u^{u} (x, z = 0, t) (continuity condition) . \end{matrix}$ (57) These are the equations to be solved.

Since we solve a Cauchy problem for a heat-like equation, we have an initial condition and so we must use the time-Laplace transform. In [12], the authors do not take an initial condition, and uses a time-Fourier transform. In all his articles, Liu and coauthors (e.g. [18]), quote [23] (pp. 153–159) in which a sine-transform (in γ) is used, but the initial condition is set to zero. In a separate calculation, not reproduced here, we used the same sine-transform in γ and paid attention to the initial condition. We were led to the very same result as the one stated hereafter.

We solve the system (55)–(57) in the following lemma.

Lemma 6.

If the initial conditions $u^{b, 0} (x, γ)$ and $u^{u, 0} (x, z = 0)$ are uniformly continuous in γ and satisfy $\begin{matrix} \int_{0}^{\infty} | u^{b, 0} (x, γ) - u^{u, 0} (x, z = 0) | d γ < \infty, \\ \int_{0}^{\infty} | u_{x}^{b, 0} (x, γ) - u_{x}^{u, 0} (x, z = 0) | d γ < \infty, \end{matrix}$ (58) then the solution to ( 55 )–( 57 ) is $\begin{array}{rclr} u (x, γ, t) & = & u^{u} (x, z = 0, t) + \frac{\sqrt{R \sqrt{b}}}{2} \int_{0}^{+ \infty} f_{0} (x, γ^{'}) \frac{e^{- (R \sqrt{b} {(γ^{'} - γ)}^{2}) / (4 t)}}{\sqrt{π t}} d γ^{'} \\ - u^{u} (x, 0, \cdot) * L^{- 1} (e^{- σ γ}) \\ - \frac{\sqrt{R \sqrt{b}}}{2} \int_{0}^{+ \infty} f_{0} (x, γ^{'}) \frac{e^{(- R \sqrt{b} {(γ^{'} + γ)}^{2}) / (4 t)}}{\sqrt{π t}} d γ^{'} + O (ε), & (59) \end{array}$ where $f_{0} (x, γ) = u^{b, 0} (x, γ) - u^{u, 0} (x, z = 0)$ , $u^{u}$ is the horizontal velocity in the upper part that satisfies ( 46 ) and σ is the only root with a positive real part of $R \sqrt{b} p$ : $σ = σ (p) = \sqrt{R \sqrt{b} p},$ (60) where p is the dual variable of time t and the convolution is in time.

Remark 7.

The solution of (55) may be known only up to any function of x. The boundary condition (57) enables to determine this function.

Remark 8.

The compatibility of the conditions (56) and (57) forces to have, when γ tends to $+ \infty$ : $u^{b, 0} (x, γ) \to u^{u, 0} (x, z = 0),$ and, when $γ \to 0$ : $u^{b, 0} (x, γ = 0) = 0 .$

Meanwhile we also prove the following proposition.

Proposition 9.

Under the same assumptions as in Lemma 6, the bottom shear stress is $\begin{array}{rclr} τ^{b} & = & {(\frac{\partial u^{b}}{\partial γ})}_{γ = 0} = \frac{\sqrt{R \sqrt{b}} u^{u} (x, z = 0, 0)}{\sqrt{π}} p.v. \frac{1}{\sqrt{t}} \\ + \frac{\sqrt{R \sqrt{b}}}{\sqrt{π}} \int_{0}^{t} \frac{u_{t}^{u} (x, z = 0, t - s)}{\sqrt{s}} d s, & (61) \end{array}$ where p.v. denotes the principal value as defined in the theory of distributions.

First let us prove Proposition 9.

Proof of Proposition 9.

The initial condition $f_{0}$ may not make any difference (it can be seen through an explicit computation), so the corresponding term is taken off. Then a simple differentiation with respect to γ and the following formula (see [8], p. 320) $L^{- 1} (e^{- a \sqrt{p}}) = \frac{a}{2 \sqrt{π} t^{3 / 2}} e^{- a^{2} / (4 t)},$ applied to (59) for any positive a leads to $\begin{array}{rcl} τ^{b} & = & - \frac{d}{d γ} (\int_{0}^{t} u^{u} (x, z = 0, t - s) \frac{e^{- (R \sqrt{b} γ^{2}) / (4 s)} \sqrt{R \sqrt{b}} γ}{2 \sqrt{π} s^{3 / 2}} d s) + O (ε) \\ = & - \sqrt{R \sqrt{b}} \int_{0}^{t} \frac{u^{u} (x, z = 0, t - s)}{2 \sqrt{π} s^{3 / 2}} e^{- (R \sqrt{b} γ^{2}) / (4 s)} d s \\ - \frac{\sqrt{R \sqrt{b}}}{\sqrt{π}} \int_{0}^{t} \frac{u^{u} (x, z = 0, t - s)}{s^{1 / 2}} (- \frac{R \sqrt{b} γ^{2}}{4 s^{2}} e^{- (R \sqrt{b} γ^{2}) / (4 s)}) d s + O (ε) . \end{array}$ The second term may be integrated by parts to get $\begin{matrix} - \frac{\sqrt{R \sqrt{b}}}{\sqrt{π}} (\frac{u^{u} (x, z = 0, 0)}{\sqrt{t}} e^{- (R \sqrt{b} γ^{2}) / (4 t)} \\ - \int_{0}^{t} (- \frac{u_{t}^{u} (x, z = 0, t - s)}{\sqrt{s}} - \frac{u^{u} (x, z = 0, t - s)}{2 s^{3 / 2}}) e^{- (R \sqrt{b} γ^{2}) / (4 s)} d s), \end{matrix}$ which simplifies partially with the first term. At the end, there remains only $\frac{\sqrt{R \sqrt{b}}}{\sqrt{π}} \frac{u^{u} (x, z = 0, 0)}{\sqrt{t}} e^{- (R \sqrt{b} γ^{2}) / (4 t)} + \frac{\sqrt{R \sqrt{b}}}{\sqrt{π}} \int_{0}^{t} \frac{u_{t}^{u} (x, z = 0, t - s)}{\sqrt{s}} e^{- (R \sqrt{b} γ^{2}) / (4 s)} d s .$ This justifies the formula as is classical in the theory of distributions. □

The scheme of the proof of Lemma 6 is to solve (55) up to two unknown functions, then to determine these functions so as to satisfy the initial and boundary conditions. This provides a necessary formula. We check in the Appendix that the solution satisfies the boundary and initial conditions. Let us prove Lemma 6.

Proof of Lemma 6.

Let us denote $f (x, γ, t) = u (x, γ, t) - u^{u} (x, z = 0, t) .$ (62) Since $f_{t} = u_{t} + η_{x} + O (ε)$ (thanks to (46)) and $f_{γ} = u_{γ}$ , Eq. (55) writes: $f_{t} - f_{γ γ} / (R \sqrt{b}) = O (ε) .$ (63) The initial condition is $f (x, γ, t = 0) = u^{b, 0} (x, γ) - u^{u, 0} (x, z = 0) = : f_{0} (x, γ)$ (64) and the boundary conditions read $\begin{matrix} f (x, γ = 0, t) = - u^{u} (x, 0, t), \\ lim_{γ \to + \infty} lim_{ε \to 0} u (x, γ, t) - u^{u} (x, z = ε, t) = lim_{γ \to + \infty} f (x, γ, t) = 0 . \end{matrix}$ (65) The second condition is merely the continuity condition of the horizontal velocity at the border of the boundary layer. So we are driven to a heat equation in a half space with vanishing condition at infinity, and non-homogeneous initial and bottom conditions. Through a Laplace transform in time, denoted either $L (f)$ or $\hat{f}$ , (63) becomes $- f_{0} (x, γ) + p \hat{f} (p) - \frac{{\hat{f}}_{γ γ}}{R \sqrt{b}} = O (ε) .$ (66) In order to solve this non-homogeneous ODE, we start with the homogeneous one and recall that we define σ as the only root with a positive real part of $R \sqrt{b} p$ in (60). Its solutions are $\hat{f} (x, γ, p) = C_{1} (x, p) e^{+ σ γ} + C_{2} (x, p) e^{- σ γ} + O (ε) .$ By applying the method of parameters variation, we look for $C_{1} (x, γ, p)$ , $C_{2} (x, γ, p)$ such that: $- C_{1, γ} σ e^{σ γ} + C_{2, γ} σ e^{- σ γ} = R \sqrt{b} f_{0} (x, γ) + O (ε),$ and solving (66) amounts to solving the system of two equations with two unknown functions $C_{1}$ and $C_{2}$ : $\{\begin{matrix} C_{1, γ} e^{σ γ} + C_{2, γ} e^{- σ γ} = 0, \\ - C_{1, γ} e^{σ γ} + C_{2, γ} e^{- σ γ} = \frac{R \sqrt{b}}{σ} f_{0}, \end{matrix}$ which solution is (thanks to assumption (58)): $\{\begin{matrix} C_{1} (x, γ, p) = - \frac{R \sqrt{b}}{2 σ} \int_{+ \infty}^{γ} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'} + {\tilde{C}}_{1} (x, p), \\ C_{2} (x, γ, p) = + \frac{R \sqrt{b}}{2 σ} \int_{0}^{γ} f_{0} (x, γ^{'}) e^{σ γ^{'}} d γ^{'} + {\tilde{C}}_{2} (x, p) . \end{matrix}$ So, the full solution is $\begin{array}{rcl} \hat{f} (x, γ, p) & = & - \frac{R \sqrt{b}}{2 σ} \int_{+ \infty}^{γ} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'} e^{+ σ γ} + {\tilde{C}}_{1} (x, p) e^{+ σ γ} \\ + \frac{R \sqrt{b}}{2 σ} \int_{0}^{γ} f_{0} (x, γ^{'}) e^{σ γ^{'}} d γ^{'} e^{- σ γ} + {\tilde{C}}_{2} (x, p) e^{- σ γ} + O (ε) . \end{array}$ We look for ${\tilde{C}}_{1}$ first. Since $f_{0}$ is bounded, simple bounds prove that the first, third and fourth terms are bounded. So ${\tilde{C}}_{1} (x, p) = 0 .$ The unknown function ${\tilde{C}}_{2} (x, p)$ is then given by the boundary condition (65)₁ at the bottom: ${\tilde{C}}_{2} (x, p) = - u^{u} (x, z = 0, p) - \frac{R \sqrt{b}}{2 σ} \int_{0}^{+ \infty} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'} .$ In a necessary way, $\begin{array}{rclr} \hat{f} (x, γ, p) & = & + \frac{R \sqrt{b}}{2 σ} \int_{γ}^{+ \infty} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'} e^{+ σ γ} \\ + \frac{R \sqrt{b}}{2 σ} \int_{0}^{γ} f_{0} (x, γ^{'}) e^{σ γ^{'}} d γ^{'} e^{- σ γ} \\ - ({\hat{u}}^{u} (x, z = 0, p) + \frac{R \sqrt{b}}{2 σ} \int_{0}^{+ \infty} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'}) e^{- σ γ} + O (ε) . & (67) \end{array}$ From the definition of f, the existence of an inverse Laplace transform and formula (49), one knows that: $\begin{array}{rcl} f (x, γ, t) & = & \frac{R \sqrt{b}}{2} \int_{γ}^{+ \infty} f_{0} (x, γ^{'}) L^{- 1} (\frac{e^{- σ (γ^{'} - γ)}}{σ}) d γ^{'} \\ + \frac{R \sqrt{b}}{2} \int_{0}^{γ} f_{0} (x, γ^{'}) L^{- 1} (\frac{e^{σ (γ^{'} - γ)}}{σ}) d γ^{'} \\ - u^{u} (x, z = 0, \cdot) * L^{- 1} (e^{- σ γ}) \\ - \frac{R \sqrt{b}}{2} \int_{0}^{+ \infty} f_{0} (x, γ^{'}) L^{- 1} (\frac{e^{- σ (γ^{'} + γ)}}{σ}) d γ^{'} + O (ε) . \end{array}$ Owing to the formula (see [8]) $L^{- 1} (\frac{e^{- \tilde{a} \sqrt{p}}}{\sqrt{p}}) = \frac{1}{\sqrt{π t}} e^{- {\tilde{a}}^{2} / (4 t)},$ if $\tilde{a} > 0$ , one may justify the explicit form of u given in (59). Until the end of this article, we denote the function of time t: $A = A (t) = \sqrt{\frac{R \sqrt{b}}{4 t}} .$ (68) We still have to check that the initial condition (64) and remaining of the boundary conditions (65)₂ are satisfied by u given by (59). This is completed in the Appendix.

So we completed the proof of the whole Lemma 6. □

From (53) and (59), we can then compute the vertical velocity $\begin{array}{rclr} w^{b} (x, γ, t) & = & - ε β \int_{0}^{γ} u_{x}^{b} (x, γ^{'}, t) d γ^{'} \\ = & - ε β u_{x}^{u} (x, 0, t) γ + ε β u_{x}^{u} (x, 0, \cdot) * L^{- 1} (\frac{e^{- σ γ} - 1}{- σ}) \\ - ε β \frac{A (t)}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} - γ^{'})}^{2}} d γ^{″} d γ^{'} \\ + ε β \frac{A (t)}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} + γ^{'})}^{2}} d γ^{″} d γ^{'} + O (ε^{2} β) . & (69) \end{array}$ We still have to satisfy the continuity conditions of all the fields u, w, p.

3.3. Continuity conditions

In the present subsection, we need to write explicitly the superscripts u and b for the upper part and bottom regions, respectively. We write the computed fields at the same height ε that is the common frontier of both subdomains.

We already used the continuity of pressure that led us to (54). So the pressure is continuous.

Regarding the horizontal velocity, we must notice that the limit when $γ \to + \infty$ of ${lim}_{ε \to 0} (u^{b} (x, γ, t) - u^{u} (x, ε γ, t)) = f (x, γ, t)$ has already been computed as vanishing (see the Appendix). So the boundary condition (65)₂ is already satisfied and the horizontal velocity is continuous.

Concerning the vertical velocity, we can use the velocity in the upper part $w^{u}$ from (43) expanded in ε: $\begin{array}{rcl} w^{u} (x, ε γ, t) & = & β (η_{t} + \int_{0}^{1} u_{x}^{u}) - β \int_{0}^{ε γ} u_{x}^{u} + ε β {(u^{u} (1) η)}_{x} + O (ε^{3}) \\ = & β (η_{t} + \int_{0}^{1} u_{x}^{u}) - β ε γ u_{x}^{u} (z = 0) + ε β {(u^{u} (1) η)}_{x} + O (ε^{3}) . \end{array}$ One may notice that as anywhere else, the $u^{u} (z = 0)$ could be replaced by $u^{u} (z = ε)$ and $\int_{0}^{1} u_{x}^{u}$ by $\int_{ε}^{1 + ε η} u_{x}^{u}$ and so on. The formula would be the same and the final result would be the same.

The velocity in the bottom $w^{b}$ is given in (69). The difference $w^{u} (x, ε γ, t) - w^{b} (x, γ, t)$ can be expanded in ε: $\begin{array}{rclr} w^{u} - w^{b} & = & β (η_{t} + \int_{0}^{1} u_{x}^{u}) - β ε γ u_{x}^{u} (z = 0) + ε β {(u^{u} (1) η)}_{x} + O (ε^{3}) \\ + ε β u_{x}^{u} (x, z = 0, t) γ - ε β u_{x}^{u} (x, 0, \cdot) * L^{- 1} (\frac{e^{- σ γ} - 1}{- σ}) \\ + ε β \frac{A}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} - γ^{'})}^{2}} d γ^{″} d γ^{'} \\ - ε β \frac{A}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} + γ^{'})}^{2}} d γ^{″} d γ^{'} + O (ε^{3}) \\ = & β (η_{t} + \int_{0}^{1} u_{x}^{u}) + ε β {(u^{u} (1) η)}_{x} - ε β u_{x}^{u} (x, 0, \cdot) * L^{- 1} (\frac{1}{σ}) \\ + ε β \frac{A}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} - γ^{'})}^{2}} d γ^{″} d γ^{'} \\ - ε β \frac{A}{\sqrt{π}} \int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) e^{- A^{2} {(γ^{″} + γ^{'})}^{2}} d γ^{″} d γ^{'} + O (ε^{3}), & (70) \end{array}$ up to functions that tend exponentially to zero when $γ \to + \infty$ .

We still must simplify the two last double integrals. This is made in the following lemma.

Lemma 10.
If $A = A (t) = \sqrt{\frac{R \sqrt{b}}{4 t}}$ , γ is positive, $f_{0} (x, γ)$ is uniformly continuous in γ and satisfies ( 58 ), then $\int_{γ^{'} = 0}^{γ} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) (e^{- A^{2} {(γ^{″} - γ^{'})}^{2}} - e^{- A^{2} {(γ^{″} + γ^{'})}^{2}}) d γ^{″} d γ^{'}$ tends to $\int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) \int_{γ^{‴} = - γ^{″}}^{γ^{″}} e^{- A^{2} γ^{''' 2}} d γ^{‴} d γ^{″},$ (71) when $γ \to + \infty$ .

The proof relies on Fubini’s theorem and changes of variables for the two integrals. The proof is only technical and left to the reader.

After simplification by β, the continuity of the vertical velocity (70) reads after making $γ \to + \infty$ thanks to Lemma 10: $\begin{array}{rclr} η_{t} + \int_{0}^{1} u_{x}^{u} + ε {(u^{u} (1) η)}_{x} - \frac{ε}{\sqrt{π R \sqrt{b}}} u_{x}^{u} (x, 0, t) * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} f_{0, x} (x, γ^{″}) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}), & (72) \end{array}$ where the convolution is in time t and the formula $L^{- 1} (\frac{1}{\sqrt{p}}) = 1 / \sqrt{π t}$ [8] is used. If one had made the more rigorous expansion according to Remark 4, assuming $u^{u}$ is defined only on $(ε, 1 + ε η)$ and $u^{b}$ is defined on $(0, 1)$ , one would have been led to $\begin{array}{rclr} u_{t} + \int_{ε}^{1 + ε η} u_{x}^{u} + ε u^{u} (1 + ε η) η_{x} - \frac{ε}{\sqrt{π R \sqrt{b}}} u_{x}^{u} (x, z = ε, t) * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{1 / \sqrt{ε}} f_{0, x} (x, γ^{″}) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}) . & (73) \end{array}$
3.4. The dependence on z of the fields

At this stage, we have reduced the equations but not as much as in the Euler case which leads to a Boussinesq system in $1 + 1$ dimension. We have derived only a $2 + 1$ dimension problem although we have eliminated w and p. The major difference with the Boussinesq derivation comes from the assumption of irrotationality of Euler flows. This assumption would provide $u_{z} = O (ε)$ . Such a condition would annihilate the dependence on z and greatly simplify the above computations.

Yet irrotationality and its corollary of a potential flow is incompatible with the number of conditions we set at the bottom, which are needed by the dissipativity of the Navier–Stokes equations. So we need to determine the dependence on z of u to have a more tractable system.

Starting from now, we drop the u superscripts for the fields in the upper part but keep the superscripts for the boundary layer. In summary, we assume $Re ≃ ε^{- 5 / 2}$ , and the assumptions of the first asymptotic stated in the Section 2.3. Up to now, the reduced equations are collected from (46) and (72): $\begin{array}{rclr} u_{t} + η_{x} + ε u u_{x} + ε u_{z} (η_{t} + \int_{z}^{1} u_{x}) - β η_{x t t} (z - 1) - β η_{x x x} {(z - 1)}^{2} / 2 = O (ε^{2}) \forall z, & (74) \\ η_{t} + \int_{0}^{1} u_{x} (z) d z + ε {(u (z = 1) η)}_{x} - \frac{ε}{\sqrt{π R \sqrt{b}}} u_{x} (x, z = 0, t) * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}) . & (75) \end{array}$ Equation (74) can be rewritten thanks to the order 0 of (75): $u_{t} + η_{x} + ε u u_{x} - ε u_{z} \int_{0}^{z} u_{x} - β η_{x t t} (z - 1) - β η_{x x x} {(z - 1)}^{2} / 2 = O (ε^{2}) \forall z .$ (76) Notice that the $η_{x x x}$ term comes from an integral of the shape $\int_{1}^{z} \int_{1}^{z^{'}} u_{x x t}$ . As an intermediate result one may see very easily that $η_{x x} = η_{t t} + O (ε)$ which is useful later.

We intend to prove the following lemma.

Lemma 11.

A localized solution of ( 75 ) , ( 76 ) is such that $\begin{array}{rclr} \int_{0}^{1} u = u (x, z, t) - β η_{x t} \frac{z^{2} - 1 / 3}{2} + O (ε^{2}), & (77) \\ u (x, 0, t) = u (x, z, t) - β η_{x t} \frac{z^{2}}{2} + O (ε^{2}), & (78) \\ u (x, 1, t) = u (x, z, t) + β η_{x t} \frac{1 - z^{2}}{2} + O (ε^{2}), & (79) \\ \frac{1}{1 + ε η - ε_{1}} \int_{ε_{1}}^{1 + ε η} u = u (x, z, t) - \frac{β}{2} η_{x t} (z^{2} - 1 / 3) + O (ε^{2}), & (80) \end{array}$ as soon as $ε_{1} = O (ε)$ .

Of course, $z = 0$ is not in the domain of u ( $= u^{u}$ ) but it is meaningful as a trace.

Proof of Lemma 11.

In a preliminary step, we prove $u_{z} (x, z, t) = β η_{x t} (x, t) z + O (ε^{2}) .$ (81) To that end, we differentiate (76) with respect to z, so as to have: $u_{z t} + ε u u_{x z} - ε u_{z z} \int_{0}^{z} u_{x} - β η_{x t t} - β η_{x x x} (z - 1) = O (ε^{2}),$ and we can integrate this equation in time using that $η_{x x} = η_{t t} + O (ε)$ : $u_{z} + ε \int_{t_{0}}^{t} (u u_{x z}) - ε \int_{t_{0}}^{t} (u_{z z} \int_{0}^{z} u_{x}) - β η_{x t} - β η_{x t} (z - 1) = C_{3} (x, z) + O (ε^{2}),$ (82) where $C_{3}$ is a function of x, z but it does not depend on t. Since the solution is localized for any x, z, there exists a time $t_{0}$ at which $u_{z}$ and $η_{x t}$ vanish or are as small as wanted (in a local norm). So $C_{3} (x, z) = O (ε),$ in a first attempt to determine $C_{3}$ . But then Eq. (82) implies $u_{z} = O (ε)$ and so the quadratic terms are all of second order in (82) since they contain at least one $u_{z}$ . Hence $u_{z} (x, z, t) - β η_{x t} z = C_{4} (x, z) + O (ε^{2}) .$ Again since for all $(x, z)$ there exists a time at which u and η vanish or are as small as wanted, then $C_{4} (x, z) = O (ε^{2})$ and this completes the proof of (81). We can then go further by integrating between $z^{'}$ and z: $u (x, z, t) = u (x, z^{'}, t) + β η_{x t} \frac{z^{2} - z^{' 2}}{2} + O (ε^{2}),$ and then, integrating in $z^{'}$ between $z^{'} = 0$ and $z^{'} = 1$ , we can state (77). By integrating between $z^{'} = ε_{1}$ and $z^{'} = 1 + ε η$ , one gets (80). Setting $z^{'} = 0$ (or $z^{'} = ε$ ), we obtain (78) and setting $z^{'} = 1$ we obtain (79). □

So the systems (75), (76) can be rewritten thanks to (77)–(79), the formula $L^{- 1} (\frac{1}{\sqrt{p}}) = 1 / \sqrt{π t}$ [8], and the fact that, as in the Euler case the wave equation is the zeroth order ( $η_{x x} = η_{t t} + O (ε)$ ): $\begin{array}{rclr} u_{t} + η_{x} + ε u u_{x} - β η_{x x x} \frac{(z^{2} - 1)}{2} = O (ε^{2}), & (83) \\ η_{t} + u_{x} (x, z, t) - \frac{β}{2} η_{x x t} (z^{2} - \frac{1}{3}) + ε {(u η)}_{x} - \frac{ε}{\sqrt{π R \sqrt{b}}} u_{x} * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}), & (84) \end{array}$ where all the fields u are evaluated at $(x, z, t)$ and the convolution is in time. This is the system stated in Proposition 3 and the proof is complete.

Alternatively, one could have stated the system after an averaging of u ( $= u^{u}$ defined only in the upper part) from $z = ε_{1}$ to $z = 1 + ε η (x, t)$ . Using (80) and defining $U (x, t) = \frac{1}{1 + ε η (x, t) - ε_{1}} \int_{ε_{1}}^{1 + ε η (x, t)} u (x, z, t) d z,$ it is a simple exercise to justify $\begin{array}{rclr} U_{t} + η_{x} + ε U U_{x} + β η_{x x x} / 3 = O (ε^{2}), & (85) \\ η_{t} + U_{x} + ε {(U η)}_{x} + \frac{ε}{\sqrt{π R \sqrt{b}}} η_{t} * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}) . & (86) \end{array}$ A careful examination of (86) may raise questions. Indeed the incompressibility equation ( $β u_{x}^{u} + w_{z}^{u} = 0$ ) can be integrated vertically between $z = ε_{1}$ (close to the bottom of the upper part and to ε) and $z = 1 + ε η$ (free boundary). This provides exactly: $β (η_{t} + U_{x} + ε {(U η)}_{x}) = w^{u} (z = ε_{1}) + ε_{1} β U_{x} .$ Usually, the no-flow condition makes that $w (ε_{1}) = 0$ and the height at the bottom is zero. So the RHS vanishes. The discrepancy with (86) stems from the fact that $ε_{1} \neq 0$ and that some particles may cross the (virtual) bottom of the upper part. One may then retrieve one more information by matching the previous equation and (86). This provides the vertical velocity at $ε_{1} = O (ε)$ : $\begin{array}{rcl} w^{u} (ε_{1}) & = & ε_{1} β η_{t} - \frac{ε β}{\sqrt{π R \sqrt{b}}} η_{t} * \frac{1}{\sqrt{t}} \\ - \frac{2 ε β}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} + O (ε^{2}) . \end{array}$ One could find the same result from (69), assuming $ε_{1} = ε γ$ (large γ), using $w^{b} (ε γ) = w^{u} (ε_{1})$ and Lemma 10.

In [6], starting from the NSE for a flow down an inclined plane, two models are derived and compared. The first is Shkadov’s. It is an integration along the normal to the plane variable. The second uses generalized characteristics and a conservative form. None of them contain a dispersive third order derivative. Looking in an hyperbolic way to these two models, the authors conclude that “the more complex system is needed in the case of a shock formation” interpreted as a nonlinear wave-breaking. None of these system look like ours above since their regime is $ε < {Re}^{- 1} < \sqrt{ε}$ with our notation (see their (3)). Moreover, they even stress that if “ ${Re}^{- 1} < [ε]$ , then the […] scaling excludes […] terms with the second order derivatives” (p. 52 but we modified the quote according to our definition of ε which is the square of [6]’s). Our regime is precisely such that ${Re}^{- 1} \sim ε^{5 / 2} < ε$ and we made a detailed analysis of the boundary layer.

4. Generalization and checkings

In a first subsection, we state the 2D Boussinesq system and check we may find the classical Boussinesq systems in the inviscid case. Then, in Section 4.2 we derive rigorously the viscous KdV equation and discuss its compatibility with the equation derived by Kakutani and Matsuuchi in [12], by Liu and Orfila in [18] and by Dutykh in [9].

4.1. The full 2D Boussinesq systems family

One may start from the 3D Navier–Stokes equations and derive in a way very similar to above a generalization of (83), (84): $\{\begin{matrix} u_{t} + η_{x} + ε u u_{x} + ε v u_{y} - β (η_{x x x} + η_{x y y}) \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ v_{t} + η_{y} + ε u v_{x} + ε v v_{y} - β (η_{y x x} + η_{y y y}) \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ η_{t} + u_{x} + v_{y} - \frac{β}{2} (η_{x x t} + η_{y y t}) (z^{2} - \frac{1}{3}) + ε {(u η)}_{x} + ε {(v η)}_{y} + \frac{ε}{\sqrt{π R \sqrt{b}}} η_{t} * \frac{1}{\sqrt{t}} \\ + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} div (\begin{matrix} u^{b, 0} (x, γ^{″}) - u^{u, 0} (x, z = 0) \\ v^{b, 0} (x, γ^{″}) - v^{u, 0} (x, z = 0) \end{matrix}) \int_{γ^{'} = 0}^{A (t) γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} = O (ε^{2}) . \end{matrix}$ (87)

In case of a Euler initial condition, the last integral term vanishes, but this is not physical as is stressed in Remark 5.

One may recover the classical Boussinesq system by noticing that when $R \to + \infty$ , the convolution term vanishes. Moreover, the resolution of (55) disappears since the highest order term vanishes. This enables to get rid of the double integral term. One has then a classical Boussinesq system.

It is well known thanks to [1] that there is a family of Boussinesq systems, indexed by three free parameters. All these systems are equivalent in the sense that up to order 1, they can be derived one from the other by using their own $O (ε^{0})$ part and by replacing partially $η_{t}$ , $η_{x}$ and $η_{y}$ by $u_{x}$ , $u_{t}$ , $v_{t}$ . We are going to prove the same for our system. Namely, the order 0 of (87) enables to interpolate with $a_{int}$ , $b_{int}$ , $c_{int}$ : $\{\begin{matrix} η_{x} = a_{int} η_{x} - (1 - a_{int}) u_{t} + O (ε), \\ η_{y} = b_{int} η_{y} - (1 - b_{int}) v_{t} + O (ε), \\ η_{t} = c_{int} η_{t} - (1 - c_{int}) (u_{x} + v_{y}) + O (ε) . \end{matrix}$ These formulas are reported in the full 2D system (87), where we drop the convolution term and the integral on the initial condition: $\{\begin{matrix} u_{t} + η_{x} + ε u u_{x} + ε v u_{y} - a_{int} β Δ η_{x} \frac{(z^{2} - 1)}{2} \\ + (1 - a_{int}) β Δ u_{t} \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ v_{t} + η_{y} + ε u v_{x} + ε v v_{y} - b_{int} β Δ η_{y} \frac{(z^{2} - 1)}{2} \\ + (1 - b_{int}) β Δ v_{t} \frac{(z^{2} - 1)}{2} = O (ε^{2}), \\ η_{t} + u_{x} + v_{y} - c_{int} \frac{β}{2} Δ η_{t} (z^{2} - \frac{1}{3}) \\ + (1 - c_{int}) \frac{β}{2} Δ (u_{x} + v_{y}) (z^{2} - \frac{1}{3}) + ε {(u η)}_{x} + ε {(v η)}_{y} = O (ε^{2}), \end{matrix}$ (88) where we denote Δ the x, y Laplacian.

This is the general Boussinesq system as can be seen in [1] (p. 285, Eq. (1.6)). Indeed if we denote $a_{BCS}$ , $b_{BCS}$ , $c_{BCS}$ and $d_{BCS}$ the interpolation parameters of this article, we can identify the 1D version of our interpolated (88) with $\begin{matrix} a_{BCS} = \frac{β}{2} (1 - c_{int}) (z^{2} - \frac{1}{3}), b_{BCS} = \frac{β}{2} c_{int} (z^{2} - \frac{1}{3}), \\ c_{BCS} = - β a_{int} \frac{z^{2} - 1}{2}, d_{BCS} = - β (1 - a_{int}) \frac{z^{2} - 1}{2} . \end{matrix}$ The meaning of our height z is the same as the θ of [1] and the relation between $a_{BCS}$ , $b_{BCS}$ , $c_{BCS}$ and $d_{BCS}$ (see (1.8) of this article) is satisfied.

4.2. About the KdV-like equation

Various authors have derived either a viscous Boussinesq system or a viscous KdV equation.

One may wonder what is the viscous KdV equation derived from our viscous Boussinesq system and compare it with what may be found in the literature. First, we state and prove the following proposition.

Proposition 12.
If the initial flow is localized, the KdV change of variables applied to the system ( 83 ) , ( 84 ) leads to $2 {\tilde{η}}_{τ} + 3 \tilde{η} {\tilde{η}}_{ξ} + \frac{b}{3} {\tilde{η}}_{ξ ξ ξ} - \frac{1}{\sqrt{π R \sqrt{b}}} \int_{ξ^{'} = 0}^{τ / ε} \frac{{\tilde{η}}_{ξ} (ξ + ξ^{'}, τ)}{\sqrt{ξ^{'}}} d ξ^{'} = O (ε)$ (89) for not too small times τ, where we set $Re = R ε^{- 5 / 2}$ and $β = b ε$ .

In the formula (89), since it has been proved in [16] that KdV is a good approximation of Euler for times up to $1 / ε^{2}$ , and that the velocity is localized, it is a strong temptation to replace the integral term by $- \frac{1}{\sqrt{π R \sqrt{b}}} \int_{ξ^{'} = 0}^{+ \infty} \frac{{\tilde{η}}_{ξ} (ξ + ξ^{'}, τ)}{\sqrt{ξ^{'}}} d ξ^{'} .$

This is the term found in [12].
Proof of Proposition 12.
We start from the most general form of (83), (84) (same as (39)) and use the KdV change of variables $(ξ = x - t, τ = ε t) \Leftrightarrow (x = ξ + τ / ε, t = τ / ε),$ (90) and change of fields $Φ (x, z, t) = \tilde{Φ} (x - t, z, ε t) \Rightarrow Φ_{t} = - {\tilde{Φ}}_{ξ} + ε {\tilde{Φ}}_{τ} (x - t, z, ε t),$ (91) where the generic field Φ is tilded when it depends on the $(ξ, z, τ)$ variables.

There are only two difficult terms in the system (83), (84) (equivalent to (39)). The first is the convolution which we denote $T_{1}$ : $\begin{array}{rcl} T_{1} (x, z, t) & = & - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{t^{'} = 0}^{t} \frac{u_{x} (x, z, t - t^{'})}{\sqrt{t^{'}}} d t^{'} \\ = & - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{t^{'} = 0}^{t} \frac{{\tilde{u}}_{ξ} (x - t + t^{'}, z, ε t - ε t^{'})}{\sqrt{t^{'}}} d t^{'}, \end{array}$ because of (90). But then it suffices to recognize the function of $(x - t, ε t) = (ξ, τ)$ in the last equation to have the term after the KdV change of variables: $\begin{matrix} {\tilde{T}}_{1} (ξ, z, τ) & = - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{t^{'} = 0}^{τ / ε} \frac{{\tilde{u}}_{ξ} (ξ + t^{'}, z, τ - ε t^{'})}{\sqrt{t^{'}}} d t^{'} \\ = - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{ξ^{'} = 0}^{τ / ε} \frac{{\tilde{u}}_{ξ} (ξ + ξ^{'}, z, τ)}{\sqrt{ξ^{'}}} d ξ^{'} + O (ε^{2}) . \end{matrix}$ (92) Since the $t^{'}$ variable is in place of a ξ, we changed the notation to $ξ^{'}$ . This term is odd because it has an integration variable ( $ξ^{'}$ ) that has a physical meaning but bounds depending on time $τ / ε$ . We discuss it below.

The second difficult term is the one that keeps the initial conditions and writes: $T_{2} (x, z, t) = + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (x, γ^{″}) - u_{x}^{u, 0} (x, z = 0)) \times \int_{γ^{'} = 0}^{\sqrt{(R \sqrt{b}) / (4 t)} γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} .$ The change of variables (90) gives: $\begin{array}{rcl} {\tilde{T}}_{2} (ξ, z, τ) & = & + \frac{2 ε}{\sqrt{π}} \int_{γ^{″} = 0}^{+ \infty} (u_{x}^{b, 0} (ξ + \frac{τ}{ε}, γ^{″}) - u_{x}^{u, 0} (ξ + \frac{τ}{ε}, z = 0)) \\ \times \int_{γ^{'} = 0}^{\sqrt{(R \sqrt{b} ε) / (4 τ)} γ^{″}} e^{- γ^{' 2}} d γ^{'} d γ^{″} . \end{array}$ If the initial boundary layer is localized, for τ not too small, $u_{x}^{b, 0} (ξ + \frac{τ}{ε}, γ^{″}) - u_{x}^{u, 0} (ξ + \frac{τ}{ε}, z = 0)$ will be small in $L_{γ^{″}}^{1}$ and so ${\tilde{T}}_{2}$ will be negligible in comparison with ε and so can be dropped. In addition, the inner integral’s upper bound is very close to the lower bound.

Then, we can claim that the Boussinesq system after the KdV change of variables and fields is $\{\begin{matrix} - {\tilde{u}}_{ξ} + ε {\tilde{u}}_{τ} + {\tilde{η}}_{ξ} + ε \tilde{u} {\tilde{u}}_{ξ} - β {\tilde{η}}_{ξ ξ ξ} (\frac{z^{2} - 1}{2}) = O (ε^{2}), \\ - {\tilde{η}}_{ξ} + ε {\tilde{η}}_{τ} + {\tilde{u}}_{ξ} + \frac{β}{2} {\tilde{η}}_{ξ ξ ξ} (z^{2} - \frac{1}{3}) + ε {(\tilde{u} \tilde{η})}_{ξ} \\ - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{ξ^{'} = 0}^{τ / ε} \frac{{\tilde{u}}_{ξ} (ξ + ξ^{'}, z, τ)}{\sqrt{ξ^{'}}} d ξ^{'} = O (ε^{2}) . \end{matrix}$ (93) We may notice that at the first order, and as in the derivation of the KdV equation, ${\tilde{u}}_{ξ} = {\tilde{η}}_{ξ} + O (ε) \Rightarrow \tilde{u} = \tilde{η} + O (ε),$ thanks to a simple and classical integration (and a localized solution). But then the sum of the two equations of (93) gives: $ε {\tilde{u}}_{τ} + ε {\tilde{η}}_{τ} + ε \tilde{u} {\tilde{u}}_{ξ} + ε {(\tilde{u} \tilde{η})}_{ξ} + \frac{β}{3} {\tilde{η}}_{ξ ξ ξ} - \frac{ε}{\sqrt{π R \sqrt{b}}} \int_{ξ^{'} = 0}^{τ / ε} \frac{{\tilde{u}}_{ξ} (ξ + ξ^{'}, z, τ)}{\sqrt{ξ^{'}}} d ξ^{'} = O (ε^{2}) .$ Using now the fact that $\tilde{u} = \tilde{η} + O (ε)$ , dividing by ε, one states exactly Eq. (89). The convolution that used to be on time is now on $ξ^{'}$ and the proof is complete. □

What can be found in the literature?

As stated in the Introduction, various authors already derived either a viscous Boussinesq system or a viscous KdV equation. Yet, none of them have the very same equation as us. We must clarify why there are such differences.

The first article is [24] in which Ott and Sudan obtained formally in nine lines: $+ α_{3} \int_{ξ^{'} = - \infty}^{+ \infty} \frac{{\tilde{u}}_{ξ} (ξ^{'}, τ) sgn (ξ - ξ^{'})}{\sqrt{| ξ - ξ^{'} |}} d ξ^{'},$ but the authors used a Fourier transform [12] (p. 243) and they made an error pointed by [12]. Our formula differs from Ott and Sudan’s by the sign and the bound!

Later, Kakutani and Matsuuchi [12] derived rather rigorously the KdV equation from Navier–Stokes and we set the same regime as them. Yet, they did not raise the problem of the initial condition. As a consequence, they used a time-Fourier transform to solve the heat-like equation. They proposed: $\frac{1}{4 \sqrt{π R}} \int_{ξ^{'} = - \infty}^{+ \infty} \frac{{\tilde{η}}_{ξ} (ξ^{'}, τ) (1 - sgn (ξ - ξ^{'}))}{\sqrt{| ξ - ξ^{'} |}} d ξ^{'} .$ Their half order derivative term differs from ours only by the bound of the integral which is $τ / ε$ for us and $+ \infty$ for them. This is not relevant since this term is multiplied by an ε. So one may say that their term is identical to ours.

Liu and Orfila in [18] (and subsequent articles) derived a Boussinesq system for a regime different from ours ( $Re = R ε^{- 7 / 2}$ ). They solved their heat equation with a sine-transform in the vertical coordinate by quoting [23] where a vanishing initial condition is assumed. Given their regime, their Boussinesq system is right. But when they derived a KdV equation (see [18], p. 89), they did not make explicit their change of variables in the term equivalent to our $T_{1} (x, z, t)$ . With the change of variable $ξ_{LO} = x - t$ , $τ_{LO} = (α_{LO} / μ_{LO}) t$ , they exhibit (see their (3.19) or (3.21)): $- \frac{1}{2 \sqrt{π}} \int_{0}^{t} \frac{η_{ξ_{LO}}}{\sqrt{t - T}} d T,$ where there remains the former variable t inside the integral and in the bounds. Moreover, the dependence of $η_{ξ_{LO}}$ on the variables ( $t, τ_{LO}, \dots$ ?) is not written. Is the T variable in the integral a time variable? One may wonder whether they did notice that the time convolution transforms into a space one.

Dutykh derived a Boussinesq system by a Helmholtz–Leray decomposition from a Linearized Navier–Stokes [9]. In order to derive the associated KdV (see his Section 3.2), he assumed $u = η + ε P + β Q + \dots$ and found P and Q. In this process, he used only the assumption that waves go right ( $η_{t} + η_{x} = O (ε)$ ). So he did not use the change of time ( $τ = ε t$ ) and wrote a formula with unscaled time t (his (14)): $- \sqrt{\frac{ν}{π} \frac{g}{h}} \int_{0}^{t} \frac{η_{x}}{\sqrt{t - τ}} d τ .$ Similar criticisms to Liu’s can be said on this formula in which the integral seems to be on time while it should be on the shifted space ξ.
5. Conclusion

In this article, we derive the viscous Boussinesq system for surface waves from Navier–Stokes equations with non-vanishing initial conditions (see Proposition 3). We also provide the Boussinesq system with vertically-integrated horizontal velocity (85)–(86). We justify (Remark 5) that assuming a Euler-type initial velocity and viscous evolution equations is incompatible. One of our by-product is the bottom shear stress as a function of the velocity (cf. Proposition 9) and the decay rate for shallow water (see Proposition 1). We also state the system in 3D in (87), and derive the viscous KdV equation from our viscous Boussinesq system (cf. Proposition 12). The differences of our viscous KdV with other equations, already derived in the literature, are highlighted and explained.

Footnotes

Acknowledgements

The author wants to thank Professor Jean-Claude Saut for initiating and following this research and Professor David Gérard-Varet for a fruitful discussion. Lastly, the author wants to thank a referee who made valuable comments that improved the initial manuscript.

Boundary and initial conditions in Lemma 6

As is said in the proof of Lemma 6, we must check that u, given by the necessary Eq. (59), satisfies the initial condition (64) and the remaining of the boundary conditions (65)₂.

Concerning the initial condition (64). We try to find the limit when t tends to $0^{+}$ and then $A = A (t)$ tends to $+ \infty$ . Since one assumes below $γ > 0$ , $- u^{u} (x, 0, \cdot) * L^{- 1} (e^{- σ γ}) = - u^{u} (x, 0, t) * \frac{\sqrt{R \sqrt{b}}}{2 \sqrt{π} t^{3 / 2}} e^{- (R \sqrt{b}) / (4 t)}$ tends to zero exponentially (the convolution is the Laplace one and on time t). Then, one can come back to the formula of f (67) and make one change of variables in every integral: $\begin{array}{rcl} f (x, γ, t) & = & \frac{A (t)}{\sqrt{π}} \int_{- γ}^{+ \infty} f_{0} (x, Γ^{'} + γ) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} \\ - \frac{A (t)}{\sqrt{π}} \int_{γ}^{+ \infty} f_{0} (x, Γ^{'} - γ) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} + O (ε), \end{array}$ up to an exponentially tending to zero function when t tends to 0 thanks to $A (t)$ . This can be rewritten $\begin{array}{rcl} f (x, γ, t) & = & \frac{A (t)}{\sqrt{π}} \int_{γ}^{+ \infty} (f_{0} (x, Γ^{'} + γ) - f_{0} (x, Γ^{'} - γ)) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} \\ + \frac{A (t)}{\sqrt{π}} \int_{- γ}^{γ} f_{0} (x, Γ^{'} + γ) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} + O (ε), \end{array}$ where we denote $I_{2}$ the second integral. The first integral may be bounded by $\begin{array}{rcl} \frac{2 A (t)}{\sqrt{π}} sup_{γ > 0} | f_{0} (x, γ) | \int_{γ}^{+ \infty} e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} \\ ⩽ \frac{2}{π} sup_{γ > 0} | f_{0} (x, γ) | \int_{A (t) γ}^{+ \infty} e^{- Γ^{'' 2}} d Γ^{″}, \end{array}$ which clearly tends to zero when t tends to zero ( $A = A (t) \to + \infty$ ).

For the second integral denoted $I_{2}$ , one may compute a similar integral where the integration variable of $f_{0}$ is frozen: $\begin{matrix} I_{2}^{'} & = \frac{A (t)}{\sqrt{π}} \int_{- γ}^{γ} f_{0} (x, γ) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} \\ = f_{0} (x, γ) \frac{1}{\sqrt{π}} \int_{- A (t) γ}^{A (t) γ} e^{- Γ^{'' 2}} d Γ^{″}, \end{matrix}$ which clearly tends to $f_{0} (x, γ)$ if $γ > 0$ when $t \to 0^{+}$ . So one may make the difference of the second integral $I_{2}$ with the previous integral (which tends to $f_{0} (x, γ)$ ) and find: $I_{2} - I_{2}^{'} = \frac{A (t)}{\sqrt{π}} \int_{- γ}^{γ} (f_{0} (x, Γ^{'} + γ) - f_{0} (x, γ)) e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} + o_{t \to 0^{+}} (1) .$ Here we must use the assumption of uniform continuity of the initial data: $\forall ϵ > 0 \exists δ > 0 / | γ^{'} - γ | < δ \Rightarrow | f_{0} (x, γ^{'}) - f_{0} (x, γ) | < ϵ .$ Then, for any $ϵ > 0$ , there exists a δ such that $I_{2} - I_{2}^{'}$ can be splitted into two parts and bounded by $\begin{matrix} \frac{2 A (t)}{\sqrt{π}} sup_{γ > 0} | f_{0} (x, γ) | \int_{| Γ^{'} | > δ \cap | Γ^{'} | < γ} e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} + \frac{A (t)}{\sqrt{π}} ϵ \int_{| Γ^{'} | < δ \cap | Γ^{'} | < γ} e^{- A {(t)}^{2} Γ^{' 2}} d Γ^{'} \\ ⩽ \frac{2 A (t)}{\sqrt{π}} sup_{γ > 0} | f_{0} (x, γ) | 2 γ e^{- A {(t)}^{2} δ^{2}} + ϵ . \end{matrix}$ So we have proved that the f given by (67) or u given by (59) satisfies the initial condition ( $A (t) \to + \infty$ when $t \to 0^{+}$ ).

Concerning the boundary condition (65)₂. Now we look for the limit when γ tends to $+ \infty$ . The formula (67) can be written: $\hat{f} (x, γ, p) = \frac{R \sqrt{b}}{2 σ} \int_{γ}^{+ \infty} f_{0} (x, γ^{'}) e^{- σ γ^{'}} d γ^{'} e^{+ σ γ} + \frac{R \sqrt{b}}{2 σ} \int_{0}^{γ} f_{0} (x, γ^{'}) e^{σ γ^{'}} d γ^{'} e^{- σ γ},$ up to some exponentially tending to zero functions of γ. In this formula, the first integral is bounded by $\begin{matrix} \frac{R \sqrt{b}}{2 σ} sup_{γ^{'} ⩾ γ} | f_{0} (x, γ^{'}) | \int_{γ}^{+ \infty} e^{- σ γ^{'}} d γ^{'} e^{σ γ} \\ ⩽ \frac{R \sqrt{b}}{2 σ^{2}} sup_{γ^{'} ⩾ γ} | f_{0} (x, γ^{'}) |, \end{matrix}$ which clearly tends to zero when γ tends to $+ \infty$ because $f_{0} (x, γ)$ tends to zero when γ tends to $+ \infty$ .

For the second integral, one needs to cut it at a value Γ given by the definition of $f_{0} \to 0$ when γ tends to $+ \infty$ ( $\forall ϵ > 0$ $\exists Γ > 0 / | γ | > Γ \Rightarrow | f_{0} | < ϵ$ ). We can bound it with: $\frac{R \sqrt{b}}{2 σ} \int_{0}^{Γ} | f_{0} (x, γ^{'}) | e^{σ γ^{'}} d γ^{'} e^{- σ γ} + \frac{R \sqrt{b}}{2 σ} ϵ \int_{Γ}^{γ} e^{σ γ^{'}} d γ^{'} e^{- σ γ} .$ Since the first term tends to zero when γ tends to $+ \infty$ (Γ fixed) and the second term is less than $R \sqrt{b} ϵ / (2 σ^{2})$ , the whole can be made smaller than any ϵ.

So the proof that (65)₂ is satisfied is complete.

References

[1]

J.L.

Bona,

Chen and

J.-C.

Saut, Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I. Derivation and linear theory, J. Nonlinear Sci. 12(4) (2002), 283–318.

[2]

Boussinesq, Lois de l’extinction d’une houle simple en haute mer, C. R. Math. Acad. Sci. Paris 121(1) (1895), 15–19.

[3]

J.G.B.

Byatt-Smith, The effect of laminar viscosity on the solution of the undular bore, J. Fluid Mech. 48(1) (1971), 33–40.

[4]

Chen,

Dumont,

Dupaigne and

Goubet, Decay of solutions to a water wave model with a nonlocal viscous dispersive term, Discrete Contin. Dyn. Syst. 27(4) (2010), 1473–1492.

[5]

Chen,

Dumont and

Goubet, Decay of solutions to a viscous asymptotical model for waterwaves: Kakutani–Matsuuchi model, Nonlinear Anal. 75(5) (2012), 2883–2896.

[6]

A.A.

Chesnokov and

P.V.

Kovtunenko, Weak discontinuities in solutions of long-wave equations for viscous flow, Studies in Applied Mathematics 132(1) (2014), 50–64.

[7]

Dias,

A.I.

Dyachenko and

V.E.

Zakharov, Theory of weakly damped free-surface flows: A new formulation based on potential flow solutions, Physics Letters A 372 (2008), 1297–1302.

[8]

Doetsch, Introduction to the Theory and Application of the Laplace Transformation, Springer-Verlag, New York, 1974. (Translated from the second German edition by W. Nader.)

[9]

Dutykh, Visco-potential free-surface flows and long wave modelling, Eur. J. Mech. B Fluids 28(3) (2009), 430–443.

10.

[10]

Dutykh and

Dias, Viscous potential free-surface flows in a fluid layer of finite depth, C. R. Math. Acad. Sci. Paris 345(2) (2007), 113–118.

11.

[11]

R.S.

Johnson, A Modern Introduction to the Mathematical Theory of Water Waves, Cambridge Univ. Press, Cambridge, 1997.

12.

[12]

Kakutani and

Matsuuchi, Effect of viscosity of long gravity waves, J. Phys. Soc. Japan 39(1) (1975), 237–246.

13.

[13]

G.H.

Keulegan, Gradual damping of solitary waves, J. Res. Natl. Bur. Stand. 40 (1948), 607–614.

14.

[14]

G.A.

Khabakhpashev, Effect of bottom friction on the dynamics of gravity perturbations, Fluid Dynamics 22(3) (1987), 430–437.

15.

[15]

Lamb, Hydrodynamics, 3rd edn, Cambridge Univ. Press, Cambridge, 1906.

16.

[16]

Lannes, The Water Waves Problem: Mathematical Analysis and Asymptotics, Vol. 188, AMS, 2013.

17.

[17]

P.L.-F.

Liu and

I.-C.

Chan, On long-wave propagation over a fluid-mud seabed, J. Fluid Mech. 579 (2007), 467–480.

18.

[18]

P.L.-F.

Liu and

Orfila, Viscous effects on transient long-wave propagation, J. Fluid Mech. 520 (2004), 83–92.

19.

[19]

P.L.-F.

Liu,

Y.S.

Park and

E.A.

Cowen, Boundary layer flow and bed shear stress under a solitary wave, J. Fluid Mech. 574 (2007), 449–463.

20.

[20]

P.L.-F.

Liu,

Simarro,

Vandever and

Orfila, Experimental and numerical investigations of viscous effects on solitary wave propagation in a wave tank, Coastal Engineering 53 (2006), 181–190.

21.

[21]

Matsuuchi, Numerical investigations on long gravity waves under the influence of viscosity, J. Phys. Soc. Japan 41(2) (1976), 681–687.

22.

[22]

C.C.

Mei, The Applied Dynamics of Ocean Surface Waves, John Wiley & Sons, Inc., New York, 1983.

23.

[23]

C.C.

Mei, Mathematical Analysis in Engineering: How to Use the Basic Tools, Cambridge Univ. Press, Cambridge, 1995.

24.

[24]

Ott and

R.N.

Sudan, Damping of solitary waves, Phys. Fluids 13(6) (1970), 1432–1434.

25.

[25]

Y.S.

Park,

P.L.-F.

Liu and

S.J.

Clark, Viscous flows in a muddy seabed induced by a solitary wave, J. Fluid Mech. 598 (2008), 383–392.