Near-critical reflection of internal waves

TL;DR

This paper investigates the reflection behavior of internal waves at critical angles on sloping boundaries, demonstrating that solutions can be accurately approximated by incident, reflected, and boundary layer components.

Contribution

It provides a rigorous mathematical proof that solutions exist and are well approximated in critical reflection geometry, confirming prior predictions.

Findings

Solution exists and is well approximated by wave components

Reflected second harmonic and boundary layers are significant

Validates previous theoretical predictions

Abstract

Internal waves describe the (linear) response of an incompressible stably stratified fluid to small perturbations. The inclination of their group velocity with respect to the vertical is completely determined by their frequency. Therefore the reflection on a sloping boundary cannot follow Descartes' laws, and it is expected to be singular if the slope has the same inclination as the group velocity. In this paper, we prove that in this critical geometry the weakly viscous and weakly nonlinear wave equations have actually a solution which is well approximated by the sum of the incident wave packet, a reflected second harmonic and some boundary layer terms. This result confirms the prediction by Dauxois and Young, and provides precise estimates on the time of validity of this approximation.

Tables

Table 1. Table 1: Sizes of the boundary layers in different regimes

$\left|\zeta \right|\gtrsim 1$	${\nu }^{\frac{1}{4}}\lesssim \left|\zeta \right|\ll 1$	${\nu }^{\frac{1}{3}}\ll \left|\zeta \right|\ll {\nu }^{\frac{1}{4}}$	$\left|\zeta \right|\lesssim {\nu }^{\frac{1}{3}}$
• One reflected wave • One BL of size ${\nu }^{1/2}$ (lifting 2 conditions)	• One highly oscillating reflected wave (with slow decay) • One BL of size ${\left(\nu /\zeta \right)}^{1/2}$ • One BL of size ${\nu }^{1/2}$	• One BL of size ${\left|\zeta \right|}^4/\nu$ • One BL of size ${\left(\nu /\zeta \right)}^{1/2}$ • One BL of size ${\nu }^{1/2}$	• One BL of size ${\nu }^{1/3}$ • One BL of size ${\nu }^{1/2}$ (lifting 2 conditions)

Table 2. Table 2: List of quadratic interactions

	type of interaction	size in $L^2$	typical decay rate
(a1)	$Q({𝒲}_{BL,{\epsilon }^2}^0,{𝒲}_{BL,{\epsilon }^2}^0)$	$O\left(1\right)$	${\epsilon }^{-2}$
(a2)	$Q({𝒲}_{\text{inc}}^0,{𝒲}_{BL,{\epsilon }^2}^0)$	$O\left(1\right)$	${\epsilon }^{-2}$
(b1)	$Q({𝒲}_{BL,{\epsilon }^2}^0,{𝒲}_{BL,{\epsilon }^3}^0)$	$O\left({\epsilon }^{1/2}\right)$	${\epsilon }^{-3}$
(b2)	$Q({𝒲}_{\text{inc}}^0,{𝒲}_{BL,{\epsilon }^3}^0)$	$O\left({\epsilon }^{1/2}\right)$	${\epsilon }^{-3}$
(b3)	$Q({𝒲}_{BL,{\epsilon }^3}^0,{𝒲}_{BL,{\epsilon }^2}^0)$	$O\left({\epsilon }^{3/2}\right)$	${\epsilon }^{-3}$
(c1)	$Q({𝒲}_{BL,{\epsilon }^2}^0,{𝒲}_{\text{inc}}^0)$	$O\left({\epsilon }^2\right)$	${\epsilon }^{-2}$
(c2)	$Q({𝒲}_{\text{inc}}^0,{𝒲}_{\text{inc}}^0)$	$O\left({\epsilon }^2\right)$	no decay
(c3)	$Q({𝒲}_{BL,{\epsilon }^3}^0,{𝒲}_{BL,{\epsilon }^3}^0)$	$O\left({\epsilon }^{5/2}\right)$	${\epsilon }^{-3}$
(c4)	$Q({𝒲}_{BL,{\epsilon }^3}^0,{𝒲}_{\text{inc}}^0)$	$O\left({\epsilon }^{7/2}\right)$	${\epsilon }^{-3}$

Equations

$${\partial }_{t}v+\delta (v\cdot \nabla )v+\nabla P+b{e}_3=\nu \Delta v,$$

$${\partial }_{t}b+\delta (v\cdot \nabla )b+\frac{g}{{\rho }_0}{v}_3{\bar{\rho }}^{\prime }=\kappa \Delta b,$$

$$\nabla \cdot v=0,$$

$${\partial }_t\left(\begin{array}{c}{\hat{v}}_{p,3}\\ {\hat{b}}_p\end{array}\right)+\left(\begin{array}{cc}0& 1-\frac{p_3^2}{|p{|}^2}\\ \frac{g{\bar{\rho }}^{\prime }}{{\rho }_0}& 0\end{array}\right)\left(\begin{array}{c}{\hat{v}}_{p,3}\\ {\hat{b}}_p\end{array}\right)=0.$$

$$\omega =\pm N\frac{|p_h|}{|p|},$$

$$\frac{|p_1|}{|p|}=\frac{1}{N}\omega$$

$$x_1\sin \gamma -x_3\cos \gamma =0.$$

$${\partial }_{t}u+\delta \left(u{\partial }_x+w{\partial }_y\right)u-\sin \gamma b+{\partial }_{x}p=0,$$

$${\partial }_{t}w+\delta \left(u{\partial }_x+w{\partial }_y\right)w-\cos \gamma b+{\partial }_{y}p=0,$$

$${\partial }_{t}b+\delta \left(u{\partial }_x+w{\partial }_y\right)b+N^2(\sin \gamma u+\cos \gamma w)=0,$$

$${\partial }_{x}u+{\partial }_{y}w=0,$$

$$-i\omega \hat{u}-\sin \gamma \hat{b}+ik\hat{P}=0,$$

$$-i\omega \hat{w}-\cos \gamma \hat{b}+im\hat{P}=0,$$

$$-i\omega \hat{b}+N^2\sin \gamma \hat{u}+N^2\cos \gamma \hat{w}=0,$$

$$ik\hat{u}+im\hat{w}=0.$$

$$\frac{1}{N^2}{\omega }^2=\frac{(k\cos \gamma -m\sin \gamma {)}^2}{k^2+m^2}={\sin }^2\beta .$$

$$m^2({\sin }^2\beta -{\sin }^2\gamma )+2km\cos \gamma \sin \gamma +({\sin }^2\beta -{\cos }^2\gamma )k^2=0.$$

$$m=k\frac{\cos \gamma \sin \gamma -\cos \beta \sin \beta }{{\sin }^2\beta -{\sin }^2\gamma },m^{\prime }=k\frac{\cos \gamma \sin \gamma +\cos \beta \sin \beta }{{\sin }^2\beta -{\sin }^2\gamma }.$$

$$\left(\begin{array}{c}\hat{u}\\ \hat{w}\\ \hat{b}\end{array}\right)$$

$$\displaystyle+B\left(\begin{array}[]{c}1\\ -\frac{k}{m^{\prime}}\\ N\frac{i(k\cos\gamma-m^{\prime}\sin\gamma)}{m^{\prime}\sin\beta}\end{array}\right)e^{-iN\sin\beta t+ikx+im^{\prime}y}+c.c.,$$

$$A\frac{k}{m}+B\frac{k}{m^{\prime }}=0,$$

$${\partial }_t\left(\begin{array}{c}{\hat{u}}_p\\ {\hat{w}}_p\\ {\hat{b}}_p\end{array}\right)+\left(\begin{array}{ccc}\nu (k^2+m^2)& 0& \frac{m(k\cos \gamma -m\sin \gamma )}{k^2+m^2}\\ 0& \nu (k^2+m^2)& \frac{k(m\sin \gamma -k\cos \gamma )}{k^2+m^2}\\ N^2\sin \gamma & N^2\cos \gamma & \kappa (k^2+m^2)\end{array}\right)\left(\begin{array}{c}{\hat{u}}_p\\ {\hat{w}}_p\\ {\hat{b}}_p\end{array}\right)=0,$$

$$\Vert (\mathcal{W}-{\mathcal{W}}_{app})(t){\Vert }_{L^2}^2\le$$

$$+{\int }_0^t\Vert R_{app}(s){\Vert }_{L^2}^2\exp (\delta \Vert \nabla {\mathcal{W}}_{app}{\Vert }_{L^{\infty }}^2(t-s))ds.$$

$${\partial }_{t}u-b\sin \gamma +{\partial }_{x}p+\delta (u{\partial }_{x}u+w{\partial }_{y}u)={\nu }_0{\epsilon }^6\Delta u,$$

$${\partial }_{t}w-b\cos \gamma +{\partial }_{y}p+\delta (u{\partial }_{x}w+w{\partial }_{y}w)={\nu }_0{\epsilon }^6\Delta w,$$

$${\partial }_{t}b+u\sin \gamma +w\cos \gamma +\delta (u{\partial }_{x}b+w{\partial }_{y}b)={\kappa }_0{\epsilon }^6\Delta b,$$

$${\partial }_{x}u+{\partial }_{y}w=0,$$

$$u_{|y=0}=w_{|y=0}={\partial }_y{b}_{|y=0}=0.$$

$${\mathbb{V}}_{\sigma }:=\{(u,w,b)\in H^1({\mathbb{R}}_{+}^2{)}^3,u_{|y=0}=w_{|y=0}=0,{\partial }_{x}u+{\partial }_{y}w=0\},$$

$${\mathbb{V}}_{\sigma }^{\prime }=\text{dual space of }{\mathbb{V}}_{\sigma }.$$

$$\mathcal{W}\in C({\mathbb{R}}^{+};{\mathbb{V}}_{\sigma }^{\prime })\cap L^{\infty }({\mathbb{R}}^{+};L^2({\mathbb{R}}_{+}^2))\cap L_{loc}^2({\mathbb{R}}^{+};{\mathbb{V}}_{\sigma })$$

$$\displaystyle\|\mathcal{W}(t)\|^{2}_{L^{2}(\mathbb{R}^{2}_{+})}+2\varepsilon^{6}\nu_{0}\int_{0}^{t}\Bigg{(}\|\nabla u(s)\|^{2}_{L^{2}(\mathbb{R}^{2}_{+})}+\|\nabla w(s)\|_{L^{2}(\mathbb{R}^{2}_{+})}^{2}\Bigg{)}\,ds$$

$$+2{\epsilon }^6{\kappa }_0{\int }_0^T\Vert \nabla b(s){\Vert }_{L^2({\mathbb{R}}_{+}^2)}^{2}ds\le \Vert {\mathcal{W}}_0{\Vert }_{L^2({\mathbb{R}}_{+}^2)}^2.$$

$${\mathcal{W}}_{app}:=(u_{app},w_{app},b_{app})={\mathcal{W}}_{\text{inc}}^0+{\mathcal{W}}_{BL}+{\mathcal{W}}_{II}^1+{\mathcal{W}}_{corr},$$

$$\begin{array}{r}{\partial }_t{u}_{app}-b_{app}\sin \gamma +{\partial }_x{p}_{app}+\delta (u_{app},w_{app})\cdot \nabla u_{app}-{\nu }_0{\epsilon }^6\Delta u_{app}=O(\delta {\epsilon }^2),\\ {\partial }_t{w}_{app}-b_{app}\cos \gamma +{\partial }_y{p}_{app}+\delta (u_{app},w_{app})\cdot \nabla w_{app}-{\nu }_0{\epsilon }^6\Delta w_{app}=O(\delta {\epsilon }^2),\\ {\partial }_t{b}_{app}+u_{app}\sin \gamma +w_{app}\cos \gamma +\delta (u_{app},w_{app})\cdot \nabla b_{app}-{\kappa }_0{\epsilon }^6\Delta b_{app}=O(\delta {\epsilon }^2),\\ {\partial }_x{u}_{app}+{\partial }_y{w}_{app}=0,\end{array}$$

$${u_{app}}_{|y=0}={w_{app}}_{|y=0}={\partial }_y{b_{app}}_{|y=0}=0.$$

$${\mathcal{W}}_0={\mathcal{W}}_{app}(t=0),$$

Full text

1 Main results

Internal waves are of utmost importance in oceanic flows. They describe small departures from equilibrium in an incompressible fluid under the combined effect of density stratification and gravity. These waves are very well described when the effect of boundaries is neglected, assuming for instance that the fluid is contained in a parallelepipedic box with zero flux condition at the boundary (or equivalently in a periodic box).

Assume that at equilibrium the stratification is given by the stable profile $\bar{\rho}(x)=\bar{\rho}(x_{3})$ with $\bar{\rho}^{\prime}(x_{3})<0$. Note that, in most physical systems, the variations of $\bar{\rho}$ are very small compared to its average $\rho_{0}$, and count only for the buoyancy effect (not for inertia). Small perturbations of this equilibrium will create both a (zero divergence) velocity field $\delta v$, and a fluctuation of the density $\rho=\bar{\rho}+\delta\dfrac{\rho_{0}}{g}b$ of order $\delta\ll 1$, where $b$ is the buoyancy. The dynamics of the system is then governed by the Boussinesq equations :

[TABLE]

see [21, 4, 1, 8] for further details on the derivation of the system. Note that the viscous term $\kappa\Delta b$ in the buoyancy equation comes from the thermal dissipation (combined with the Boussinesq approximation connecting the density and the temperature, see [10]). One often considers in addition that the stratification is locally affine so that $\bar{\rho}^{\prime}$ is a constant, see again [10]. For the sake of simplicity, we will focus on this case. Note however that this assumption is not consistent with the approximation $\bar{\rho}\sim\rho_{0}$ in a large domain, especially in the whole space or in a half space.

Keeping only the leading order terms (the linear inviscid approximation) in (1.1), and taking the Fourier transform of this linear system with constant coefficients, we obtain the linear inviscid wave equation

[TABLE]

The solution in $\mathbb{R}^{3}$ can be expressed as a sum of plane waves with dispersion relation

[TABLE]

where $N=(-g\bar{\rho}^{\prime}/\rho_{0})^{1/2}$ denotes the Brunt-Väisälä frequency (see [21] for further details). By analogy with geometric optics for electromagnetic or acoustic waves (see for instance [13, 14, 16]), one can use this dispersion relation to study the “propagation” of internal waves (which actually makes sense only for wave packets and not for single plane waves, see [17]). The crucial difference here is that the frequency $\omega$ prescribes the direction of the propagation instead of the modulus of the wavelength : in 2D, on the energy level $\omega$, the wavenumber $p$ satisfies indeed

[TABLE]

and the group velocity $\nabla_{p}\omega$, which is orthogonal to $p$, makes an angle $\beta=\pm\arcsin(\omega/N)$ with respect to the horizontal (see also [6] for a mathematical description of attractors for waves with homogeneous dispersion relation of degree 0).

What we would like to understand in this paper is how these waves behave in presence of a sloping boundary. This is a major challenge from the physical point of view. Sandstrom [23] proposed indeed in 1966 the oceanic internal wave field as a possible source of the energy which is needed to activate strong mixing near sloping boundaries.

1.1 Physical predictions for the near-critical reflection

$\bullet$ The reflection of inviscid internal waves off a uniformly sloping bottom was first investigated by Philipps [20] in 1966. Because the wave frequency is related to the direction of propagation by $\omega=N\sin\beta$, preservation of $\omega$ in the reflection implies preservation of the angle $\beta$ (which is the angle between the group velocity and the horizontal, or equivalently between the phase velocity and the vertical). A simple geometric picture (see Figure 1) shows then that this will generate a focusing mechanism.

More precisely, to determine the reflection laws, we look at the solutions of the linear inviscid wave equation (1.2) in the half space delimited by the slope

[TABLE]

We seek these solutions in the form of an incident wave propagating with an angle $\beta$ with respect to the horizontal (recall that the direction of propagation is orthogonal to the wavenumber) plus a reflected wave. In order that the zero flux condition (which is the only admissible boundary condition in the inviscid regime) is satisfied on the slope, we then obtain necessary conditions on the wavenumber of the reflected wave, as well as some polarization conditions to determine the amplitude of the reflected wave. To simplify our study, we work in a 2D setting from now on.

An appropriate system of coordinates to express these conditions is $(x,y)$ where $x$ is the abscissa along the boundary and $y$ the distance to the boundary (see Figure 2).

We will denote by $(k,m)$ the corresponding Fourier variables, and by $(u,w)$ the tangential and normal components of the velocity. In these coordinates, the fully nonlinear system (1.1) becomes

[TABLE]

while the linear inviscid wave system (1.2) can be restated

[TABLE]

Since it has to lift the boundary condition $w_{|y=0}=0$, the reflected wave should have the same horizontal wave number $k$ and time frequency $\omega$ as the incident wave, but a different vertical wave number $m^{\prime}$, which also satisfies the dispersion relation

[TABLE]

The equation for $m$ and $m^{\prime}$ is therefore

[TABLE]

The roots are

[TABLE]

Note in particular that, as $\beta\to\gamma$, $|m^{\prime}|\to\infty$.

The elementary solution of (1.4) supplemented with the zero flux condition is therefore of the form

[TABLE]

where $c.c.$ denotes the complex conjugate, with the polarization condition

[TABLE]

meaning that the amplitude of the reflected wave is $O(\frac{m^{\prime}}{m})$.

The most effective situation for boundary mixing arises therefore when an oncoming wave reflects off a bottom slope with angle $\gamma$ which nearly matches the angle of wave propagation, namely $\beta=\gamma-{\varepsilon}^{2}$ with ${\varepsilon}\ll 1$. At the critical angle $\beta=\gamma$, the analytic theory of internal waves reflecting off a uniformly sloping bottom predicts that the reflected wave has infinite amplitude and infinitesimal wavelength. These unphysical results signal the failure of the idealizations (namely, linear waves and inviscid fluid).

$\bullet$ The viscous and nonlinear effects associated to the near-critical reflection of internal waves have been studied by Dauxois and Young [10] in 1999. The main idea is that, even though the original perturbation generating the incident wave is small, the amplitude of the reflected wave is enhanced by the focusing mechanism, so that the nonlinear coupling between the incident and the reflected wave might no longer be negligible. Moreover, the dissipation is enhanced because of the fact that the normal wavenumber of the reflected wave is much bigger than the original wavenumber.

Actually, when one considers the viscous case, the analysis is much more delicate as the system of equations has to be supplemented with more boundary conditions : since the equations for $u$, $w$ and $b$ are parabolic, one should impose one boundary condition for each one of these quantities. From the physical point of view, it is natural to prescribe a Dirichlet boundary condition on both components of the velocity (no slip condition), and a Neumann condition on the buoyancy (no diffusive flux through the slope, see [10] for a discussion on this point). Even at the linear level, i.e. for the viscous wave system

[TABLE]

with $p=(k,m)$, it is therefore much more complicated to compute elementary solutions satisfying the suitable boundary conditions since there are three matching conditions instead of one. In their paper, Dauxois and Young actually discarded the boundary condition on the buoyancy (which is justified for instance if $\kappa=0$), and chose the distinguished scaling $\nu={\varepsilon}^{6}$. They found that elementary solutions can be then decomposed as the sum of

• •

the incident wave;

• •

two boundary layer terms with exponential decay in $y/{\varepsilon}^{2}$.

In other words, this means that one can find two values $m_{1},m_{2}$ (distinct from $m$, and having a non negative imaginary part) such that the matrix in (1.7) has $i\omega=iN\sin\beta$ as an eigenvalue. These two eigenvalues $m_{1}$, $m_{2}$ bifurcate from $m^{\prime}$ when we add the small viscosity $\nu=O({\varepsilon}^{6})$, and they are of order $O(1/{\varepsilon}^{2})$ like $m^{\prime}$. The matching conditions for $u$ and $w$ at $y=0$ provide then polarization conditions for the amplitudes of these two damped reflected waves.

In the boundary layer, for $\delta$ small enough, we expect the nonlinear term to be weak so that it should not destroy the wave structure. This means that the nonlinearity can be dealt with as a weak perturbation of the linear evolution. At leading order, we then expect it will produce two new modes, one which is oscillating with frequency $2\omega$ and a mean flow. The amplitudes of these modes should be a small correction to the boundary layer, but which has non zero trace on the boundary! In their paper, Dauxois and Young lifted this (small) remaining trace by adding

• •

a second harmonic propagating in the outer domain provided that $\frac{2\omega}{N}\leq 1$;

• •

a mean flow (which is not written explicitly as it is expected to stay localized close to the boundary).

As a result of this construction, they obtain an approximate solution in the sense that it satisfies the equation (1.1) with suitable scalings for the parameters of the system, together with the boundary conditions $u_{|y=0}=w_{|y=0}=0$, up to remainders which should be either smaller in energy, or of same order but fast oscillating (thus converging weakly to 0). We will say that the approximation is consistent. This Ansatz shows in particular that there is a critical amplitude where the nonlinearity, although small, is no longer negligible even in the outer domain. This is made evident by the generation of the second harmonic.

$\bullet$ This scenario has been validated by experimental visualization of the reflection process. The lab experiments conducted by physicists, especially in the group of Dauxois [9, 12], have indeed confirmed that the nonlinearity plays a key role in the boundary layer.

Dauxois, Didier, and Falcon introduced in [9] the Schlieren technique to study the spatiotemporal evolution of the internal waves reflection close to the critical reflection. The internal wave, producing density disturbances, causes lines to distort, this distorting line pattern being recorded by a camera. Note that this experimental technique is sensitive to the index gradient, and therefore to the density gradient. The dynamics of isopycnals is then found in good qualitative agreement with the theory. Moreover, this experiment confirms the theoretically predicted scenario for the transition to boundary-layer turbulence responsible for boundary mixing: the growth of a density perturbation produces a statically unstable density field which then overturns with small-scale fluctuations inside.

Peacock and Tabaei then presented in [18] the first set of experimental visualizations (using also the digital Schlieren method) that confirm the existence of radiated higher-harmonic beams. For arrangements in which the angle of propagation of the second harmonic exceeds the slope angle, radiated beams are visualized. When the propagation angle of the second harmonic deceeds the slope angle no radiated beams are detected, as the associated density gradient perturbations are too weak for the experimental method. The case of a critical slope is also reported.

Quantitative results have finally been obtained in [12]. Experiments were carried out in the Coriolis platform, in Grenoble, filled with salted water. The large scale of the facility allows to strongly reduce the viscous dissipation along wave propagation and quantitative results are obtained thanks to high-resolution Particle Image Velocimetry measurements. Generation of the second and third harmonic frequencies is clearly demonstrated in the impact zone. Although these harmonics are almost invisible from the instantaneous velocity field, they are very clearly apparent after the filtering procedure. These experiments also provide evidence that harmonics with frequency higher than $N$ cannot propagate and remain trapped near the slope.

1.2 Mathematical description of the near-critical reflection

The mathematical analysis of the near-critical reflection follows essentially the same lines, but requires a more careful treatment of some delicate points which we will explain now.

First of all, the consistency of the approximation is not sufficient to deduce that the solution of the original problem will be close to the approximate solution. This issue is related to the possible instabilities of the system : a small error on the equation (in the form of a source term $R_{app}$) could then generate important deviations of the dynamics. This means that, in order to prove that the matched asymptotic expansion $\mathcal{W}_{app}$ provides a good approximation of what happens in reality, one needs to prove some stability for the original weakly nonlinear system (1.1). The classical way of doing so is to establish energy inequalities : for the Navier-Stokes equations, we typically expect that the growth of the $L^{2}$ norm $\|\mathcal{W}-\mathcal{W}_{app}\|_{L^{2}}$ should be controlled by the Lipschitz norm $\delta\|\nabla\mathcal{W}_{app}\|_{L^{\infty}}$, provided that there is no error (at all) on the boundary conditions.

[TABLE]

However this has many implications. More details will be provided in Section 4.2.

• (i)

We need to work with solutions of finite energy (or at least such that $\mathcal{W}-\mathcal{W}_{app}$ has finite energy), which is not the case of (a finite sum of) plane waves!

• (ii)

The remainder has to be small (in energy), and actually smaller than the second harmonic if we would like to prove that this term is physically relevant.

• (iii)

The time interval on which the approximation is accurate depends on the Lipschitz (or at least the $L^{\infty}$) norm of the approximate solution.

Points (ii) and (iii) are somehow technical, they can be tackled by tracking the dependency of $\mathcal{W}_{app}$ and $R_{app}$ with respect to the different parameters. We will also need to construct some additional correctors to lift the boundary conditions. We refer to Section 3 for these technical details.

Point (i) is more tricky. It is related to the fact that a “plane wave” is not a good physical object (although it is considered as a very basic object in physics!) This has already been mentioned in the introduction when talking about propagation of wave and group velocity. These notions do not make sense for a single plane wave, but rather for a wave packet. Plane waves just provide a tool to decompose these wave packets on elementary objects (just like the Fourier transform). A very pedagogical discussion on these topics can be found for instance in [17]. In all the sequel of this paper, we will thus consider wave packets, with energy density concentrated close to the frequency $\omega_{0}=\sin\gamma$ in order to see the effects of criticality. Of course, as long as we study linear equations, looking at a wave packet (i.e. at a superposition of plane waves) does not introduce any additional difficulty for the matched asymptotic expansion. But this is no longer true when the nonlinearity enters into the game.

Let us then go back to the construction of the approximate solution $\mathcal{W}_{app}$, as there are also here some points which need to be clarified.

• (iv)

In the construction by Dauxois and Young [10], one boundary condition (the one on the buoyancy $b$) is discarded, supposedly because it should produce effects of higher order. This is however not completely clear from the computations given in the paper (see the comments p. 283). Section 2 of the present paper will provide a very systematic construction for linear boundary layers. We refer to [11, 7] for a presentation of the method. We will see in particular that the solution to the viscous wave equation (1.7) involves a superposition of boundary layers of different sizes, and that there is no overdetermination of the problem if we take into account all these contributions.

• (v)

Since we will consider wave packets, we will need to understand how to deal with the nonlinear interactions in the boundary layer for a general superposition of waves (see Section 3). Section 3 will provide a careful analysis of the generation of both the rectified (or mean) flow and the second harmonic.

• (vi)

The rectified flow in [10] (defined as the component with $\omega=0$ generated by the weak nonlinearity) is neglected since it has to vanish far from the boundary (see p. 282). Although in [10] there is no quantitative estimate of the decay rate, nor of the energy contained in this rectified flow (which could be more energetic than the second harmonic for instance), this is in agreement with our construction. We have indeed two boundary layer mean flows (vanishing far from the boundary), and the third boundary condition can be lifted in a very crude way (by adding a small corrector which has no physical relevance).

1.3 Notation

Our construction of an approximate solution will involve many terms, which all play a different role. For the reader’s convenience, we have gathered here the notations that will be used throughout the paper, and by doing so we roughly sketch the main steps of the construction.

• •

We denote by $\mathcal{W}$ the vector $\begin{pmatrix}u\\ w\\ b\end{pmatrix}$;

• •

The first step of our construction is the definition of a solution of the linear problem. This solution will be the main order term in the approximate solution, and therefore we will denote it by $\mathcal{W}^{0}$;

• •

The first order correctors, which lift the nonlinearity $\delta Q(\mathcal{W}^{0},\mathcal{W}^{0})$, will be denoted by $\mathcal{W}^{1}$.

The notations $\mathcal{W}^{j}$ allow us to distinguish between terms of different orders within the approximate solution. But we also need to specify the nature of each term within a given order. For instance, $\mathcal{W}^{0}$ will be the sum of an incident wave packet, of a boundary layer term localized within a region of width ${\varepsilon}^{2}$ of the boundary $y=0$, and of a boundary layer term localized within a region of width ${\varepsilon}^{3}$ of the boundary $y=0$. Therefore we adopt the following notation:

• •

The incident wave packet will be denoted $\mathcal{W}^{0}_{\text{inc}}$;

• •

The boundary layer terms localized within a region of width ${\varepsilon}^{2}$ (resp. ${\varepsilon}^{3}$) of the boundary will be denoted $\mathcal{W}_{BL,{\varepsilon}^{2}}^{j}$ (resp. $\mathcal{W}_{BL,{\varepsilon}^{3}}^{j}$), where the superscript $j$ refers to the order of the term they belong to;

• •

The mean flow will be denoted by $\mathcal{W}_{MF}^{1}$;

• •

The second harmonic will be denoted by $\mathcal{W}_{II}^{1}$.

1.4 The stability result

From now on, we choose $N=1$ in order to alleviate the notation.

We consider the Cauchy problem associated to the Boussinesq equations in $\mathbb{R}^{2}_{+}=\mathbb{R}\times\mathbb{R}^{+}$,

[TABLE]

endowed with the boundary conditions

[TABLE]

Note that this is a just rewriting of the Boussinesq system (1.1) in the slope coordinates, nothing but the viscous version of system (1.3). As discussed at the beginning of Section 3, after a complete linear analysis of all the possible regimes induced by the different parameters of the system in Section 2, the nonlinear system will be treated according to the scalings by Dauxois and Young in [10], where, in particular, the size of the viscosity is $\nu=\nu_{0}\varepsilon^{6},\,\kappa=\kappa_{0}\varepsilon^{6}$.

Our main result is a stability estimate for the approximate solution $\mathcal{W}_{app}=(u_{app},w_{app},b_{app})$, i.e. an $L^{2}$ estimate of the difference between $\mathcal{W}_{app}$ and the weak solutions $\mathcal{W}=(u,w,b)$ to (1.9)-(1.10). Before providing the stability theorem, the global existence of weak solutions to the Cauchy problem for system (1.9) with boundary conditions (1.10) is stated below.

We will use the following notation.

[TABLE]

Proposition 1.1**.**

Let $\mathcal{W}_{0}=(u_{0},w_{0},b_{0})$ a divergence free $L^{2}$ initial data. Then there exists a unique global weak solution

[TABLE]

to system (1.9)-(1.10), which satisfies the following energy inequality for all $t\geq 0$,

[TABLE]

The proof of this theorem follows by adapting the result of global weak (Leray) solutions to the incompressible Navier-Stokes equation in general domains in [5]. The crucial point is the conservation of energy, represented by the energy inequality (1.12), which can be simply obtained in the usual way, by taking the scalar product of (1.9) with $\mathcal{W}$.

For the sake of completeness, we will provide a sketch of the argument in the Appendix.

We now state our main result.

Theorem 1.1**.**

*Consider the Boussinesq equations (1.9) in the scaling by Dauxois and Young in $\mathbb{R}_{+}^{2}=\mathbb{R}\times\mathbb{R}_{+}$, with boundary conditions (1.10).

Then there exists a vector field*

[TABLE]

where $\mathcal{W}^{0}_{\text{inc}},\,\mathcal{W}_{BL},\,\mathcal{W}^{1}_{II},\,\mathcal{W}_{corr}$ are respectively an incident wave packet, a boundary layer, a second harmonic wave packet and an additional correction term, which is an approximate solution, in the sense that

[TABLE]

where the remainders $O(\delta\varepsilon^{2})$ have to be understood in the sense of the $L^{2}(\mathbb{R}^{2}_{+})$ norm, and endowed with the boundary conditions

[TABLE]

Furthermore, denoting by $\mathcal{W}$ the unique weak solution to the Cauchy problem associated with system (1.9)-(1.10), with initial data

[TABLE]

we have the following stability estimate:

[TABLE]

Remark 1.2**.**

Let us add some comments on Theorem 1.1.

• •

The first part of the statement could be reformulated saying that $\mathcal{W}_{app}$ is a consistent approximate solution, which means that it satisfies the equations of system (1.9)-(1.10) modulo a remainder (of order $\delta\varepsilon^{2}$ in $L^{2}$). In the second part, we infer that the consistent approximate solution is also stable in the sense of the $L^{2}$ norm, meaning that the difference in $L^{2}$ between this approximate solution and the unique weak solution to the original system (1.9)-(1.10) with initial data $\mathcal{W}_{app}(t=0)$ is smaller than a remainder of size $\delta\varepsilon^{2}$ in $L^{2}$. Note that, as widely discussed in the rest of the paper, the boundary conditions (1.10) need to be exactly satisfied by the approximate solution in order to establish the stability inequality (1.13). A consistent approximate solution is indeed already given by the sum of first three terms $\mathcal{W}^{0}_{\text{inc}}+\mathcal{W}_{BL}+\mathcal{W}_{II}^{1}$, while to achieve stability we need to add another corrector $\mathcal{W}_{corr}$, whose explicit expression is given in Proposition 3.1.

• •

We point out that far from the boundary, the leading order term of the approximate solution $\mathcal{W}_{app}$ is represented by the incident wave packet $\mathcal{W}^{0}_{\text{inc}}$.

• •

Theorem 1.1 could also be stated in a more general fashion, by taking into account any incident wave packet (not only the ones of the critical regime by Dauxois and Young). The construction of the approximate solution can be done exactly in the same spirit and it should be easier, the critical case being the most difficult one to handle.

More details on the terms involved in the expression of the approximate solution $\mathcal{W}_{app}$ are provided in Lemma 2.7 and Proposition 3.1.

The proof of Theorem 1.1 relies on a quite accurate quantification of all the sizes of the terms and the remainders of the approximate solution in $L^{2}$ and $L^{\infty}$. These computations will be provided in details in Section 3. A general construction of the profiles and sizes of the boundary layers, depending on their time frequency and tangential wave number, is provided in Section 2. The stability inequality is established in Section 4.

2 Linear viscous boundary layers

The purpose of this section is to provide a systematic description of boundary layers in the linear case. In particular, we will explain how boundary layer sizes and profiles can be computed, which boundary conditions can be lifted, and we will derive the asymptotic behavior of the boundary layer sizes in different regimes, including the case of critical reflection. In this case, we will also derive the expansions of the linear boundary layer operators, that will allow us to perform the weakly nonlinear analysis of Section 3.

As explained in the introduction, the first item in the construction of an approximate solution is an incident wave $(u_{i},w_{i},b_{i})$ (or an incident wave packet, which is an infinite linear superposition of incident waves). This incident wave does not solve the boundary conditions

[TABLE]

and it is expected that boundary layers take place close to the wall $y=0$ to lift the traces of $(u_{i},w_{i},{\partial}_{y}b_{i})$. Therefore the purpose of our analysis in the present section is the following: we seek an exact solution of the linear system

[TABLE]

endowed with the boundary conditions

[TABLE]

where $k,\omega\in\mathbb{R}$ are prescribed, and $\mathfrak{u},\mathfrak{w},\mathfrak{b}\in\mathbb{C}$. We refer to the work of Gérard-Varet and Paul [11] for a presentation of the methodology developed here, and to the paper [7] by the last two authors for an application to a case in which several boundary layers with different sizes co-exist. Since equation (2.2) is linear and has constant coefficients, it is natural to look for modal solutions, i.e. (linear combinations of) functions of the form

[TABLE]

with $\Re(\lambda)>0$. When we plug this expression into system (2.2), we obtain

[TABLE]

where

[TABLE]

As a consequence, a modal function of the form (2.4) is a non trivial solution of (2.2) if an only if $\det A_{\nu,\kappa,\omega,k}(\lambda)=0$ and $(U,W,B,P)\in\ker A_{\nu,\kappa,\omega,k}(\lambda)\setminus\{0\}$. The corresponding solution is a boundary layer mode if $\Re(\lambda)\gg 1$, and the size of the boundary layer is then $(\Re(\lambda))^{-1}$. Furthermore, the number of boundary conditions that can be lifted by the boundary layer is equal to the dimension of the vector space

[TABLE]

where $(U,W)$ lift the Dirichlet boundary conditions for the velocity field, while $-\lambda B$, which is the amplitude of the derivative in $y$ of expression (2.4), counts for the no flux condition on the buoyancy $\partial_{y}b=0$ at $y=0$. Now, it can be easily checked that $\det A_{\nu,\kappa}(\omega,k,\lambda)$ is a polynomial of degree 6 in $\lambda$, that has therefore 6 complex roots, counted with multiplicity. We now investigate the asymptotic values of these roots depending on the viscosity and the slope criticality. We will prove in particular that we can usually lift three boundary conditions thanks to this linear boundary layer.

2.1 Sizes of the boundary layers in terms of the viscosity and slope criticality

A straightforward computation leads to

[TABLE]

Note that in the case $\nu=\kappa=0$ and $\lambda=-im$, we retrieve the dispersion relation (1.5). In all regimes considered below, we will make the following assumptions: $0<\nu\ll 1$, $0<\kappa\ll 1$, $|k|\leq c_{0}$, $c_{0}^{-1}\leq|\cos^{2}\gamma-\omega^{2}|\leq c_{0}$ for some fixed constant $c_{0}$. We introduce a criticality parameter $\zeta:=\omega^{2}-\sin^{2}\gamma$, and we will investigate both the cases $|\zeta|\ll 1$ and $|\zeta|\gtrsim 1$. To simplify the discussion, we will also assume that $\nu$ and $\kappa$ are of the same order of magnitude. However, we will not systematically assume that we are in a critical setting. In other words, we will not assume any relationship between $\nu$ and $\zeta$. It follows from the above assumptions that when $\omega$ and $k$ are bounded away from zero, the roots of $\det A_{\nu,\kappa,\omega,k}$ are asymptotically equivalent to the ones of

[TABLE]

However, for further purposes, we will also need to consider the case $|\omega|\ll 1$, $|k|\ll 1$, which we will do in a separate paragraph.

Let us now explain our general strategy to find the boundary layer sizes. We first look for roots with the smallest possible absolute value. In the case $|k|\geq c_{0}^{-1}$, since all coefficients of $P_{\nu,\kappa,\omega,k}$ are bounded and the coefficient of order zero is bounded away from zero, the smallest roots will be of size 1, and all coefficients of order $o(1)$ can be neglected. We then look for roots $\lambda$ with a large absolute value, balancing the higher order terms, whose coefficients are smaller, with (some of) the lower order ones. This leads to an Ansatz of the form $|\lambda|\propto\nu^{-\alpha}$ for some $\alpha>0$. We then plug this Ansatz into $P_{\nu,\kappa,\omega,k}$, check that the terms that have not been taken into account can indeed be neglected, and derive an effective equation for $\lambda$. While looking for possible asymptotic behaviors for $\lambda$, it is also useful to keep in mind that the product of all roots is $-(\nu\kappa)^{-1}(\cos^{2}\gamma-\omega^{2})$.

2.1.1 Case $|k|\geq c_{0}>0$, $|\omega|\geq c_{0}>0$

Case $|\zeta|\gtrsim 1$ (non-critical reflection)

In this case, it can be easily checked that the roots of $P_{\nu,\kappa,\omega,k}$ can be classified in the following way:

• •

There are two roots of size 1, which therefore do not correspond to boundary layers, and which are approximate roots of the equation

[TABLE]

These two roots converge towards pure imaginary values in the limit $\nu,\kappa\to 0$ when $|\omega|$ remains below the critical value $\omega_{c}=1$. They then correspond to the incident and the reflected wave in the non-critical case (compare the above equation with the dispersion relation (1.5) or with (1.6)). When $|\omega|>1$, the two roots have a non-zero real part: one of them (say $\lambda_{1}$) has a negative real part (and must therefore be discarded), and the other (say $\lambda_{2}$) has a positive real part. It corresponds to an exponentially decaying mode, which is however not confined to a boundary layer.

• •

There are four roots $\lambda$ such that $|\lambda|=O(\nu^{-1/2})$ (remember that we have assumed that $\nu$ and $\kappa$ are of the same order), which are approximate roots of the equation

[TABLE]

Among these roots, two have a strictly positive real part, and two have a negative real part.

Case $\nu^{1/3}\ll|\zeta|\ll 1$ (critical reflection with small diffusivity)

In this case, there are four types of roots:

• •

There is one root of size 1, say $\lambda_{1}$, such that

[TABLE]

This root corresponds to the incident wave.

• •

There is one root of size $\zeta^{-1}$, say $\lambda_{2}$, which satisfies

[TABLE]

Plugging this expression back into the equation $A_{\nu,\kappa,{\omega},k}(\lambda_{2})=0$, it can be proved that $\Re(\lambda_{2})$ is negative and of size $\nu/\zeta^{4}$.

• •

There are two roots of size $|\zeta/\nu|^{1/2}$, say $\lambda_{3}$ and $\lambda_{4}$, which are approximate solutions of

[TABLE]

Among these two roots, exactly one (say $\lambda_{3}$) has a strictly positive real part of size $|\zeta/\nu|^{1/2}$, and the other one has a strictly negative real part.

• •

There are two roots of size $\nu^{-1/2}$, which are approximate solutions of

[TABLE]

Among these two roots, exactly one (say $\lambda_{5}$) has a strictly positive real part of size $\nu^{-1/2}$, and the other one has a strictly negative real part.

Case $\nu^{1/3}\sim|\zeta|$ (critical reflection in the scaling of Dauxois and Young)

This case is very similar to the one above. The only difference lies in the fact that the roots $\lambda_{2},\lambda_{3}$ and $\lambda_{4}$ are now all of size $\nu^{-1/3}$, and are approximate solutions of

[TABLE]

One can check that this equation has exactly two roots with positive real part (see Remark 2.1 below).

Case $|\zeta|\ll\nu^{1/3}$ (critical reflection with large diffusivity)

Once again, this case is similar to the two cases above. The three roots $\lambda_{2},\lambda_{3}$ and $\lambda_{4}$ are all of size $\nu^{-1/3}$, and are now approximate solutions of

[TABLE]

The conclusion remains the same.

Remark 2.1** (Number of roots with positive real part).**

In the analysis above, we have counted by hand the number of roots of $\det A_{\nu,\kappa,\omega,k}$ with positive real part in some specific cases. Les us now prove that this number remains constant as the parameter $\zeta$ varies. Assume that the parameters $\nu,\kappa,\gamma$ and $k$ are fixed, with $\nu$ and $\kappa$ small, and consider the roots $\lambda_{i}$, $i=1,...6$ as functions of the criticality parameter $\zeta$. Then it can be easily checked, using the previous arguments, that the curves $\zeta\mapsto\lambda_{i}(\zeta)$ do not cross. Indeed, if there were a crossing, it would obviously occur for two eigenvalues within the same regime, and a quick look at each of the regimes described above ensures that this cannot happen. As a consequence, each function $\zeta\mapsto\lambda_{i}(\zeta)$ is continuous (and even $\mathcal{C}^{1}$). Furthermore, the curves do not cross the imaginary axis: indeed, if a root $\lambda$ is imaginary, then we deduce that

[TABLE]

and therefore $\lambda=\pm|k|$ and $\lambda\in\mathbb{R}$: contradiction. As a consequence, each curve $\zeta\mapsto\lambda_{i}(\zeta)$ remains in the same complex half-plane. In particular, $\lambda_{2},\lambda_{3}$ and $\lambda_{5}$ always have a positive real part, and $\lambda_{4}$ and $\lambda_{6}$ always have a negative real part. The eigenvalue $\lambda_{1}$ corresponds to the incident wave and is therefore discarded.

Note that an eigenvector corresponding to a root $\lambda$ is

[TABLE]

We now go back to each of the cases above, and check that the vectors

[TABLE]

are linearly independent. As a consequence, we can always lift three boundary conditions.

2.1.2 Case $|k|\ll 1$, $|\omega|\ll 1$

In the following sections, we will also need to investigate the case when $|k|$ and $|\omega|$ are small. More specifically, we will be interested in the scaling $|k|\lesssim\nu^{1/3}$, $|\omega|\lesssim\nu^{1/3}$. Going back to the expression of $\det A_{\nu,\kappa,\omega,k}$, we obtain the following asymptotic behaviors for the roots:

• •

There are two roots of size $\nu^{1/3}$, say $\lambda_{1}$ and $\lambda_{2}$, such that if $|k|\gtrsim\nu^{1/3}$

[TABLE]

Note that one of these roots (say $\lambda_{1}$) has a negative real part, the other one (say $\lambda_{2}$) a positive real part. However, the root $\lambda_{2}$ does not correspond to a boundary layer because $\Re(\lambda_{2})\ll 1$. Computations lead to

[TABLE]

in the scaling $|k|\gtrsim\nu^{1/3}$.

• •

There are four roots of size $O(\nu^{-1/2})$, say $\lambda_{i}$ for $i\in\{3,\cdots,6\}$ which satisfy

[TABLE]

Among these, two have a strictly positive real part (say $\lambda_{3}$ and $\lambda_{5}$), and two have a strictly negative real part.

Once again, we can theoretically lift three boundary conditions. Note however that the vector $X_{\lambda_{2}}\exp(-\lambda_{2}z)$ has a very slow decay, which leads to several complications: first, if $|k|\ll\nu^{1/3}$ and $|\omega|\gtrsim\nu^{1/3}$, the asymptotic given for $\lambda_{2}$ above may become invalid, and computing an exact rate of decay in all cases becomes quite technical. Furthermore, in the rest of the paper, we will consider wave packets, i.e. superposition of plane waves. This entails that the boundary layer part of the solution corresponding to the rate $\lambda_{2}$ will behave like

[TABLE]

where $\varphi\in\mathcal{C}^{\infty}_{0}(\mathbb{R})$. The Plancherel theorem entails that the $L^{2}_{x,y}$ norm of this quantity is bounded by

[TABLE]

It is not clear that the right-hand side is finite, unless we add further assumptions on $\varphi$. *Therefore, in the rest of this paper, we will discard the part of the solution that should be handled by the eigenvalue $\lambda_{2}$ in this regime. * As a consequence, we will only be able to lift two boundary conditions in this case.

Let us sum up the results of this paragraph in the following table.

Case $|k|\gtrsim 1$:

• •

If $|\zeta|\gtrsim 1$, there is one reflected wave (lifting the vertical boundary condition) and one boundary layer of size $\nu^{1/2}$, lifting the two remaining boundary conditions;

• •

If $\nu^{1/3}\ll|\zeta|\ll 1$, there are three boundary layers, of sizes $\zeta^{4}/\nu$, $(\nu/\zeta)^{1/2}$, and $\nu^{1/2}$;

• •

If $|\zeta|\lesssim\nu^{1/3}$, there are two boundary layers, of sizes $\nu^{1/3}$ and $\nu^{1/2}$. 2. 2.

Case $|k|\lesssim\nu^{1/3},\ |\omega|\lesssim\nu^{1/3}$:

There is one boundary layer of size $\nu^{1/2}$, lifting two (out of three) boundary conditions.

2.2 Boundary operators

We now gather the results and observations from the preceding paragraph.

$\bullet$ Assume that we are in the regime $|k|\gtrsim 1$. Let $(\mathfrak{u},\mathfrak{w},\mathfrak{b})\in\mathbb{C}^{3}$, and consider the linear system with unknowns $(a_{2},a_{3},a_{5})$

[TABLE]

Define the function

[TABLE]

Then by construction, we have the following result:

Lemma 2.2** (Critical reflection for an oscillating flow).**

Assume that there exists a constant $c_{0}>0$ such that $c_{0}^{-1}\leq|k|\leq c_{0}$, $c_{0}^{-1}\leq|\cos^{2}\gamma-\omega^{2}|\leq c_{0}$. Assume that $|\zeta|\ll\nu^{1/4}$ and that $\nu\ll 1$, $\kappa\ll 1$ with $\nu/\kappa\propto 1$.

Let $\mathfrak{u},\mathfrak{w},\mathfrak{b}\in\mathbb{C}$. Then there exists an exact solution $(u,w,b)$ of the linear system (2.2) with:

[TABLE]

This solution is given by

[TABLE]

where the coefficients $a_{j}$ are determined by the system (2.8). It is exponentially small outside a small boundary layer localized in the vicinity of $y=0$.

We have a similar result in the case of a non-critical reflection; the main difference lies in the fact that the solution is the sum of a boundary layer and of a flow that is not confined to a small region.

Lemma 2.3** (Non-critical reflection).**

Assume that there exists a constant $c_{0}>0$ such that $c_{0}^{-1}\leq|k|\leq c_{0}$, $c_{0}^{-1}\leq|\cos^{2}\gamma-\omega^{2}|\leq c_{0}$. Assume that $|\zeta|\gtrsim 1$ and that $\nu\ll 1$, $\kappa\ll 1$ with $\nu/\kappa\propto 1$.

Let $\mathfrak{u},\mathfrak{w},\mathfrak{b}\in\mathbb{C}$. Then there exists an exact solution $(u,w,b)$ of the linear system (2.2) with:

[TABLE]

This solution is given by

[TABLE]

where the coefficients $a_{j}$ are determined by the system (2.8). It is composed of a boundary layer of size $\nu^{\frac{1}{3}}$ and of a (possibly decaying) reflected wave.

$\bullet$ We now address the case of a non-oscillating flow, i.e. $|k|\lesssim\nu^{1/3}$, $|\omega|\lesssim\nu^{1/3}$. As stressed in the previous subsection, in this case we only lift two boundary conditions. Looking at system (2.8) and recalling that $\lambda_{2}=O(k)=O(\nu^{1/3})$, $\lambda_{3},\lambda_{5}=O(\nu^{-1/2})$, we see that the best approximation of the actual linear solution is to choose $a_{2},a_{3}$ and $a_{5}$ so that $a_{2}=0$ and $a_{3}$ and $a_{5}$ satisfy

[TABLE]

We obtain the following result:

Lemma 2.4** (Non-critical reflection for a non-oscillating flow).**

Assume that $\nu\ll 1$, $\kappa\ll 1$ with $\nu/\kappa\propto 1$, and that $|k|\lesssim\nu^{1/3}$, $|\omega|\lesssim\nu^{1/3}$.

Let $\mathfrak{u},\mathfrak{b}\in\mathbb{C}$. Then there exists an exact solution $(u,w,b)$ of the linear system (2.2) such that

[TABLE]

This solution is given by

[TABLE]

where the coefficients $a_{3},a_{5}$ are determined by system (2.10). Furthermore, the remaining trace satisfies

[TABLE]

2.3 Superposition of boundary layers in the Dauxois-Young scaling

In the rest of this section, we focus on the case of a critical reflection in the scaling of Dauxois and Young [10], i.e. we assume that $\nu^{1/3}\sim|\zeta|$. We recall that in this case, there is a boundary layer of size $\nu^{1/3}$, which is described in [10], and another one much smaller, of size $\nu^{1/2}$, which is absent from [10]. We now describe the interplay between these two boundary layers. More precisely, we compute asymptotic values for the coefficients $a_{j}$ defined in (2.8). We also switch to the notation that will be used in the rest of this paper, and that was introduced in [10]: more precisely, we take $\nu=\nu_{0}{\varepsilon}^{6}$, $\kappa=\kappa_{0}{\varepsilon}^{6}$, where ${\varepsilon}\ll 1$ is a small parameter, and $\nu_{0},\kappa_{0}>0$ are constant and independent of ${\varepsilon}$. We therefore refer to the boundary layer of size $\nu^{1/3}$ (resp. $\nu^{1/2}$) as the ${\varepsilon}^{2}$ (resp. ${\varepsilon}^{3}$) boundary layer.

Definition 2.5**.**

In the rest of this section, we denote by $\bar{\Lambda}_{2},\bar{\Lambda}_{3},\bar{\Lambda}_{5}$ the complex numbers with positive real parts defined by

[TABLE]

A straight-forward linear analysis shows that there exists $\bar{A}_{2},\bar{A}_{3},\bar{A}_{5}\in\mathbb{C}\setminus\{0\}$ such that

[TABLE]

where the coefficients $\bar{A}_{2},\bar{A}_{3},\bar{A}_{5}$ satisfy the linear system

[TABLE]

where $c:=\kappa_{0}/\nu_{0}>0$. As a consequence, if we isolate the ${\varepsilon}^{2}$ boundary layer part, namely

[TABLE]

then

[TABLE]

while

[TABLE]

Remark 2.6**.**

Note that the boundary condition on the normal component of the velocity is mostly handled by the boundary layer of size ${\varepsilon}^{2}$, which is not the case for the tangential component. There is indeed an intricate relationship between the ${\varepsilon}^{2}$ and the ${\varepsilon}^{3}$ boundary layers, and the boundary term $\mathfrak{u}$ is in fact distributed among the two boundary layers.

An (almost) exact solution to the linear system (2.2)-(2.1) in the case $c_{0}^{-1}\leq|k|\leq c_{0}$, $|\zeta|\ll{\varepsilon}^{3/2}$ is provided by the sum of incident wave and the boundary layers, i.e.

[TABLE]

where the incident wave is $\mathcal{U}_{\text{inc}}=(\mathfrak{u},\,\mathfrak{w},\,\mathfrak{b})\exp(i(kx+my-\omega t))$, and the boundary layer is

[TABLE]

Note that the boundary conditions are satisfied without any error, since the trace of the incident wave is balanced by the boundary layer. The only source of error comes from the action of the diffusion operator on the incident wave.

2.4 Definition and size of the linear solution

Since we wish to work with finite energy solutions, we will consider wave packets rather than plane waves, i.e. functions of the form

[TABLE]

where $\widehat{A}\in\mathcal{S}(\mathbb{R}^{2})$, and $X_{k,m}:=\begin{pmatrix}1\\ -\frac{k}{m}\\ \frac{i(k\cos\gamma-m\sin\gamma)}{m\omega_{k,m}}\end{pmatrix}$. The boundary layer is then

[TABLE]

where $\mathcal{B}_{\omega,k}$ can be $\mathcal{B}_{\omega,k}^{c}$, $\mathcal{B}_{\omega,k}^{nc,BL}$, or $\mathcal{B}_{\omega,k}^{no}$, depending on the regime under consideration.

In order that $\mathcal{W}^{0}_{\text{inc}}$ is an incident wave packet, we also need that the group velocity is oriented towards the slope. Since the time frequency is given by

[TABLE]

then the group velocity reads

[TABLE]

We thus choose the sign of $\omega_{k,m}$ in such a way that

[TABLE]

where $e_{y}$ is the unit vector along the normal axis $y$.

In order to focus on the critical case, we will choose functions $\widehat{A}$ which are supported close to $\pm(k_{0},m_{0})$, where $(k_{0},m_{0})$ are such that

[TABLE]

More precisely, we will take

[TABLE]

where $\chi\in\mathcal{C}^{\infty}_{0}(\mathbb{R})$ is a real valued function. Note that in this case, for all $(k,m)\in\mathrm{Supp}\widehat{A}$, we have $\zeta_{k,m}:=\omega_{k,m}^{2}-\sin^{2}\gamma=O({\varepsilon}^{2})$.

The previous analysis shows that we can decompose $\mathcal{W}_{BL}^{0}$ as $\mathcal{W}_{BL,{\varepsilon}^{2}}^{0}+\mathcal{W}_{BL,{\varepsilon}^{3}}^{0}$, where

[TABLE]

In the rest of the paper, we set

[TABLE]

We can then estimate the sizes of $\mathcal{W}^{0}_{\text{inc}},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0}$ and $\mathcal{W}_{BL,{\varepsilon}^{3}}^{0}$ in $L^{\infty}$ and $L^{2}$:

Lemma 2.7**.**

Assume that $\widehat{A}$ is given by (2.16), where $\chi\in\mathcal{C}^{\infty}_{0}(\mathbb{R})$ is a real valued function. We then have the following properties, for all $t\in\mathbb{R}$:

• •

Size of the incident wave packet:

[TABLE]

• •

Size of the ${\varepsilon}^{2}$ boundary layer:

[TABLE]

Note that the normal component of the velocity is much smaller (by a factor $O({\varepsilon}^{2})$).

• •

Size of the ${\varepsilon}^{3}$ boundary layer:

[TABLE]

Note that the normal component of the velocity is much smaller (by a factor $O({\varepsilon}^{3})$).

As a consequence, $\mathcal{W}^{0}$ is of size $O(1)$ in $L^{\infty}$ and $L^{2}$.

Proof.

$\bullet$ We start with the incident wave packet. Note that definition (2.14) allows us to extend the definition of $\mathcal{W}^{0}_{\text{inc}}$ to the whole space $\mathbb{R}^{2}$. According to the Plancherel formula, we have

[TABLE]

by our choice of scaling. Furthermore,

[TABLE]

$\bullet$ We now turn towards the estimates on $\mathcal{W}_{BL,{\varepsilon}^{2}}^{0}$. We start with the $L^{\infty}$ estimate, which is simpler. We recall that $a_{2},a_{3}=O({\varepsilon}^{-2})$, and that $\begin{pmatrix}U_{\lambda_{j}}\\ W_{\lambda_{j}}\\ B_{\lambda_{j}}\end{pmatrix}=O(1)$ for all $k,m$ in the support of $\widehat{A}$. We deduce that

[TABLE]

The $L^{2}$ norm is slightly more involved. Using the Fubini theorem, we have

[TABLE]

Note that if $(k,m)\in\mathrm{Supp}\widehat{A}$, there exist constants $C,c>0$ (independent of ${\varepsilon}$) such that for $j=2,3$,

[TABLE]

Therefore, using the Plancherel theorem,

[TABLE]

Integrating this inequality with respect to $y$, we obtain $\|\mathcal{W}_{BL,{\varepsilon}^{2}}^{0}(t,y)\|_{L^{2}_{x,y}}\lesssim 1$.

$\bullet$ The ${\varepsilon}^{3}$ boundary layer is treated with similar arguments. We recall that $a_{5}=O({\varepsilon}^{-1})$, and that there exists a constant $c$ such that $\lambda_{5}\gtrsim c{\varepsilon}^{-3}$ for all $(k,m)\in\mathrm{Supp}A$. Therefore

[TABLE]

and

[TABLE]

The result follows. ∎

Remark 2.8**.**

• •

The intuition behind the sizes of the ${\varepsilon}^{2}$ and ${\varepsilon}^{3}$ boundary layers is the following. Because of the wave packet, both boundary layers are localized in $x$ in a band of width ${\varepsilon}^{-2}$. As a consequence, the ${\varepsilon}^{2}$ (resp. ${\varepsilon}^{3}$) boundary layer is localized in a region of size ${\varepsilon}^{-2}\times{\varepsilon}^{2}=1$ (resp. ${\varepsilon}^{-2}\times{\varepsilon}^{3}={\varepsilon}$). Hence the $L^{2}$ and $L^{\infty}$ sizes of the ${\varepsilon}^{2}$ boundary layer are of the same order, while the $L^{2}$ size of the ${\varepsilon}^{3}$ boundary layer is

[TABLE]

• •

The same bounds hold for all derivatives with respect to $x$. Concerning the $y$ derivatives, the same bounds hold for the incident wave packet, but a power ${\varepsilon}^{j}$ is lost with each derivation for the ${\varepsilon}^{j}$ boundary layer, $j=2,3$.

3 Weakly nonlinear system

This section aims at providing a description of our approximate solution when the weak nonlinearity comes into play. As anticipated in the introduction, the nonlinear solution can be formally seen as a small perturbation of the linear one.

More precisely, we will start from the (almost) exact solution to the linear system $\mathcal{W}^{0}$, and we will construct an approximate solution to the weakly nonlinear system in the form

[TABLE]

where the corrector $\mathcal{W}^{1}$ compensates the self-interactions of $\mathcal{W}^{0}$. In other words, the nonlinear term will be treated as a source term for the linear equation. One difficulty here is to deal with nonlinear interactions in the context of the wave packet approximation, which is an infinite superposition of plane waves.

In the case of near-critical reflection with the scalings suggested by Dauxois and Young [10], where $\delta$ is the parameter of the weak nonlinearity and the viscosity is of order $O(\varepsilon^{6})$, the weakly nonlinear system is the one considered in (1.9). Let us rewrite this system in the following more compact way. Let $\mathbb{P}$ be the orthogonal projection onto the divergence free vector fields in $L^{2}(\mathbb{R}_{+})^{3}$, i.e. fields of the form $\mathcal{W}=(u,w,b)\in L^{2}(\mathbb{R}^{2}_{+})^{3}$ with $\partial_{x}u+\partial_{y}w=0$ . For $\mathcal{W},\,\mathcal{W}^{\prime}$ also belonging to $H^{1}(\mathbb{R}_{+})^{3}$, define

[TABLE]

and

[TABLE]

Then system (1.9) becomes

[TABLE]

The main result of this section is the following.

Proposition 3.1**.**

Let $\mathcal{W}^{0}$ be given by (2.18). There exists a function $\mathcal{W}^{1}$, which satisfies

[TABLE]

together with the boundary conditions $u^{1}_{|y=0}=w^{1}_{|y=0}=0$, ${\partial}_{y}b^{1}_{|y=0}=0$, where the remainder $r^{1}$ satisfies

[TABLE]

and which can be decomposed as

[TABLE]

The terms of the above decomposition satisfy the following properties:

• •

The ${\varepsilon}^{2}$ boundary layer satisfies

[TABLE]

Moreover, the normal component of the velocity is smaller by a factor $O({\varepsilon}^{2})$. Derivatives with respect to $t$ and $x$ are bounded, while each derivative with respect to $y$ has a cost $O({\varepsilon}^{-2})$.

• •

The ${\varepsilon}^{3}$ boundary layer satisfies

[TABLE]

Moreover, the normal component of the velocity is smaller by a factor $O({\varepsilon}^{3})$. Derivatives with respect to $t$ and $x$ are bounded, while each derivative with respect to $y$ has a cost $O({\varepsilon}^{-3})$.

• •

The second harmonic oscillates around frequencies $\pm 2\omega_{0}$, and satisfies

[TABLE]

• •

The mean flow is much smaller

[TABLE]

for any $s\in\mathbb{N}$. For these second harmonic and mean flow contributions, derivatives with respect to $t,x,y$ are bounded.

Note that the terms do not all have the same size. We need all of them to get a small remainder $r^{1}$ and to lift exactly all boundary conditions. However our convergence result will not be precise enough to give a physical meaning to all correctors (in particular to the mean flow). There are several tools that are important in our construction:

• •

First, since $Q$ is a bilinear operator, we can use definition (2.18) to write $Q(\mathcal{W}^{0},\mathcal{W}^{0})$ as a sum of nine terms. According to Lemma 2.7, the largest of these terms is $Q(\mathcal{W}_{BL,{\varepsilon}^{2}}^{0},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})$. Hence we will need to carefully write the expression of this quadratic term, and in particular, to outline its dependency in time (i.e. its time-frequency) and in space (frequency in the tangential variable, decay in $y$).

• •

By linearity, it is also clear that every term in $Q(\mathcal{W}^{0},\mathcal{W}^{0})$ (except for $Q(\mathcal{W}^{0}_{\text{inc}},\mathcal{W}^{0}_{\text{inc}})$) is a linear superposition of terms of the type

[TABLE]

with $\Re(\lambda)\gg 1$, and $k$ (resp. $\omega$) is in a neighborhood of size ${\varepsilon}^{2}$ of [math] or $2k_{0}$ (resp. [math] or $2{\omega}_{0}$). Therefore it will be sufficient to understand the behavior of solutions of

[TABLE]

with homogeneous boundary conditions. We will perform an almost explicit computation of an approximate solution, relying on an asymptotic expansion of $\mathcal{L}_{\varepsilon}$ and on the boundary layer analysis of the previous section.

Let us now list the nine possible type of interactions by decreasing size in $L^{2}$, and give their typical decay (based on the estimates in Lemma 2.7).

Remark 3.2**.**

Note that boundary layer terms of size ${\varepsilon}^{2}$ (type (a)) will give birth to boundary layer terms of size ${\varepsilon}^{2}$ AND of size ${\varepsilon}^{3}$ (through the boundary condition). Note also that we will not lift terms of type (c) as they are already of size $O({\varepsilon}^{2})$ and can be directly incorporated to the remainder $r^{1}$.

3.1 Interactions of type (a)

In this paragraph, we treat the nonlinear terms which are of size $O(1)$ in $L^{2}$, with an exponential decay of rate $O({\varepsilon}^{-2})$. More precisely, we will prove the following result:

Proposition 3.3**.**

There exists a corrector $\mathcal{W}^{1}_{(a)}$ which satisfies the following properties:

• (i)

$\mathcal{W}^{1}_{(a)}$* is a solution of the evolution equation*

[TABLE]

with

[TABLE]

• (ii)

$\mathcal{W}^{1}_{(a)}$* satisfies exactly the homogeneous boundary conditions (1.10);*

• (iii)

$\mathcal{W}^{1}_{(a)}$* can be decomposed as the sum of a mean flow $\mathcal{W}^{1}_{MF;(a)}$, a second harmonic $\mathcal{W}^{1}_{II;(a)}$, an ${\varepsilon}^{2}$ boundary layer $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(a)}$, and an ${\varepsilon}^{3}$ boundary layer $\mathcal{W}^{1}_{BL,{\varepsilon}^{3};(a)}$*

[TABLE]

• –

The ${\varepsilon}^{2}$ boundary layer satisfies

[TABLE]

Moreover, the normal component of the velocity is smaller by a factor $O({\varepsilon}^{2})$. Derivatives with respect to $t$ and $x$ are bounded, while each derivative with respect to $y$ has a cost $O({\varepsilon}^{-2})$.

• –

The ${\varepsilon}^{3}$ boundary layer satisfies

[TABLE]

Moreover, the normal component of the velocity is smaller by a factor $O({\varepsilon}^{3})$. Derivatives with respect to $t$ and $x$ are bounded, while each derivative with respect to $y$ has a cost $O({\varepsilon}^{-3})$.

• –

The second harmonic satisfies

[TABLE]

• –

The mean flow is much smaller

[TABLE]

For these second harmonic and mean flow contributions, derivatives with respect to $t,x,y$ are bounded.

Our strategy will be the following: we will first find an approximate solution of the equation

[TABLE]

without taking into account the boundary condition. We then lift the remaining trace by using the boundary operator of the previous section. Note that our construction is very reminiscent from the construction of Ekman boundary layers for instance.

We start with a computation of $Q(\mathcal{W}_{BL,{\varepsilon}^{2}}^{0},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})$ and $Q(\mathcal{W}^{0}_{\text{inc}},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})$. Using the formulas of the previous section, we have

[TABLE]

where we wrote in a condensed manner $\widehat{A}$ for $\widehat{A}(k,m)$ defined in (2.16), $a$ for $a_{j}(k,m)$, $(U,W,B)$ for $(U_{\lambda_{j}(k,m)},W_{\lambda_{j}(k,m)},B_{\lambda_{j}(k,m)})$, see (2.17), and used the same convention for $\widehat{A}^{\prime},a^{\prime},(U^{\prime},W^{\prime},B^{\prime})$ (depending on $j^{\prime},k^{\prime},m^{\prime}$). Similarly,

[TABLE]

where $(U_{inc},W_{inc},B_{inc})$ denotes here the Fourier transform of the incident wave, see (2.14).

Note that, in this representation, the Leray projection is simply a multiplier which does not change the general form of the source term. More precisely, since $\mu$ is such that $\Re(\mu)\gtrsim{\varepsilon}^{-2}$, it is given at leading order by

[TABLE]

In particular the second component of $Q(\mathcal{W}_{BL,{\varepsilon}^{2}}^{0},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})+Q(\mathcal{W}^{0}_{\text{inc}},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})$ is smaller than the other components by a factor $O({\varepsilon}^{2})$.

Since $\widehat{A}$ in (2.16) is supported in the vicinity of $\pm(k_{0},m_{0})$, it follows that in the expression above, for $(k,m)\in\mathrm{Supp}\widehat{A}$, $(k^{\prime},m^{\prime})\in\mathrm{Supp}\widehat{A}^{\prime}$, $(k+k^{\prime},\omega+\omega^{\prime})$ is localized either in a ball of size ${\varepsilon}^{2}$ around $(0,0)$, or in a ball of size ${\varepsilon}^{2}$ around $\pm(2k_{0},2\omega_{0})$. Hence, by linear superposition, it is sufficient to consider the equation

[TABLE]

where $\mu$ is such that $\Re(\mu)\gtrsim{\varepsilon}^{-2}$, $U^{\prime},B^{\prime}=O(1)$, $W^{\prime}=O({\varepsilon}^{2})$, $l$ is close to $0,\,\pm 2k_{0}$, $\alpha$ is close to $0,\,\pm 2\omega_{0}$.

Step 1 : approximation of the linear operator. Let us now compute an approximate operator for $\mathcal{L}_{\varepsilon}$ in this scaling. Since the operator $\mathcal{L}_{\varepsilon}$ has constant coefficients, we look for a solution of (3.5) of the form $\mathcal{W}=\delta\exp(ilx-i\alpha t-\mu y)X$, for some vector $X\in\mathbb{C}^{3}$ to be determined. Because of the divergence free condition, we expect the second component of $X$ to be $O({\varepsilon}^{2})$. Therefore, looking at the equation for the second component and taking into account the exponential decay, we infer that the pressure is also $O({\varepsilon}^{2})$, and it is legitimate to keep only the equation on the $u$ and $b$ components of $\mathcal{W}$ (which will be checked a posteriori anyway). The equation then becomes

[TABLE]

where $\mathcal{V}:\mathbb{R}_{+}\times\mathbb{R}^{2}_{+}\to\mathbb{C}^{2}$, $\mathcal{V}^{\prime}=(U^{\prime},B^{\prime})^{T}$, and

[TABLE]

It is well-known that the eigenvalues of $L$ are $\pm i\sin\gamma$, and that the corresponding normalized eigenvectors are $V_{\pm}=\frac{1}{\sqrt{2}}(1,\,\mp i)$. We therefore have that

[TABLE]

with

[TABLE]

We then have the following lemma :

Lemma 3.4**.**

There exists a solution to (3.6) which is given explicitly by

[TABLE]

In particular, it has the same decay and the same size as $\mathcal{V}^{\prime}$.

Proof.

Conjugating the equation (3.6) by $\exp(Lt)$, we get

[TABLE]

Now, by assumption, recalling that $\zeta=\omega_{0}^{2}-\sin^{2}\gamma\approx O(\varepsilon^{2})$ and that $\omega+\omega^{\prime}$ in formulas (3.3)-(3.4) is localized in a neighborhood of $0,\pm 2\omega_{0}$, then

[TABLE]

for ${\varepsilon}$ small enough. Therefore the right-hand side does not generate any secular growth, and we have

[TABLE]

which is the expected formula. ∎

We then define $\mathcal{V}^{1}_{BL,{\varepsilon}^{2};(a1)},\mathcal{V}^{1}_{BL,{\varepsilon}^{2};(a1)}$ by superposition :

[TABLE]

where we recall that in the notation $\mathcal{V}^{1}_{BL,{\varepsilon}^{2};(aj)}$, the subscript $(aj$) refers to the classification of quadratic interactions in Table 2. By Lemma 2.7, it is clear that these quantities are ${\varepsilon}^{2}$ boundary layer terms, whose size in $L^{2}$ and in $L^{\infty}$ are $O(\delta)$.

Step 2 : restoring the divergence-free condition. We then define $\mathcal{W}_{BL,{\varepsilon}^{2};(a1)}^{1}$ (resp. $\mathcal{W}_{BL,{\varepsilon}^{2};(a2)}^{1}$) by adding the normal velocity of $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(a1)}$ (resp. $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(a2)})$ obtained simply by integrating the divergence free condition

[TABLE]

The sizes of $W^{1}_{BL,{\varepsilon}^{2};(a1)},W^{1}_{BL,{\varepsilon}^{2};(a2)}$ in $L^{2}$ and in $L^{\infty}$ are $O(\delta{\varepsilon}^{2})$, since $\lambda,\lambda^{\prime}\approx O(\varepsilon^{-2})$.

Note that $\mathcal{W}_{BL,{\varepsilon}^{2};(a1)}^{1}$ and $\mathcal{W}_{BL,{\varepsilon}^{2};(a2)}^{1}$ are superpositions of vectors of respective form

[TABLE]

By construction, $\sum_{j=1,2}\mathcal{W}_{BL,{\varepsilon}^{2};(aj)}^{1}$ is divergence free, and it can be easily checked that it is a solution of

[TABLE]

where the remainder $r^{1}_{(a),L}$ (coming from the terms which have been neglected in $\mathcal{L}_{\varepsilon}$, i.e. viscous terms, terms involving the normal component $W^{1}_{BL,{\varepsilon}^{2};(aj)}$ and correctors of order $O({\varepsilon}^{2})$ in the Leray projection) satisfies

[TABLE]

Step 3 : Lifting the boundary conditions. Now, $\mathcal{W}_{BL,{\varepsilon}^{2};(a)}^{1}$ has a non-zero trace on the boundary, which must be lifted. In order to do so, we use the boundary operators of the previous section. Let us recall that $\mathcal{W}_{BL,{\varepsilon}^{2};(a)}^{1}$ oscillates around the frequencies [math] and $\pm 2\omega_{0}$, so that we are never in the critical case. More precisely, we will use either the construction in the case $|\zeta|\gtrsim 1$ (this corresponds to the frequencies $\pm 2\omega_{0}$, and will give rise to the second harmonic), or the one in the case $|k|\lesssim\nu^{1/3}={\varepsilon}^{2},|\omega|\lesssim{\varepsilon}^{2}$ (this corresponds to the frequencies close to zero, and will give rise to the mean flow).

Before doing the explicit computation, we need to identify the parts of $\mathcal{W}_{BL,{\varepsilon}^{2};(a)}^{1}$ that will give rise to a second harmonic or to a mean flow. To that end, we write

[TABLE]

where

[TABLE]

If $(k,k^{\prime},m,m^{\prime})\in\mathrm{Supp}\mathcal{A}_{0}$, then $(k+k^{\prime},m+m^{\prime},\omega+\omega^{\prime})=O({\varepsilon}^{2})$, while if $(k,k^{\prime},m,m^{\prime})\in\mathrm{Supp}\mathcal{A}_{II}$, then $(k+k^{\prime},m+m^{\prime},\omega+\omega^{\prime})=\pm 2(k_{0},m_{0},\omega_{0})+O({\varepsilon}^{2})$.

Using the expressions above, we infer that

[TABLE]

where

[TABLE]

We recall that $a,a^{\prime},\lambda,U,U^{\prime}$, etc. are condensed notations for $a_{j},a_{j^{\prime}}$, $\lambda_{j}(k,m,\omega)$, $U_{\lambda},U_{\lambda^{\prime}}$ respectively, see (2.17), while we refer to (2.14) for the notations $U_{inc},W_{inc}$.

We now lift the contributions of $\mathcal{A}_{0}$ and $\mathcal{A}_{II}$ using the results of the preceding section. Note that the time frequency of $\mathcal{A}_{II}$ is $\pm 2\omega_{0}$, therefore it fits the framework of the non-critical case discussed in Lemma 2.3. On the other hand, $\mathcal{A}_{0}$ is a non oscillating term in time and we will apply the results of Lemma 2.4.

$\bullet$ Using Lemma 2.3, we define the boundary layer corrector associated with $\mathcal{A}_{II}$

[TABLE]

and the second harmonic flow associated with the term $(a)$

[TABLE]

where $\Lambda_{2},\Lambda_{3},\Lambda_{5}$ denote the roots with positive real parts of the determinant of the matrix $A_{\nu,\kappa}(\omega+\omega^{\prime},k+k^{\prime},\Lambda_{j})$ defined in (2.5), associated to $\omega+\omega^{\prime},k+k^{\prime}$. The coefficients $\alpha_{2},\alpha_{3},\alpha_{5}$ satisfy (2.8), i.e.

[TABLE]

with $\tilde{\Lambda}_{2}\sim\frac{\Lambda_{2}}{\omega+\omega^{\prime}}=O(1)$, $\Lambda_{3},\Lambda_{5}=O({\varepsilon}^{-3})$ and

[TABLE]

We can prove the following result:

Lemma 3.5**.**

The additional boundary corrector satisfies the following estimate

[TABLE]

If $\sin\gamma>1/2$, the second harmonic flow $\mathcal{W}^{1}_{II;(a)}$ is evanescent far from the boundary.

[TABLE]

If $\sin\gamma<1/2$, the second harmonic flow $\mathcal{W}^{1}_{II;(a)}$ is a wave packet propagating in the outer domain, and satisfying the following estimates

[TABLE]

Proof.

We need first to estimate the coefficients $\alpha_{2},\alpha_{3},\alpha_{5}$, which means that we have to compute a good approximation of $C^{-1}$ for

[TABLE]

The determinant of the matrix $C$ is of the form

[TABLE]

We then obtain the inverse matrix by computing the comatrix. At leading order, we get

[TABLE]

Combining this formula with the estimates on $\mathfrak{u}_{(a)},\mathfrak{w}_{(a)},\mathfrak{b}_{(a)}$, we obtain that

[TABLE]

Integrating with respect to $m,m^{\prime},k-k^{\prime}$ as already done in the proof of Lemma 2.7, we obtain that the Fourier transform with respect to $x$ of $\mathcal{W}_{BL,{\varepsilon}^{3};(aII)}^{1}$ is $O(\delta{\varepsilon}^{-2})$ in $L^{\infty}_{k+k^{\prime},y}$ with decay rate $O({\varepsilon}^{-3})$ in $y$, and supported in a small domain of size $O({\varepsilon}^{2})$ in $k+k^{\prime}$. We then have the following estimates

[TABLE]

Recall that $\Lambda_{2}=O(1)$ has a non negative real part and satisfies approximately

[TABLE]

Since $k+k^{\prime}=\pm 2k_{0}+O({\varepsilon}^{2})$ and $\omega+\omega^{\prime}=\pm 2\omega_{0}+O({\varepsilon}^{2})$, we have that

[TABLE]

More precisely, we can expand $\omega+\omega^{\prime}-2\omega_{0}$, and then $\Lambda_{2}-\Lambda_{0}$ as a (linear) function of $m-m_{0}$ and $k-k_{0}$.

• •

If $\sin\gamma>1/2$ then $\tan\gamma>1/3$ and $\Lambda_{0}$ has a nonnegative real part. Integrating with respect to $m,m^{\prime},k-k^{\prime}$ as done in the proof of Lemma 2.7, we obtain that the Fourier transform with respect to $x$ of $\mathcal{W}^{1}_{II;(a)}$ is $O(\delta)$ in $L^{\infty}_{k+k^{\prime},y}$, and supported in a small domain of size $O({\varepsilon}^{2})$ in $k+k^{\prime}$. We then have the following estimates

[TABLE]

• •

If $\sin\gamma<1/2$, then $\tan\gamma<1/3$ and $\Lambda_{0}$ is purely imaginary. Since there is no decay in $y$ at main order, we have to use the integral with respect to $m,m^{\prime}$ to get some integrability. Going back to the expression of $\det A_{\nu,\kappa,\omega+\omega^{\prime},k+k^{\prime}}$ in Section 2, we find that there exist real constants $\mu_{1},\mu_{2}$ such that

[TABLE]

In other words, up to a change of variable (normal form), $m+m^{\prime}$ is the Fourier variable of $y$. Integrating with respect to $m-m^{\prime},k-k^{\prime}$, we obtain that the Fourier transform with respect to $x,y$ of $\mathcal{W}^{1}_{II;(a)}$ is $O(\delta{\varepsilon}^{-2})$ in $L^{\infty}$, and supported in a small domain of size $O({\varepsilon}^{4})$. We then have the following estimates

[TABLE]

∎

$\bullet$ The contribution of $\mathcal{A}_{0}$ is slightly more complicated, because the boundary operator only lifts two of the components. We first consider the boundary layer part due to $\mathcal{A}_{0}$, namely

[TABLE]

Then as observed in Lemma 2.4, there remains a non-zero trace for the $w$ component, namely

[TABLE]

where $\bar{\alpha}_{3},\bar{\alpha}_{5}$ solve the system

[TABLE]

with $\Lambda_{3},\Lambda_{5}=O(\varepsilon^{-3})$ and $k+k^{\prime},\omega+\omega^{\prime}=O(\varepsilon^{2})$.

In particular, $\bar{\alpha}_{3},\bar{\alpha}_{5}=O({\varepsilon}^{-4})$ and we obtain as previously that

[TABLE]

In addition, we can prove the following result:

Lemma 3.6**.**

The boundary term $\mathfrak{W}^{1}_{(a)}$ can be lifted by a mean flow $\mathcal{W}^{1}_{MF;(a)}$ with divergence free velocity, satisfying the following estimates:

[TABLE]

In particular

[TABLE]

with

[TABLE]

Proof.

Using the estimates on $(\mathfrak{u}_{(a)},\mathfrak{w}_{(a)},\mathfrak{b}_{(a)})$ and the fact that $k+k^{\prime}=O({\varepsilon}^{2})$, we can check that ${\bar{\alpha}_{j}/\Lambda_{j}}=O({\varepsilon}^{-1})$ for $j=3,5$, and $\mathfrak{w}_{(a)}=O(1)$. Integrating with respect to $k,k^{\prime},m,m^{\prime}$, we deduce that the boundary term to be lifted $\mathfrak{W}^{1}_{(a)}$ satisfies

[TABLE]

Furthermore we can build a function $G\in L^{2}\cap L^{\infty}(\mathbb{R})$ such that $\mathfrak{W}^{1}_{(a)}={\partial}_{x}G$, given by

[TABLE]

We then have the following estimates :

[TABLE]

We then define

[TABLE]

for some function $\theta\in\mathcal{C}^{\infty}_{0}(\mathbb{R})$ such that $\theta\sim 1$ in a neighbourhood of zero, with derivative $\theta^{\prime}$. As a consequence, the first two-components of $\mathcal{W}^{1}_{MF;(a)}$ are divergence free, and by definition $\mathcal{W}^{1}_{MF;(a)}$ lifts the remaining trace :

[TABLE]

The explicit formula for $\mathcal{W}^{1}_{MF;(a)}$ leads immediately to

[TABLE]

We conclude the proof by observing that each derivative (with respect to $x,y$) actually improves the bound by a factor $O({\varepsilon}^{2})$. We could also differentiate with respect to time, which would not change the estimates. ∎

Step 4 : consistency of the approximation.

Define

[TABLE]

It satisfies all the properties stated in Proposition 3.3, with the remainder

[TABLE]

Interactions of type (b)

We will now turn to the nonlinear terms which are of size $O({\varepsilon}^{1/2})$ in $L^{2}$, with an exponential decay of rate $O({\varepsilon}^{-3})$. More precisely, we will prove the following result:

Proposition 3.7**.**

There exists a corrector $\mathcal{W}^{1}_{(b)}$ which satisfies the following properties:

• (i)

$\mathcal{W}^{1}_{(b)}$* is a solution of the evolution equation*

[TABLE]

with

[TABLE]

• (ii)

$\mathcal{W}^{1}_{(b)}$* satisfies exactly the homogeneous boundary conditions (1.10);*

• (iii)

$\mathcal{W}^{1}_{(b)}$* can be decomposed as the sum of a mean flow $\mathcal{W}^{1}_{MF;(b)}$, a second harmonic $\mathcal{W}^{1}_{II;(b)}$, and an ${\varepsilon}^{3}$ boundary layer $\mathcal{W}^{1}_{BL,{\varepsilon}^{3};(b)}$*

[TABLE]

• –

The ${\varepsilon}^{3}$ boundary layer satisfies

[TABLE]

Moreover, the normal component of the velocity is smaller by a factor $O({\varepsilon}^{3})$. Derivatives with respect to $t$ and $x$ are bounded, while each derivative with respect to $y$ has a cost $O({\varepsilon}^{-3})$.

• –

The second harmonic satisfies

[TABLE]

• –

The mean flow is much smaller

[TABLE]

For these second harmonic and mean flow contributions, derivatives with respect to $t,x,y$ are bounded.

We start with a computation of $Q(\mathcal{W}_{BL,{\varepsilon}^{2}}^{0},\mathcal{W}_{BL,{\varepsilon}^{3}}^{0})$, $Q(\mathcal{W}^{0}_{\text{inc}},\mathcal{W}_{BL,{\varepsilon}^{3}}^{0})$ and $Q(\mathcal{W}_{BL,{\varepsilon}^{3}}^{0},\mathcal{W}_{BL,{\varepsilon}^{2}}^{0})$. Using the formulas of the previous section and the leading order approximation of the Leray projector in this scaling, we obtain that

[TABLE]

[TABLE]

and

[TABLE]

where we recall that the quadratic term $Q(\mathcal{W},\mathcal{W}^{\prime})$ has been defined in (3.1) and the notations $\hat{A},a_{j},U,W,U_{inc},W_{inc},\mathcal{W}^{0}_{\text{inc}},\mathcal{W}_{BL,\varepsilon^{2}}^{0},\mathcal{W}_{BL,\varepsilon^{3}}^{0}$ have been introduced in (2.14), (2.17), while $U_{5},W_{5}$ are the tangential and the normal component of the velocity field associated to the boundary layer of decay $\exp(-y/\varepsilon^{3})$ in (2.17).

We have thus to study the solutions to the equation

[TABLE]

where $\mu$ is such that $\Re(\mu)\gtrsim{\varepsilon}^{-3}$, $U^{\prime},B^{\prime}=O(1)$, $W^{\prime}=O({\varepsilon}^{3})$, and $l$ (resp. $\alpha$) so located in an ${\varepsilon}^{2}$ neighborhood of [math], $2k_{0}$, or $-2k_{0}$ (resp. of [math], $2\omega_{0}$, or $-2\omega_{0}$).

The strategy will be very similar to the one used in the previous paragraph, with one major difference : since the scaling is different, the linear operator contains viscous dissipation in the $y$ direction at leading order.

Step 1 : approximation of the linear operator. At leading order, the equation then becomes

[TABLE]

where $\mathcal{V}:\mathbb{R}_{+}\times\mathbb{R}^{2}_{+}\to\mathbb{C}^{2}$, and $V^{\prime}=(U^{\prime},B^{\prime})^{T}$.

Lemma 3.8**.**

There exists a matrix $M$ depending on $\alpha,\nu_{0},\kappa_{0}$ and $m=\mu{\varepsilon}^{3}$, but uniformly bounded with respect to all these parameters, such that

[TABLE]

is a solution to (3.12).

Proof.

We look at a solution to (3.12) of the form

[TABLE]

Plugging this Ansatz in the equation, we obtain the following system

[TABLE]

where $m$ is an approximate solution with positive real part to

[TABLE]

which is the main order equation for the eigenvalues of order $\nu^{-\frac{1}{2}}=\varepsilon^{-3}$ in the critical case with small diffusivity, see section 2.

The matrix $M$ is therefore defined by

[TABLE]

We have

[TABLE]

In particular, since $\alpha=\pm 2\omega_{0}+O({\varepsilon}^{2})$ or $\alpha=O({\varepsilon}^{2})$, $\text{det}(M^{-1})$ is bounded away from 0. More precisely, by considering the three different values of $\alpha$, namely $\pm 2\omega_{0}+O({\varepsilon}^{2})$, $O(\varepsilon^{2})$, we can check that $|\text{det}(M^{-1})|\geq C\omega_{0}^{2}$, where the constant value $C$ is independent of $\nu_{0},\kappa_{0}$. We therefore obtain a uniform bound on $M$ in all regimes we will consider. ∎

We then define $\mathcal{V}^{1}_{BL,{\varepsilon}^{3};(b)}$ by superposition :

[TABLE]

[TABLE]

[TABLE]

This quantity is an ${\varepsilon}^{3}$ boundary layer term, whose size is $O(\delta{\varepsilon}^{1/2})$ in $L^{2}$ and $O(\delta)$ in $L^{\infty}$ (as the source term (b)). We indeed recall that $a,a^{\prime}=O({\varepsilon}^{-2})$ and $a_{5},a^{\prime}_{5}=O({\varepsilon}^{-1})$.

Step 2 : restoring the divergence-free condition.

We then define the normal velocity of $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(b1)}$, $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(b2)}$ and $\mathcal{W}^{1}_{BL,{\varepsilon}^{2};(b3)}$ simply by integrating the divergence free condition. Denoting by $M_{1}$ (resp. $M_{2}$) the first (resp. the second) line of the matrix $M$, we have

[TABLE]

[TABLE]

[TABLE]

The size of $W^{1}_{BL,{\varepsilon}^{3};(bj)}$ is $O(\delta{\varepsilon}^{7/2})$ in $L^{2}$ and $O(\delta{\varepsilon}^{3})$ in $L^{\infty}$.

Note that $\mathcal{W}_{BL,{\varepsilon}^{3};(b)}^{1}$ is a superposition of vectors of the form

[TABLE]

By construction, $\sum_{j=1}^{3}\mathcal{W}_{BL,{\varepsilon}^{3};(bj)}^{1}$ is divergence free, and it can be easily checked that it is a solution of

[TABLE]

where the remainder $r^{1}_{(b),L}$ (coming from the terms which have been neglected in $\mathcal{L}_{\varepsilon}$, i.e. viscous terms ${\partial}_{xx}$, terms involving $W^{1}_{BL,{\varepsilon}^{3};(b)}$ and correctors of order $O({\varepsilon}^{3})$ in the Leray projection) satisfies

[TABLE]

Step 3 : Lifting the boundary conditions. Now, $\mathcal{W}_{BL,{\varepsilon}^{3};(bj)}^{1}$ has a non-zero trace on the boundary, which must be lifted. Note that this trace has exactly the same structure as in the previous case, namely

[TABLE]

where

[TABLE]

We now lift the contributions of $\mathcal{A}_{0}$ and $\mathcal{A}_{II}$ exactly as in the previous paragraph. Using Lemma 2.3, we first define the boundary layer corrector associated with $\mathcal{A}_{II}$

[TABLE]

and the second harmonic flow associated with the term $(b)$

[TABLE]

where $\Lambda_{2},\Lambda_{3},\Lambda_{5}$ denote the roots with positive real parts of the determinant of the matrix $A_{\nu,\kappa}(\omega+\omega^{\prime},k+k^{\prime},\Lambda_{j})$ defined in (2.5), associated to $\omega+\omega^{\prime},k+k^{\prime}$. The coefficients $\beta_{j}$ are defined by

[TABLE]

where the matrix $C$ is defined in (3.8), and satisfies

[TABLE]

As a consequence,

[TABLE]

Note that $\beta_{3}$ and $\beta_{5}$ are of the same size as $\alpha_{3}$ and $\alpha_{5}$. However, $\beta_{2}=O({\varepsilon}^{-1})$ is a priori smaller than $\alpha_{2}=O({\varepsilon}^{-2}).$ As a consequence, the estimate on $\mathcal{W}_{BL,{\varepsilon}^{3};(bII)}^{1}$ is exactly the same as the one on $\mathcal{W}_{BL,{\varepsilon}^{3};(aII)}^{1}$, but the estimate on $\mathcal{W}^{1}_{II;(b)}$ is smaller than the one on $\mathcal{W}^{1}_{II;(a)}$ by a factor ${\varepsilon}$.

We then define the boundary layer part due to $\mathcal{A}_{0}$, namely

[TABLE]

We finally lift the remaining trace on the $w$ component by the mean flow term

[TABLE]

where

[TABLE]

and

[TABLE]

denoting

[TABLE]

We then have exactly the same estimates as for the interactions of type (a).

Step 4 : consistency of the approximation.

Define

[TABLE]

It satisfies all the properties stated in Proposition 3.7, with the remainder

[TABLE]

4 Accuracy of the approximation

The aim of this section is to quantify the accuracy of our construction, by providing an $L^{2}$ estimate of the difference between the approximate solution and the solution to the Cauchy problem for (3.2)-(1.10) with initial data

[TABLE]

where $\mathcal{W}^{0},\mathcal{W}^{1}$ have been defined in (2.18) and Proposition 3.1 respectively.

4.1 The error estimate

The first step is to estimate the size of the error in the approximation of the evolution equation, coming from terms of different types  :

$\bullet$ the viscosity for the incident wave packet $\mathcal{W}^{0}_{\text{inc}}$;

$\bullet$ the remainder $r^{1}$ in Proposition 3.1 (corresponding to the error coming from $\mathcal{W}^{1}$);

$\bullet$ the quadratic interactions (c) in Table 2 that have not been corrected;

$\bullet$ the interactions between $\mathcal{W}^{0}$ and $\mathcal{W}^{1}$ and the self-interactions of $\mathcal{W}^{1}$.

Proposition 4.1**.**

Consider $\mathcal{W}_{app}=\mathcal{W}^{0}+\mathcal{W}^{1},$ with $\mathcal{W}^{0},\mathcal{W}^{1}$ defined in (2.18) and Proposition 3.1 respectively. Then $\mathcal{W}_{app}$ solves

[TABLE]

and satisfies exactly the boundary conditions (1.10), where

[TABLE]

Moreover,

• •

$\|\mathcal{W}_{app}\|_{L^{2}(\mathbb{R}_{+}^{2})}=O(1),\quad\|\mathcal{W}_{app}\|_{L^{\infty}(\mathbb{R}^{2}_{+})}=O(1);$

• •

$\|\nabla\mathcal{W}_{app}\|_{L^{2}(\mathbb{R}_{+}^{2})}=O(\varepsilon^{-2}),\quad\|\nabla\mathcal{W}_{app}\|_{L^{\infty}(\mathbb{R}^{2}_{+})}=O(\varepsilon^{-2});$

• •

$\|R_{app}\|_{L^{2}(\mathbb{R}_{+}^{2})}=O(\delta\varepsilon^{2})+O(\delta^{2})+O(\varepsilon^{6})$.

Proof.

The consistency of $\mathcal{W}_{app}$ can be easily checked by recalling the equations satisfied by $\mathcal{W}^{0},\mathcal{W}^{1}$ in Lemma 2.7 and Proposition 3.1 respectively. Again, the sizes of $\mathcal{W}_{app}$ and $\nabla\mathcal{W}_{app}$ are a direct consequence of Lemma 2.7 and Proposition 3.1. Let us estimate the size of the remainder.

$\bullet$ From Proposition 3.1, $\|r^{1}\|_{L^{2}(\mathbb{R}^{2}_{+})}=O(\delta\varepsilon^{2})$;

$\bullet$ From Table 2, $\|\delta\times(c)\|_{L^{2}(\mathbb{R}^{2}_{+})}=O(\delta\varepsilon^{2})$;

$\bullet$ From Lemma 2.7, the size of the remainder due to the viscosity applied to the incident wave packet is $O(\varepsilon^{6})$;

$\bullet$ It remains then to estimate the interactions between $\mathcal{W}_{app}$ and $\mathcal{W}^{1}$:

[TABLE]

[TABLE]

[TABLE]

∎

Remark 4.2**.**

It is natural to wonder whether we could iterate our construction in order to define an approximate solution up to any order. In other words, given $N\in\mathbb{N}$ arbitrary, is it possible to define a solution $\mathcal{W}_{app}$ (different from the solution defined above) such that

[TABLE]

There are several limitations preventing us to do so:

• •

First, because of the degeneracy of the root $\lambda_{2}(k+k^{\prime},\omega+\omega^{\prime})$ in the region $(k+k^{\prime},m+m^{\prime})\in\mathrm{Supp}\mathcal{A}_{0}$, we have chosen to lift the traces $\mathfrak{W}^{1}_{(a)}$ and $\mathfrak{W}^{1}_{(b)}$ by correctors that have no physical relevance, namely $\mathcal{W}^{1}_{MF;(a)}$ and $\mathcal{W}^{1}_{MF;(b)}$. In order to reach a better approximation, we need to understand how to lift these traces in a more specific way, which possibly entails tracking down the exact degeneracy of $\lambda_{2}$.

• •

In a similar fashion, in the above construction, we have treated $(I-\mathbb{P})((u^{0}{\partial}_{x}+w^{0}{\partial}_{y})\mathcal{W}^{0})$ as a remainder. If we want to go further in the definition of the approximate solution, we need to understand the precise structure of this remainder and to correct it by additional terms.

• •

The next step in the construction of an approximate solution would be to consider the remainder terms stemming from the interactions between $\mathcal{W}^{0}$ and $\mathcal{W}^{1}$, namely $\delta Q(\mathcal{W}^{0},\mathcal{W}^{1})$ and $\delta Q(\mathcal{W}^{1},\mathcal{W}^{0})$. One could try and apply the same method as above in order to add a corrector $\mathcal{W}^{2}$, satisfying

[TABLE]

with $\|r^{2}\|_{L^{2}}\ll\delta^{2}$.

However, there is a major difference with the case of the self-interactions of $\mathcal{W}^{0}$ (i.e. the term $Q(\mathcal{W}^{0},\mathcal{W}^{0})$). Indeed, the term $\mathcal{W}^{0}$ oscillates in time around the frequencies $\pm\omega_{0}$; whence $Q(\mathcal{W}^{0},\mathcal{W}^{0})$ oscillates around the frequencies $\pm 2\omega_{0}$ and [math]. Since the resonance frequencies of the rotation operator $L$ are $\pm\omega_{0}$, the definition of $\mathcal{W}^{1}$ does not create any resonance, see Lemma 3.4. But the terms $Q(\mathcal{W}^{0},\mathcal{W}^{1})$ and $Q(\mathcal{W}^{1},\mathcal{W}^{0})$ contain the frequencies $\pm 3\omega_{0},\pm\omega_{0}$, and therefore resonances may occur.

A more careful look at the quadratic terms shows that in fact, the resonant interactions between $\mathcal{W}^{0}_{\text{inc}}+\mathcal{W}_{BL,{\varepsilon}^{2}}^{0}$ and $\mathcal{W}_{BL,{\varepsilon}^{2}}^{1}$ cancel at main order. This is due to the fact that the eigenvectors $(U_{inc},B_{inc})$ and $(U_{\lambda},B_{\lambda})$ are both $(1,\pm i)+O({\varepsilon}^{2})$ when $|\lambda|\propto{\varepsilon}^{-2}$. In more abstract terms, this cancellation comes from the Jacobi identity. However, this algebraic cancellation no longer holds for the terms involving the ${\varepsilon}^{3}$ boundary layer. As a consequence, we cannot lift these interactions with the tools developed above, and new ideas must be introduced if one wants to push the iteration further.

4.2 The stability inequality

Here we will establish the stability inequality (1.12) leading to Theorem 1.1.

Recalling the definitions of $\mathcal{L}_{\varepsilon}$ and $Q(\mathcal{W},\mathcal{W}^{\prime})$ in (3.2), the approximate solution satisfies the following equation

[TABLE]

where

[TABLE]

We point out that the strong formulation of the system, which is required in order to get the energy estimate leading to the stability inequality, is actually satisfied by a sequence approximate solutions, which are smooth by Friedrichs approximation, see the Appendix. Then the energy inequality for the weak solution $\mathcal{W}$ is obtained by passing to the limit. With a slight abuse of notation, here we omit this passage and we consider directly the equation for $\mathcal{W}$.

We write the equation for the difference $\mathcal{W}_{app}-\mathcal{W}$ and take the scalar product against $\mathcal{W}_{app}-\mathcal{W}$. This yields:

[TABLE]

We analyse each term separately.

[TABLE]

since $\mathbb{P}$ is symmetric by definition and $\mathbf{L}$ is skew-symmetric.

[TABLE]

where $\mathbb{P}$ is bounded in $L^{2}$ by definition and the last inequality follows from Proposition 4.1.

$\bullet$ Because of the symmetry of $\mathbb{P}$, by integrating by parts the third integral reads

[TABLE]

which is zero thanks to the zero trace of the velocity field and the divergence free condition.

$\bullet$ Integrating by part the term with the Laplacian, one gets

[TABLE]

where the last equality follows from the boundary conditions (1.10).

$\bullet$ It remains to deal with the remainder

[TABLE]

The Gronwall inequality yields the result.

Remark 4.3**.**

We point out that there is an alternative way to establish the stability inequality, by estimating the last term in the following way:

[TABLE]

This yields

[TABLE]

*Note that this last version of the stability estimate allows us to have a stable approximate solution for longer times, of order $\varepsilon^{-1}$, provided that $\delta\leq\varepsilon^{3}$, with a remainder of size $\delta^{\frac{1}{2}}\varepsilon^{3}$ ($\delta^{\frac{1}{2}}\varepsilon^{2}$ in case of secular growths).

However, in the present paper we consider interval of times $[0,t]$ with $t=O(1)$, then we rather rely on the first version of the stability estimate (1.13), which provides a smaller remainder for $\delta<\varepsilon^{2}$.*

Appendix

We present a general approach to proving the global existence of weak solutions to system (1.9)-(1.10) in $\mathbb{R}^{2}_{+}=\mathbb{R}\times\mathbb{R}^{+}$. We adapt the result of existence of global weak (Leray) solutions for the incompressible Navier-Stokes equations in [5], in the case of a general domain of $\mathbb{R}^{d},\,d=2,3$. Since there are only minor modifications, we shortly sketch the argument and we refer to the book [5] for a complete proof.

We point out that another more explicit approach to prove the global existence existence of (local energy) weak solutions to the 3D incompressible Navier-Stokes equations has recently been developed in [15].

We start by recalling the definition of weak solutions to the Cauchy problem associated with system (1.9)-(1.10). Some of the notations introduced below have been defined in Section 1.4.

[TABLE]

Definition 4.4**.**

A weak solution to the Boussinesq system (1.9)-(1.10) on $\mathbb{R}^{+}\times\mathbb{R}^{2}_{+}$ with initial data $\mathcal{W}_{0}\in\mathbb{K}$ is

[TABLE]

such that, for any function $\phi\in C^{1}(\mathbb{R}^{+};\mathbb{V}_{\sigma})$,

[TABLE]

Moreover, it satisfies the energy inequality in (1.12).

We now provide a sketch of the proof of Proposition 1.1.

There are essentially two main difficulties in dealing with the half plane $\mathbb{R}^{2}_{+}$.

$\bullet$ The domain is unbounded: this prevents us from using some compactness results for bounded domains;

$\bullet$ The construction of an approximation for the Leray projection (in other words, the construction of an approximation with divergence free velocity field at every step) requires some work to do. In the case of bounded domains, an Hilbertian basis of $\mathbb{K}$ is indeed provided by the spectral theorem for self-adjoint compact operators in Hilbert spaces, since the inverse of the Stokes operator is compact, see [5]. In the whole space, the Leray projector has an explicit expression in terms of the symbol of a pseudodifferential operator homogeneous of degree 0, see [2]. In the case of the half space, the embedding of $H_{0}^{1}\subset L^{2}$ is only continuous and we cannot use the Fourier transform in the vertical direction, then a more implicit approach is needed.

The strategy presented in [Chapther 2, [5]] consists in the following main steps.

Construct a family of orthogonal projectors $\mathbb{P}_{\eta},\;\eta\in\mathbb{R},$ on $\mathbb{K}$. This is done as in [[5], Chapter 2], by proving that the inverse of the Stokes operator

[TABLE]

is continuous, self-adjoint, one-to-one, contractive and with range $R(B)$ dense in $\mathbb{K}$. This allows us to define the inverse of $B$ as an unbounded operator $A$ with a dense domain of definition. Since $A$ is self-adjoint, the Stokes operator $A-Id$ is also self-adjoint and then one can apply the spectral theorem to get the following result.

Proposition 4.5**.**

There exists a family of orthogonal projectors on $\mathbb{K}$, denoted by $\mathbb{P}_{\eta}$, with $\eta\in\mathbb{R}$, which commutes with the Stokes operator $A-Id$ and satisfies the following properties.

• •

$\mathbb{P}_{\eta}\mathbb{P}_{\eta^{\prime}}=\mathbb{P}_{\inf(\eta,\eta^{\prime})}$* for any $(\eta,\eta^{\prime})\in\mathbb{R}^{2}$;*

• •

$\mathbb{P}_{\eta}=0$* for $\eta<0$ and for any $\mathcal{W}\in\mathbb{K}$*

[TABLE]

• •

The family $\mathbb{P}_{\eta}$ satisfies the right continuity, then for every $\mathcal{W}\in\mathbb{K}$

[TABLE]

• •

For any $\mathcal{W}\in\mathbb{K}$, the function $\eta\rightarrow(\mathbb{P}_{\eta}\mathcal{W},\mathcal{W})_{\mathbb{K}}=\|\mathbb{P}_{\eta}\mathcal{W}\|_{\mathbb{K}}^{2}$ is increasing, and

[TABLE]

• •

$\mathbb{P}_{\eta}\mathcal{W}\in\mathcal{V}_{\sigma}$* for any $\mathcal{W}\in\mathbb{K}$ and*

[TABLE] 2. 2.

Show that the family of smoothing operators $\mathbb{P}_{\eta}$ can be extended to the dual space $\mathbb{V}^{\prime}$, see [Corollary 2.2, [5]] and converges to the Leray projector $\mathbb{P}$ in the sense of the $\mathbb{K}$ norm: for any $\mathcal{W}\in\mathbb{K}$,

[TABLE]

as in [Proposition 2.4, [5]]. 3. 3.

Denoting by $\mathbb{K}_{k}$ the space $\mathbb{P}_{k}\mathbb{K}$ for any integer $k$, define the approximate system of ODEs (Friedrichs approximation)

[TABLE]

where $\mathcal{W}_{k}(t)\in\mathbb{K}_{k}$ for every $k$.

Then the rest of the proof and the uniqueness result follow in the usual way, as for the case of Leray solutions to the 2D incompressible Navier-Stokes equations in bounded domains, see [5] for further details.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program Grant agreement No 637653, project BLOC “Mathematical Study of Boundary Layers in Oceanic Motion”. This work was supported by the SingFlows project, grant ANR-18-CE40-0027 of the French National Research Agency (ANR). The authors would like to thank Thierry Dauxois for several helpful discussions.