Long time behavior of the two-dimensional Boussinesq equations without buoyancy diffusion

Accepted Manuscript Long time behavior of the two-dimensional Boussinesq equations without buoyancy diffusion Charles R. Doering, Jiahong Wu, Kun Zhao, Xiaoming Zheng

PII: DOI: Reference:

S0167-2789(17)30362-7 https://doi.org/10.1016/j.physd.2017.12.013 PHYSD 31994

To appear in:

Physica D

Please cite this article as: C.R. Doering, J. Wu, K. Zhao, X. Zheng, Long time behavior of the two-dimensional Boussinesq equations without buoyancy diffusion, Physica D (2018), https://doi.org/10.1016/j.physd.2017.12.013 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG A BSTRACT. We study the global well-posedness and stability / instability of perturbations near a special type of hydrostatic equilibrium associated with the 2D Boussinesq equations without buoyancy (a.k.a. thermal) diffusion on a bounded domain subject to stress-free boundary conditions. The boundary of the domain is not necessarily smooth and may have corners such as in the case of rectangles. We achieve three goals. First, we establish the global-in-time existence and uniqueness of large-amplitude classical solutions. Efforts are made to reduce the regularity assumptions on the initial data. Second, we obtain the large-time asymptotics of the full nonlinear perturbation. In particular, we show that the kinetic energy and the first order derivatives of the velocity field converge to zero as time goes to infinity, regardless of the magnitude of the initial data, and the flow stratifies in the vertical direction in a weak topology. Third, we prove the linear stability of the hydrostatic equilibrium T (y) satisfying T 0 (y) = α > 0, and the linear instability of periodic perturbations when T 0 (y) = α < 0. Numerical simulations are supplemented to corroborate the analytical results and predict some phenomena that are not proved. The authors are pleased to dedicate this paper to Professor Edriss Saleh Titi on the occasion of his 60th birthday. Professor Titi’s myriad research contributions—including contributions to the problem constituting the focus of this work—and leadership in the mathematical fluid dynamics community serve as an inspiration to his students, collaborators and colleagues.

1. I NTRODUCTION 1.1. Background. The two-dimensional (2D) incompressible Boussinesq equations describing the motion of buoyancy driven fluid flows are (1.1)

 ∂ u + u · ∇u + ∇P = νx ∂xx u + νy ∂yy u + θe2 ,   t ∂t θ + u · ∇θ = κx ∂xx θ + κy ∂yy θ,   ∇ · u = 0,

where x = (x, y) ∈ R2 , t > 0 and u = (u, v)T , P, θ denote, respectively, the velocity, pressure and buoyancy fields. The vertical unit vector is e2 = (0, 1)T and νx , νy , κx , κy are the (real non-negative) viscosity and buoyancy diffusion coefficients. The acceleration term θe2 in the momentum equation, the

first equation in (1.1), models the buoyancy due to density variations in presence of gravity. Such density variations are often due to thermal expansion so we will interchangeably refer to θ as a temperature field and its evolution, the second equation in (1.1) as the temperature equation. This system is routinely used to model flows across a tremendous range of length and time scales from microfluidics and biophysics to geodynamics and astrophysics. It plays an important role in the study of turbulence and, suitably modified, in atmospheric and oceanographic situations where rotation and stratification play a dominant role [46, 47, 48]. It is also known for its close connection with fundamental models in mathematical fluid mechanics: 2010 Mathematics Subject Classification. 35Q35, 35B40. Key words and phrases. 2D Boussinesq equations, initial-boundary value problem, classical solution, long-time behavior. 1

2

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

• Analogy to 3D Euler Equations. One important feature of the non-dissipative (νx = νy = κx = κy = 0) 2D Boussinesq system is its analogy to the 3D incompressible Euler equations. Indeed, the vorticity formulation of the 3D incompressible Euler equations for axisymmetric swirling flow in cylindrical coordinates, (r, φ, x3 ), read   φ  φ  φ   1 ω ω  r 3 ∂t ω + v ∂r + v ∂x3 = − 4 ∂x3 (rv φ )2 , r r r r  ∂ (rv φ )2  + v r ∂ (rv φ )2  + v 3 ∂ (rv φ )2  = 0, t

r

x3

where ω φ denotes the angular vorticity in cylindrical coordinates. On the other hand, the vorticity formulation of the 2D non-dissipative Boussinesq system in Cartesian coordinates takes the form ( ∂t ω + v∂y ω + u∂x ω = ∂x θ, ∂t θ + v∂y θ + u∂x θ = 0,

where ω = ∂x v − ∂y u denotes the 2D vorticity in Cartesian coordinates. In view of the two systems we see the correspondence: ωφ ←→ ω, r

v r ←→ v,

v 3 ←→ u,

r ←→ y,

x3 ←→ x,

(rv φ )2 ←→ −θ.

Therefore, away from the line of singularity (r = 0), one should expect that the qualitative behavior of the 3D incompressible Euler equations for axisymmetric swirling flow is identical to that of the 2D non-dissipative Boussinesq equations. • Vortex Stretching Effect . A second important feature of the 2D Boussinesq system is that it shares a similar vortex stretching effect as that in 3D flows. As a matter of fact, by letting W = ∇θ = (∂x θ, ∂y θ)T , one can show that for classical solutions, (1.1) with νx = νy = κx = κy = 0 is equivalent to  ∂ ω + u · ∇ω = ∂x θ,   t ∇ · u = 0,   ∂t W + u · ∇W = −(∇uT )W,

where (∇uT ) = (∇u ∇v) is a 2 × 2 matrix. The first equation indicates that qualitatively the growth of the vorticity depends on the temporal accumulation of ∂x θ which is the first component of W. From the third equation we see that W has the same degree of regularity as ∇u. Therefore, qualitatively the growth of the vorticity depends on the temporal accumulation of ∇u – a scenario that is similar to the vortex stretching effect in the 3D incompressible Euler equations. In fact, in view of the vorticity equation of the 3D incompressible Euler equations: ∂t Σ + U · ∇Σ = Σ · ∇U,

where U is the 3D velocity and Σ = ∇ × U, we see that the growth of the vorticity depends on the temporal accumulation of the gradient of the velocity field, as well. 1.2. Literature Review. Due to its physical background and mathematical features, the 2D Boussinesq system has attracted considerable attention from the community of mathematical fluid mechanics in recent years. The following results concerning the qualitative behavior (global/local well-posedness, blowup criteria, regularity, explicit solutions, finite-time singularities, et al) of (1.1) have been well documented in the literature: • global well-posedness, blowup criteria and regularity when νx , νy , κx , κy are not all zeros [1, 2, 3, 4, 11, 12, 13, 18, 20, 21, 22, 27, 28, 29, 30, 32, 34, 35, 40, 41, 42, 45, 64], • local well-posedness, blowup criteria, explicit solutions and finite-time singularities when all four parameters are zeros [5, 14, 15, 16, 17, 19, 23, 33, 50, 54],

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION

3

• well-posedness and regularity with critical and supercritical dissipation [37, 38, 43, 52, 60, 61, 62]. Comparing with the magnitude of research conducted on the well-posedness of the model, the longtime asymptotic behavior of (1.1) has been investigated relatively little. To the authors’ knowledge, the following results have been achieved so far: • exponential decay of θ to constant states and uniform-in-time boundedness of kinetic energy for initial-boundary value problems on 2D bounded smooth domains for large data when νx = νy = 0, κx = κy = κ > 0 [64], • uniform-in-time boundedness of kinetic energy for initial-boundary value problems on 2D bounded smooth domains for large data when νx = νy = ν > 0, κx = κy = 0 [40], • algebraic decay of (u, θ) to constant ground states for the Cauchy problem on R3 for small data when νx = νy = ν > 0, κx = κy = κ > 0 [9]. • long time averaged heat transport sustained by thermal boundary conditions, i.e., bounds for Rayleigh-B´enard convection [59]. • existence of a global attractor containing infinitely many invariant manifolds on 2D periodic domains when νx = νy = ν > 0, κx = κy = 0 [7]. 1.3. Motivation and Goal. The 2D Boussinesq system describes the motion of buoyancy driven flows. The hydrostatic equilibria associated with buoyancy driven flows take the form [0, Phe (y), θhe (y)] where the vertically stratified pressure and temperature satisfy ∂y Phe (y) = θhe (y) obtained from the first equation of (1.1) by setting u = 0. Moreover, when the thermal diffusion coefficients κx , κy > 0, by formally substituting the ansatz into the second equation of (1.1), we see that ∂yy θhe (y) = 0 which suggests that the equilibrium temperature is an affine function of the depth. Mathematically, it is then natural to ask whether such kind of ansatz is stable under appropriate initial and boundary conditions. When the dissipation parameters satisfy νx = νy = ν > 0 and κx = κy = κ > 0, it is not difficult to check that when (1.1) is set on a bounded domain with sufficiently smooth boundary and subject to the boundary conditions: θ|∂Ω = αy + θ¯ & u|∂Ω = 0 (no-flow), or θ|∂Ω = αy + θ¯ & u · n|∂Ω = 0, ω|∂Ω = 0 ¯ is globally (stress-free), where ω = ∂x v − ∂y u denotes the 2D vorticity, the ansatz (u∞ , θ∞ ) = (0, αy + θ) asymptotically stable, provided α > 0. Indeed, the global asymptotic stability follows from the energy conservation law: (1.2)


d ¯ 2 2 = 0, ¯ 2 2 + ανkωk2 2 + κk∇(θ − αy − θ)k αkuk2L2 + kθ − αy − θk L L L dt

and the Poincar´e inequality. On the other hand, in certain circumstances when the thermal diffusion is insignificant, (1.1) becomes (1.3)

 ∂ u + u · ∇u + ∇P = νx ∂xx u + νy ∂yy u + θe2 ,   t ∂t θ + u · ∇θ = 0,   ∇ · u = 0,

which arises naturally as a relevant system in geophysics [31, 49, 56]. Comparing with (1.1), the system (1.3) is only partially dissipative with missing dissipation in the second major function of interest. Such a degeneracy, together with the nonlinear coupling through advection and gravity, makes the mathematical analysis of (1.3) considerably more challenging than that for (1.1). Although the global well-posedness of (1.3) has been proved under various initial and/or boundary conditions, the long-time behavior is still elusive. So far the only such results can be found in [7] and [40], in which it was shown that the kinetic energy associated with (1.3) on bounded domains is uniformly bounded with respect to time, which indicates that the flow is trapped inside of a finite ball in the energy space. Moreover, it was shown that (1.3) has a so-called weak sigma-attractor which contains infinitely many invariant manifolds in

4

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

which several universal properties of the Batchelor, Kraichnan, Leith theory of turbulence are potentially present. However, we note that the results in [7, 40] do not imply the unconditional “slowing down” of the flow, and therefore provide no information about the global stability of hydrostatic equilibria associated with (1.3). The objective of this paper is to study the global well-posedness of (1.3) on 2D bounded domains with non-smooth boundaries for large data, and investigate the stability/instability of hydrostatic equilibria associated with the partially-dissipative system. The spatial domain Ω ⊂ R2 under consideration is either a rectangle or more generally a connected bounded domain with Lipschitz boundary satisfying some minor constraints that we are specifying now. Since ∂Ω concerned here is not smooth and may have corners such as in the case of rectangles, we impose additional conditions on Ω (in addition to ∂Ω being Lipschitz) so that some of the classical analysis tools such as the Calder´on-Zygmund inequality are still valid for such domains. Jerison and Kenig in [36] give a striking example in which the domain Ω has C 1 boundary and the Calder´on-Zygmund inequality fails. We focus on bounded and simply connected planar domains with ∂Ω belonging to C 2,γ (γ ∈ (0, 1)) except at a finite number of points, where ∂Ω is a corner of angle in (0, π2 ]. Planar polygonal domains such as rectangles in R2 are special examples of such domains. As shown in [6, 25, 39], the Calder´on-Zygmund inequality and some other tools still hold for such domains. This paper attempts to achieve three main goals. The first is to establish the global-in-time existence and uniqueness of classical solutions to the following initial-boundary value problem:  ∂t u + u · ∇u + ∇P = ν∆u + θe2 ,       ∂ θ + u · ∇θ = 0,   t ∇ · u = 0,     (u, θ)(x, 0) = (u0 , θ0 )(x),    u · n| = 0, ω| = 0, ∂Ω ∂Ω

(1.4)

where Ω ⊂ R2 is as specified in the previous paragraph, n is the unit outward normal to ∂Ω, and ω = ∂x v−∂y u is the 2D vorticity. A particular type of such domains are the rectangles. (1.4) in the whole space R2 or in a smooth bounded domain with the Dirichlet boundary condition has been studied previously and important results on the well-posedness have been obtained (see, e.g., [7, 13, 32, ?, 41, 40]). In our study presented in this paper, we attempt to make weak regularity assumptions on the initial data. In addition, the boundary of the domain here is not smooth. Motivated by the reasoning for the fully dissipative system (see (1.2) and the paragraph therein), the second goal is to understand the large time behavior and the global stability of the perturbation near the hydrostatic equilibrium [uhe , Phe (y), θhe (y)] given by uhe = 0,

¯ θhe (y) = T (y) ≡ αy + θ,

The perturbation (1.5)

satisfies

(1.6)

e = u − uhe , u

Pe = P − Phe (y),

Phe =

Z

0

y

T (z) dz.

θe = θ − θhe (y)

 e 2,  e+u e · ∇e ∂t u u + ∇Pe = ν∆e u + θe     ∂ θe + u  e · ∇θe + ve T 0 (y) = 0,   t e = 0; ∇·u     e  (e u, θ)(x, 0) = (e u0 , θe0 )(x),    u e · n|∂Ω = 0, ω e |∂Ω = 0,

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION

5

e denotes the corresponding vorticity. where ω e =∇×u The third goal is to prove the linear stability of the hydrostatic equilibrium when the vertically stratified temperature T (y) satisfies T 0 (y) = α > 0 and the linear instability when T 0 (y) = α < 0. In order to achieve this goal, we first derive an equivalent system of equations for the perturbation near the hydrostatic equilibrium (see (4.7) below) and the corresponding linearized system (4.8).

1.4. Results. We now state our main results. Before doing so, we introduce some notations for convenience. Notation 1.1. Throughout this paper, k · kLp ,k · kL∞ and k · kW s,p denote the norms of the usual Lebesgue measurable spaces Lp (Ω), L∞ (Ω) and the usual Sobolev space W s,p (Ω), respectively. For p = 2, we denote the norms k · kL2 and k · kW s,2 by k · k and k · kH s , respectively. Unless specified, ci denote generic constants which are independent of the unknown functions and t, but may depend on ν, Ω and initial data. The values of the constants may vary line by line according to the context. We present three main results. The first one assesses the global-in-time existence and uniqueness of solutions to the initial-boundary value problem (1.4). The second one describes the large time behavior of the nonlinear perturbation while the third establishes the linear stability when the vertically stratified temperature T (y) satisfies T 0 (y) = α > 0 and the linear instability for periodic perturbations when T 0 (y) = α < 0. Theorem 1.1. Let Ω ⊂ R2 be a rectangle or a more general bounded and simply connected planar domain with ∂Ω belonging to C 2,γ (γ ∈ (0, 1)) except at a finite number of points, where ∂Ω is a corner of angle in (0, π2 ]. Consider the initial-boundary value problem (1.4). Suppose that the initial data u0 ∈ H 2 (Ω) and θ0 ∈ H 1 (Ω), ∇ · u0 = 0 and are compatible with the boundary conditions. Then (1.4) has a unique global-in-time solution (u, θ) satisfying (1.7)

u ∈ L∞ ((0, t); H 2 (Ω)) ∩ L2 ((0, t); H 3 (Ω))

and

θ ∈ L∞ ((0, t); H 1 (Ω))

for any t > 0. In addition, if (u0 , θ0 ) ∈ H 3 (Ω), then u ∈ L∞ ((0, t); H 3 (Ω)) ∩ L2 ((0, t); H 4 (Ω)) and θ ∈ L∞ ((0, t); H 3 (Ω)) for any t > 0. Our second main result focuses on the large time behavior of solutions to (1.6), the full nonlinear e. equations for the perturbation [e u, Pe, θ]

Theorem 1.2. Assume Ω ⊂ R2 is either a rectangle or a more general Lipschitz domain as specified e . Assume the previously. Consider the initial-boundary value problem (1.6) for the perturbation [e u, Pe, θ] 2 1 ∞ e 0 ∈ H (Ω), θe0 ∈ H (Ω) ∩ L (Ω), ∇ · u e 0 = 0, and are compatible with the boundary initial perturbations u e satisfying conditions. Then (1.6) possesses a unique and global-in-time solution (e u, θ) e ∈ L∞ ((0, t); H 2 (Ω)) ∩ L2 ((0, t); H 3 (Ω)), u

θe ∈ L∞ ((0, t); H 1 (Ω) ∩ L∞ (Ω)).

e obeys the following large-time behavior, as t → ∞, Furthermore, if T 0 (y) = α > 0, the solution (e u, θ) • ke u(t)kL2 → 0,

(1.8)

e • kθ(t)k → c0 ,

e (t)k → 0, k∇e u(t)kL2 → 0, k∂t u q u0 k2 + kθe0 k2 , 0 < c0 < αke

e 2 kH −1 → 0, • k∇Pe − θe

and there exists a positive constant C , which is independent of t > 0, such that (1.9)

ke u(t)kH 2 ≤ C,

∀ t > 0.

Remark 1.1. We remark that by combining the uniform estimate (1.9) with the arguments in [7], one can establish the existence of global attractors in stronger topological spaces, as was pointed out in [7].

6

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

The third main result assesses the linear stability when the vertically stratified temperature T (y) satisfies T 0 (y) = α > 0 and the linear instability for periodic perturbations when T 0 (y) = α < 0. As a preparation, we derive an equivalent system of equations for the perturbation near the hydrostatic equilibrium and the corresponding linearized system. In contrast to the periodic case or the whole space case, this process is more involved due to the nonsmoothness of the domain. We invoke a result asserting the existence and boundedness of the operator projecting onto divergence-free vector fields for bounded convex domains (see [24, 6]). Denote by D(Ω) the space of smooth vector fields compactly supported on Ω. Let V = {u ∈ D(Ω) : ∇ · u = 0}

Let Lp (Ω) with 1 < p < ∞ denote the standard Lebesgue space and set Lpσ (Ω) = {u ∈ Lp (Ω), ∇ · u = 0 in Ω, u · n = 0 on ∂Ω} .

Alternatively, Lpσ (Ω) is the closure of V in Lp (Ω). Let PΩ : L2 (Ω) → L2τ (Ω) be the orthogonal projector. According to Theorem I.1.4 of [55], for any Lipschitz domain Ω, PΩ is well defined. The following lemma [24] states that PΩ can be extended to general Lp (Ω). Lemma 1.3. Let Ω ⊂ R2 be a bounded convex domain. Let 1 < p < ∞. Then

(1) PΩ can be extended to a bounded linear operator from Lp (Ω) to Lpσ (Ω) and kPΩ ukLp ≤ C(Ω, p) kukLp .

In addition, (1.10)

PΩ u = u for any u ∈ Lpσ (Ω).

(2) For any u ∈ Lp (Ω), there exists a unique π ∈ W 1,p (Ω) (up to an additive constant) such that u = PΩ u + ∇π.

Additionally, in the weak sense, ∆π = ∇ · u, which allows us to write π = ∆−1 ∇ · u,

(1.11)

where π ’s different by an additive constant are regarded the same. In addition, k∇πkLp ≤ C(Ω, p) kukLp .

With this lemma at our disposal, we can now state the linear stability and instability results. e denote Theorem 1.4. Assume the spatial domain Ω ⊂ R2 is as specified in Theorem 1.2. Let (e u, Pe, θ) e 0 ∈ H 2 (Ω), the perturbation as stated in Theorem 1.2 with the initial perturbation (e u0 , θe0 ) satisfying u e 0 = 0 and the compatibility conditions. Let θe0 ∈ H 1 (Ω), ∇ · u PΩ :

L2 (Ω) → L2σ (Ω)

denote the projection operator as in Lemma 1.3. (1) (1.6) is equivalent to the following system without the pressure term,

(1.12)

 e 2 − ∇∆−1 ∇ · (θe e 2 ),  e + PΩ (e ∂t u u · ∇e u) = ν ∆e u + θe     ∂ θe + u  e · ∇θe + ve T 0 (y) = 0,   t e = 0, ∇·u     e  (e u, θ)(x, 0) = (e u0 , θe0 )(x),    u e · n|∂Ω = 0, ω e |∂Ω = 0,

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION

7

with the corresponding linearized system of equations given by

(1.13)

 ∂t U = ν ∆U + Θ e2 − ∇∆−1 ∇ · (Θ e2 ),       ∂ Θ + V T 0 (y) = 0,   t ∇ · U = 0,     (U, Θ)(x, 0) = (U0 (x), Θ0 (x)),     U · n|∂Ω = 0, ω|∂Ω = 0,

where we have used capital letters to distinguish the solutions of the linearized system from those of the full system, and U = (U, V ) with V being the vertical component and ω = ∇ × U. (2) Let T 0 (y) = α > 0. Assume the initial data (U0 , Θ0 ) satisfying U0 ∈ H 2 , ∇ · U0 = 0 and Θ0 ∈ H 1 . When Ω is a rectangle, we also make the natural assumption that Θ0 is zero on the top and the bottom. When Ω is a general domain, we impose that Θ is zero on the boundary ∂Ω. Then (1.13) has a unique global solution (U, Θ) that is also stable in the sense, as t → ∞, • kU(t)k,

(1.14)

k∇U(t)k,

• kΘ(t)k → C0 , • k∇Θ(t)k → C1 ,

k∆U(t)k → 0, p 0 < C0 < αkU0 k2 + kΘ0 k2 , p 0 < C1 < αk∇U0 k2 + k∇Θ0 k2 .

(3) If T 0 (y) = α < 0, then any spatially periodic solution (U, Θ) is unstable.

1.5. Idea of Proof. The core of the proof of Theorem 1.1 is the global-in-time a priori bound for (u, θ) in (1.7). We use L2 -based energy method. The main difficulties are due to the lack of dissipation in the temperature equation and the weak regularity setup. We begin with the Lyapunov functional (free energy formulation) associated with (1.4): (1.15)


d ¯ 2 + 2ανkωk2 = 0, αkuk2 + kθ − (αy + θ)k dt

where α > 0 and θ¯ are arbitrary constants. It follows from (1.15) that the kinetic energy is uniformly bounded with respect to time, and kωk2 is uniformly integrable with respect to time. Since the spatial domain is bounded, the uniform integrability of kωk2 , together with the Poincar´e inequality, implies that the kinetic energy is also uniformly integrable with respect to time. The global bound for k∇uk and the time integrability of kukH 2 are obtained via the vorticity equation. To further the estimates, we seek the global bounds for k∇ωk and k∇θk simultaneously. Due to the lack of dissipation in the temperature equation, we need to control k∇ukL∞ . This is achieved via a logarithmic Sobolev inequality for domains that may not necessarily have smooth boundaries (see Lemma 2.2). In addition, due to the regularity limitation, we also invoke an improved Gronwall inequality (see Lemma 2.3). We make no effort in this proof to optimize the bounds in terms of time. e in Theorem 1.2, we start with a lemma To prove the large-time asymptotics for the perturbation (e u, θ) stating that any nonnegative, integrable and uniformly continuous function f = f (t) on (0, ∞) must decay to zero as t → ∞ (see Lemma 3.1). When the vertically stratified temperature T (y) satisfies 2 e T 0 (y) = α > 0, we first show αke u(t)k2 + kθ(t)k decreases in time t, and k∇e u(t)k2 is integrable on (0, ∞) 2 and ke u(t)k is uniformly continuous in t ∈ [0, ∞). Poincar´e’s inequality implies that ke u(t)k2 ∈ L1 (0, ∞). Lemma 3.1 then leads to ke u(t)k2 → 0 as t → ∞. To show k∇e u(t)k → 0, we resort to the vorticity ω and 2 show kω(t)k is uniformly continuous in t ∈ (0, ∞). An application of Lemma 3.1 then yields the desired e (t)k → 0 as t → ∞ is obtained by applying Lemma 3.1 to the property for k∇e u(t)k . The proof of k∂t u R − 0t kω(τ )k2 dτ 2 e 2 kH −1 → 0 and the uniform e (t)k . The large-time behavior k∇Pe − θe quantity A(t) ≡ e k∂t u e. boundedness of ke u(t)kH 2 are consequences of the decay properties of the velocity u The proof of Theorem 1.4 starts with a derivation of an equivalent system of equations for the perturbation near the hydrostatic equilibrium and the corresponding linearized system. To derive the equivalent

8

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

systems, we need to eliminate the pressure term. Since the domain Ω may not be smooth, we make use of a recent work on the existence of the projection operator PΩ (see [24]) and show that, for the stressfree boundary conditions, PΩ ∆ = ∆, where ∆ denotes the Laplacian operator. With the suitable form of the linearized system at our disposal, the linear stability when T 0 (y) = α > 0 is obtained by careful energy estimates and the application of Lemma 3.1. The linear instability of periodic perturbations when T 0 (y) = α < 0 follows from the fact that one of the eigenvalues of the matrix (in Fourier space) associated with the linear system is always positive. The rest of the paper is organized as follows. Section 2 contains the proof of Theorem 1.1. Section 3 proves Theorem 1.2 while Section 4 proves Theorem 1.4. Section 5 is devoted to the numerical studies of (1.4), in which we try to understand some of the unresolved issues concerning the long-time behavior of the problem. The paper ends with concluding remarks. 2. P ROOF OF T HEOREM 1.1 This section is devoted to the proof of Theorem 1.1. We first state the Calder´on-Zygmund inequality which is frequently used in the subsequent sections. As aforementioned, we assume that Ω is a bounded and simply connected planar domains with ∂Ω belonging to C 2,γ (γ ∈ (0, 1)) except at a finite number of points, where ∂Ω is a corner of angle in (0, π2 ]. Planar polygonal domains such as rectangles in R2 are special examples of such domains. As shown in [6, 25, 39], the Calder´on-Zygmund inequality holds for such domains. Lemma 2.1. Assume that Ω is a bounded and simply connected planar domains with ∂Ω belonging to C 2,γ (γ ∈ (0, 1)) except at a finite number of points, where ∂Ω is a corner of angle in (0, π2 ]. Let v ∈ W s,p (Ω) be a vector-valued function satisfying ∇ · v = 0 and v · n|∂Ω = 0 where n is the unit outward normal to ∂Ω. Then there exists a constant c = c(s, p, Ω), such that kvkW s,p ≤ ck∇ × vkW s−1,p

for any s ≥ 1 and p ∈ (1, ∞). The proof of Theorem 1.1 also involves a logarithmic Sobolev inequality for domains that may not necessarily have smooth boundaries. The following lemma provides a logarithmic Sobolev interpolation inequality for a bounded convex domain with a finite number of corners. Planar polygonal domains such as rectangles in R2 are special examples of such domains. Lemma 2.2. Let Ω ⊂ R2 be a Lipschitz domain verifying that there exist constants A > 1, δ0 ∈ (0, diam(Ω)) such that for any r ∈ (0, δ0 ) and x ∈ Ω, A|Ω(x, r)| ≥ |B(x, r)| ≥ |Ω(x, r)|,

where Ω(x, r) = Ω ∩ B(x, r) with B(x, r) being the ball with center x and radius r. Then there exists a constant C depending on δ0 such that, for any smooth function f and 0 < α, β ≤ 1,   kf kLp (Ω) + 1 logα (e + kf kC 0,β (Ω) ). kf kL∞ (Ω) ≤ C β −α sup pα p≥2

A special consequence is the following logarithmic Sobolev interpolation inequality. Assume Ω is a bounded and simply connected domain of R2 , with ∂Ω belonging to C 2,γ (γ ∈ (0, 1)) except at a finite number of points, where ∂Ω is a corner of angle between 0 and π2 . Assume u is smooth vector. Then, for some constant C , k∇ukL∞ (Ω) ≤ C (1 + k∇ukH 1 (Ω) ) log(1 + k∆∇ukL2 (Ω) ).

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION

9

The first part of Lemma 2.2 can be found in [26]. The second part is a special consequence of the first part. In fact, if we take α = 1 and any β < 1, we obtain the second part by recalling the Sobolev inequalities, kf kLp (Ω) ≤ C p k∇f kL2 (Ω) ,

kf kC 0,β (Ω) ≤ C kf kH 2 (Ω) ,

which hold for Lipschitz domains. We also need an improved Gronwall type inequality (see, e.g., [44]).

Lemma 2.3. Assume that Y, Z, A and B are non-negative functions satisfying d Y (t) + Z(t) ≤ A(t) Y (t) + B(t) Y (t) ln(1 + Z(t)), dt Let T > 0. Assume A ∈ L1 (0, T ) and B ∈ L2 (0, T ). Then, for any t ∈ [0, T ],

(2.1)

Rt

Y (t) ≤ (1 + Y (0))e 0

(2.2)

and (2.3)

Z

t

0

Z(τ ) dτ ≤ Y (t)

Z

B(τ ) dτ

e

Rt 0

e

Rt s B(τ ) dτ (A(s)+B 2 (s)) ds

t

A(τ ) dτ + Y 2 (t)

0

Z

t

0

B 2 (τ ) dτ < ∞.

For the convenience of the readers, we provide a proof of this lemma. Proof of Lemma 2.3. Setting Y1 (t) = ln(1 + Y (t)),

Z1 (t) = Z(t)/(1 + Y (t)),

we have d Y1 (t) + Z1 (t) ≤A(t) + B(t) ln(1 + Z(t)) dt ≤A(t) + B(t) ln(1 + (1 + Y (t)) Z1 (t)) ≤A(t) + B(t) ln(1 + Y (t))(1 + Z1 (t))

≤A(t) + B(t) Y1 (t) + B(t) ln(1 + Z1 (t)).

Using the simple fact that, for f ≥ 0,

1

(2.4)

we obtain

ln(1 + f (t)) ≤ f 2 (t), d 1 Y1 (t) + Z1 (t) ≤ A(t) + B(t) Y1 (t) + B 2 (t) + Z1 (t). dt 4

Gronwall’s inequality then implies

Y1 (t) ≤ Y1 (0) e

Rt 0

B(τ ) dτ

+

Z

0

t R t

e

s

B(τ ) dτ

(A(s) + B 2 (s)) ds,

which yields (2.2). In addition, (2.2) allows us to obtain (2.3) by using the inequality (2.4) in (2.1) and integrating in time. This completes the proof of Lemma 2.3.  We are now ready for the proof of Theorem 1.1. Proof of Theorem 1.1. The key component of the proof is the global-in-time a priori bound for (u, θ) in the Sobolev space u ∈ L∞ ((0, t); H 2 (Ω)) ∩ L2 ((0, t); H 3 (Ω)),

θ ∈ L∞ ((0, t); H 1 (Ω))

for any t > 0. Once the global bound is established, the existence of solutions then follows from some of the standard procedures such as various approximation schemes. The uniqueness of the solutions is a simple consequence of the fact Z

0

t

k∇u(·, τ )kL∞ (Ω) dτ < ∞

10

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

due to u ∈ L2 ((0, t); H 3 (Ω)). This remark explains why it suffices to prove the global a priori bounds. In general the bounds depend on the time t and we shall not make efforts to obtain the optimal bounds. First of all, since θ is purely transported by a divergence-free vector field, it holds that (2.5)

kθ(t)kLp = kθ0 kLp , ∀ t > 0, ∀ p ∈ [2, ∞].

For any fixed constants a > 0 and b, c ∈ R, we write θ˜ = θ − ay − b and P˜ = P − a2 y 2 − by − c to convert (1.4) to  ˜ 2, ∂t u + u · ∇u + ∇P˜ = ν∆u + θe       ∂t θ˜ + u · ∇θ˜ = −au · e2 ,  ∇ · u = 0;    ˜  (u, θ)(x, 0) = (u0 , θ0 − ay − b)(x),     u · n|∂Ω = 0, ω|∂Ω = 0,

(2.6)

Taking the L2 -inner products of the first equation of (2.6) with au and the second equation with θ˜, and using the boundary conditions, we have  1 d  ˜ 2 + aνkωk2 = 0, akuk2 + kθk 2 dt

(2.7)

which implies, in terms of the original functions, that


1 d akuk2 + kθ − ay − bk2 + aνkωk2 = 0. 2 dt

(2.8)

This provides a Lyapunov functional associated with (1.4). Moreover, since d kθk2 = 0 dt

the Lyapunov functional is equivalent to d dt



d dt

and

1 kuk2 − 2

Z

Z

θ dx = 0,

Ω



yθ dx + νkωk2 = 0,

Ω

where the quantity inside the parenthesis is the total mechanical energy associated with (1.4). Our numerical simulations, see Section 5, illustrate that the total mechanical energy decreases as time evolves. Upon integrating (2.8) with respect to time, we obtain (2.9)

aku(t)k2 + kθ(t) − ay − bk2 + 2aν

which implies (2.10)

ku(t)k2 + ν

Z

Z

0

t

kω(τ )k2 dτ = aku0 k2 + kθ0 − ay − bk2 , ∀ t ≥ 0,

t

0

kω(τ )k2 dτ + ν

Z

0

t

k∇u(τ )k2 dτ ≤ c1 , ∀ t ≥ 0.

To prove the global H 1 -bound, we resort to the vorticity equation (2.11)

∂t ω + u · ∇ω = ν∆ω + ∂x θ.

Multiplying (2.11) by ω , integrating over Ω and applying the boundary conditions u · n|∂Ω = 0 and ω|∂Ω = 0, we have (2.12)

Inserting the inequality

1 d kωk2 + νk∇ωk2 = − 2 dt

Z

θ ∂x ω dx.

Ω

Z θ ∂x ω dx ≤ ν k∇ωk2 + c2 kθk2 2 Ω

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 11

in (2.12) and integrating in time yields kω(t)k2 + ν

(2.13)

Z

0

t

k∇ω(τ )k2 dτ ≤ 2c2 kθ0 k2 t.

According to Lemma 2.1, (2.10) and (2.13) especially imply Z

(2.14)

0

t

ku(τ )kH 2 dτ ≤ c3 t.

Next we prove a global bound for kukH 2 and k∇θk. This is achieved by combining the energy inequalities for k∇ωk and k∇θk and invoking Lemma 2.2 and Lemma 2.3. Multiplying the vorticity equation (2.11) by ∆ω , integrating over Ω and applying the boundary conditions u · n|∂Ω = 0 and ω|∂Ω = 0, we have Z Z 1 d k∇ωk2 + νk∆ωk2 = − 2 dt

By the Ladyzhenskaya inequality,

Ω

∂x θ ∆ω dx −

Ω

(∇ω)T ∇u ∇ω dx.

Z (∇ω)T ∇u ∇ω dx ≤k∇uk k∇ωk2 4 L Ω

≤c4 k∇uk k∇ωk k∆ωk ν ≤ k∆ωk2 + c5 k∇uk2 k∇ωk2 . 4

Therefore, (2.15)

d k∇ωk2 + νk∆ωk2 ≤ c6 k∇θk2 + c7 k∇uk2 k∇ωk2 . dt

It follows from the temperature equation that d k∇θk2 ≤ 2k∇ukL∞ k∇θk2 . dt

Invoking the logarithmic Sobolev inequality in Lemma 2.2, we have d k∇θk2 ≤ c8 (1 + kωkH 1 log(1 + k∆ωk)) k∇θk2 . dt Combining (2.15) and (2.16) and setting Y (t) = k∇ωk2 + k∇θk2 and Z(t) = k∆ωk2 , we obtain (2.16)

d Y (t) + ν Z(t) ≤ c9 (1 + k∇uk2 ) Y (t) + c10 kukH 2 Y (t) log(1 + Z(t)). dt Applying Lemma 2.3 to (2.17) and recalling the bound for kukH 2 in (2.14), we obtain the desired global bound for k∇ωk and k∇θk, Z t c13 t k∇ω(t)k, k∇θ(t)k, k∆ω(t)k2 dτ ≤ c11 ec12 e .

(2.17)

0

By Lemma 2.1 and the Sobolev embedding, we have, for 2 < p < ∞, (2.18)

Z

0

t

ku(τ )k2W 2,p dτ ≤ c14 ec12 e

c13 t

.

We have made no effort here to optimize the time growth rate. We derive further energy estimates for the solution when (u0 , θ0 ) ∈ H 3 , and complete the proof of Theorem 1.1. We start with the Lp -estimate of ∇θ. By taking the gradient ∇ of the temperature equation, then testing the result with p|∇θ|p−2 ∇θ, we have

d k∇θkpLp ≤ pk∇ukL∞ k∇θkpLp , 2 < p < ∞, dt which, together with the Sobolev embedding, L∞ ,→ W 1,p (2 < p < ∞), implies (2.19)

d k∇θkLp ≤ c15 kukW 2,p k∇θkLp , 2 < p < ∞. dt

12

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

We remark that no boundary condition for θ is necessary in the derivation of (2.19). In a similar fashion, we can show that (2.20)

d kD2 θkLp ≤ c16 k∇θkW 1,p kukW 2,p + c17 kukW 2,p kD2 θkLp , 2 < p < ∞. dt

Combining (2.19) and (2.20), we obtain d k∇θkW 1,p ≤ c18 kukW 2,p k∇θkW 1,p , 2 < p < ∞. dt

Gronwall’s inequality and (2.18) then imply (2.21)

k∇θkW 1,p ≤ ec19 e

c12 ec13 t

k∇θ0 kW 1,p , 2 < p < ∞.

Moreover, direct calculations show that kθt kLp ≤ kukL∞ k∇θkLp

≤ c20 kukH 2 k∇θkLp ≤ c21 ec12 e

c13 t

k∇θkLp ,

k∇θt kLp ≤ c22 k∇ukH 1 k∇θkL2p + c23 kukH 2 kD2 θkLp
c13 t ≤ c24 ec12 e k∇θkL2p + kD2 θkLp , 2 < p < ∞,

which, together with (2.21), shows that (2.22)

kθt kW 1,p ≤ c25 ec19 e

c12 ec13 t

, 2 < p < ∞.

In addition, by taking the temporal derivative to the equation of vorticity, then calculating the L2 inner product of the resulting equation with ∂t ω , we deduce 1 d k∂t ωk2 + νk∇∂t ωk2 = − 2 dt (2.23)

Z

Ω

(∂t u · ∇ω)∂t ω dx +

Z

(∂t ∂x θ)∂t ω dx

Ω

1 1 ≤ k∂t ukL4 k∇ωkL4 k∂t ωk + k∂t ωk2 + k∂t ∂x θk2 2 2 1 1 2 2 ≤ c26 kukW 2,4 k∂t ωk + k∂t ωk + k∂t ∂x θk2 , 2 2

where we have applied Gagliardo-Nirenberg interpolation inequality and Lemma 2.1. Applying the Gronwall inequality to (2.23) and using (2.18) and (2.22), we find (2.24)

∂t ω ∈ L∞ (0, T ; L2 (Ω)) ∩ L2 (0, T ; H 1 (Ω)), ∀ 0 < T < ∞.

Applying the classic elliptic regularity theory to the equation of vorticity and using (2.21) and (2.24), we can easily show that ω ∈ L∞ (0, T ; H 2 (Ω)) ∩ L2 (0, T ; H 3 (Ω)), ∀ 0 < T < ∞,

which, together with Lemma 2.1, implies (2.25)

u ∈ L∞ (0, T ; H 3 (Ω)) ∩ L2 (0, T ; H 4 (Ω)), ∀ 0 < T < ∞.

Finally, applying (2.25) to the equation of temperature, we can show that θ ∈ L∞ (0, T ; H 3 (Ω)), ∀ 0 < T < ∞.

This completes the proof of Theorem 1.1.



LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 13

3. P ROOF OF T HEOREM 1.2 This section proves Theorem 1.2. The following decay lemma will be used several times. Lemma 3.1. Assume f ∈ L1 (0, ∞) is a nonnegative and uniformly continuous function. Then, as t → ∞.

f (t) → 0

In particular, if f ∈ L1 (0, ∞) is nonnegative and satisfies, for a constant C and any 0 ≤ s < t < ∞, |f (t) − f (s)| ≤ C |t − s|,

then f (t) → 0 as t → ∞. The proof of this lemma is not difficult. Proof. We prove by the definition of limit. For any  > 0, the uniform continuity of f implies that there exists δ > 0 such that, for any 0 ≤ s < t < ∞, |f (t) − f (s)| <

 2

whenever

Since f ∈ L1 (0, ∞), there exists T > 0 such that Z

∞

f (t)dt <

T

Especially, for any k = 0, 1, 2, · · · ,

Z

δ . 2

T +(k+1)δ

f (t)dt <

T +kδ

|t − s| < δ.

δ , 2

which implies that there exists tk ∈ (T + kδ, T + (k + 1)δ) such that

 . 2 Then, for any t > T , there exists k0 satisfying t ∈ (T + k0 δ, T + (k0 + 1)δ] and thus   f (t) ≤ f (tk0 ) + |f (t) − f (tk0 )| ≤ + = . 2 2 f (tk ) <



This concludes the proof of Lemma 3.1.

Although primarily we focus on the case when Ω is a rectangle, our results hold for more general domains as specified in Theorem 1.1. When Ω is a rectangle, the curvature of the boundary κ is zero and no boundary term gets into play. For a more general domain, a boundary term related to κ plays a role. We will use n to denote the unit outward normal vector to ∂Ω. If s is arc length and we parameterize the dn dn curve by arc length, then dn ds · n = 0; i.e., ds is perpendicular to n. The magnitude of ds is called the dn curvature of the curve and is denoted by κ. That is to say, κ = | ds |. The unit vector in the direction of dn ds is defined as the unit tangent vector to the curve and is denoted as τ . The following lemma facilitates the process of energy estimates in the proof of Theorem 1.2. Lemma 3.2. Let Ω be a rectangle or a more general domain as specified in Theorem 1.1. Assume the vector field u and its corresponding ω = ∇ × u satisfy the boundary conditions u · n|∂Ω = 0,

Then the following identity holds

ω|∂Ω = 0.

(n · ∇u) · τ + κ u · τ = 0 on ∂Ω,

where n and τ denote the unit outward normal and tangent vector to ∂Ω, respectively. We provide a proof of Lemma 3.2.

14

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

Proof. Since u · n = 0 on ∂Ω, the directional derivative of u · n along ∂Ω should also be zero, namely d (u · n) = 0 on ∂Ω. dτ

The product rule then yields 

   d d u ·n+u· n = 0 or τ · ∇u · n + u · (τ · ∇n) = 0 dτ dτ

Due to u · n = 0 on ∂Ω, Therefore, due to κ =

dn ds

on

∂Ω.

u = (u · n)n + (u · τ )τ = (u · τ )τ on ∂Ω.     ∂n dx1 ∂n dx2 ∂n ∂n · τ = ∂x + · τ = τ + τ ∂x2 ds ∂x1 1 ∂x2 2 · τ = τ · ∇n · τ , 1 ds

τ · ∇u · n + (τ · ∇n · τ ) (u · τ ) = 0 or τ · ∇u · n + κ u · τ = 0 on

Since ω = 0 on ∂Ω,

∇u − (∇u)T =

That is, ∇u = (∇u)T . Therefore,



0 ω

−ω 0



∂Ω.

= 0 on ∂Ω.

(n · ∇u) · τ = τ · (∇u)T · n = τ · ∇u · n = − κ u · τ,

which is the desired equality.



We now prove Theorem 1.2. Proof of Theorem 1.2. Since the global-in-time existence and uniqueness can be similarly obtained as in the proof of Theorem 1.1, we focus on the asymptotic behavior. Assume T 0 (y) = α > 0. Taking the inner e , we have product of the first two equations in (1.6) with (e u, θ)  1 d  e2 αke uk2 + kθk 2 dt

(3.1)

Z e dx − α (e e dx ∆e u ·u u · ∇e u) · u ZΩ Z Ω e ·u e dx − (e e dx. −α ∇Pe · u u · ∇θ)

= να

Z

Ω

Ω

The boundary conditions allow us to eliminate the last three terms, Z

(3.2) (3.3) (3.4)

Z 1 e |e e dx = (e u · ∇e u) · u n·u u|2 dS(x) = 0, 2 ∂Ω Ω Z Z e dx = e Pe dS(x) = 0, ∇Pe · u n·u ∂Ω ZΩ Z 1 e ·u e dx = e θe2 dS(x) = 0. (e u · ∇θ) n·u 2 Ω ∂Ω

For the dissipative term, we first apply the divergence theorem to obtain Z

Ω

e dx = ∆e u ·u

e · n = 0 on ∂Ω, we can write Since u

∂Ω

Therefore, (3.5)

∂Ω

e dS(x) − (n · ∇)e u·u

e = (e u u · n)n + (e u · τ )τ = (e u · τ )τ

By Lemma 3.2, Z

Z

e dS(x) = (n · ∇)e u·u Z

Ω

Z

∂Ω

Z

Ω

|∇e u|2 dx.

on ∂Ω.

(n · ∇)e u · τ (e u · τ ) dS(x) = −

e dx = − ∆e u ·u

Z

∂Ω

κ (e u · τ )2 dS(x) −

Z

Ω

Z

∂Ω

κ (e u · τ )2 dS(x).

|∇e u|2 dx.

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 15

In the special case when Ω is a rectangle, the boundary term vanishes and Z

Ω

e dx = − ∆e u ·u

Inserting (3.2), (3.3), (3.4) and (3.5) in (3.1) yields (3.6)

 1 d  e 2 +να αke uk2 + kθk 2 dt

Z

Z

Ω

|∇e u|2 dx.

2

Ω

|∇e u| dx + ν α

Integrating in time yields, for any 0 ≤ s < t < ∞,

(3.7)

2 e αke u(t)k2 + kθ(t)k +να

2 e ≤ αke u(s)k + kθ(s)k . 2

Z

t

s

k∇e u(ρ)k2 dρ + ν α

Z

∂Ω

Z

s

t

κ (e u · τ )2 dS(x) = 0.

√ u · τ )(ρ)k2L2 (∂Ω) dρ k κ (e

2 e Especially, (3.7) implies, Y (t) ≡ αke u(t)k2 + kθ(t)k is a decreasing function of t ∈ [0, ∞) and

0 < Y (t) ≤ Y0 ≡ αke u0 k2 + kθe0 k2 .

As a consequence, there exists c0 > 0 such that, as t → ∞, (3.8)

Next we show that, as t → ∞,

2 e Y (t) ≡ αke u(t)k2 + kθ(t)k → c20 < Y02 .

ke u(t)k → 0

and

k∇e u(t)k → 0.

We recall the Poincar´e inequality, for any u ∈ H (Ω) with ∇ · u = 0 and u · n = 0 on ∂Ω, 1

kukL2 (Ω) ≤ C(Ω) k∇ukL2 (Ω) .

Invoking the Poincar´e inequality in (3.7) yields, for C > 0, αke u(t)k2 + C

Z

∞

0

which in particular implies

ke u(τ )k2 dτ ≤ αke u0 k2 + kθe0 k2 ,

ke u(t)k2 ∈ L1 (0, ∞).

(3.9)

We also need another ingredient, the uniform continuity of ke u(t)k2 as a function of t ∈ [0, ∞). Taking e in (1.6) with u e , we obtain the inner product of the equation of u 1 d ke uk2 + ν α k∇e uk2 ≤ 2 dt

Therefore, for any 0 ≤ s ≤ t < ∞ and C = (3.10)

Y0 √ , α

Z

ev dx ≤ ke e ≤ Y √0 . θe uk kθk 2 α

ke u(t)k2L2 − ke u(s)k2L2 ≤ C |t − s|.

(3.9) and (3.10), together with Lemma 3.1, lead to the desired decay (3.11)

In addition, (3.8) and (3.11) yield We now prove (3.12)

ke u(t)k → 0

as t → ∞.

e kθ(t)k → c0 . k∇e u(t)k → 0

as t → ∞.

According to (3.7), k∇e u(t)k2 ∈ L1 (0, ∞), which, in particular, implies

(3.13)

ke ω (t)k2 ∈ L1 (0, ∞).

Taking the inner product of the vorticity equation

e e · ∇e ∂t ω e+u ω = ν∆e ω + ∂1 θ.

16

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

with ω and invoking the boundary conditions yield 1 d ke ω k2 + νk∇e ω k2 = − 2 dt

or

Z

ν 1 e 2 ω k2 + kθ0 k ∂1 ω θe dx ≤ k∇e 2 2ν

d 1 ke ω k2 + νk∇e ω k2 ≤ kθe0 k2 . dt ν Integrating in time yields, for any 0 ≤ s < t < ∞, 1 ke (3.14) ω (t)k2 − ke ω (s)k2 ≤ |t − s|kθe0 k2 . ν

(3.13) and (3.14), and Lemma 3.1 allow us to conclude that (3.15)

as t → ∞.

ke ω (t)k → 0

By Lemma 2.1, k∇e u(t)k ≤ C ke ω (t)k and (3.15) implies (3.12). Next we show that, as t → ∞, (3.16)

e (t)k → 0. k∂t u

e in (1.6) yields Applying ∂t to the equation of u

e 2. e+u e · ∇∂t u e+u e t · ∇e e + ∂t θe ∂tt u u + ∇∂t Pe = ν∆∂t u

(3.17)

e and invoking the equation of θe in (1.6), we have Taking the L2 -inner product of (3.17) with ∂t u

(3.18)

1 d e (t)k2 + νk∂t ω k∂t u e k2 = 2 dt

Z

e (e2 · ∂t u e ) dx − (∂t θ)

Z

e · ∇e e ) dx (∂t u u) · (∂t u Z eu · ∇(e2 · ∂t u e ) dx − a e ) dx = θe ve(e2 · ∂t u Ω Ω Z e · ∇e e ) dx. − (∂t u u) · (∂t u ZΩ

Ω

Ω

By the H¨older inequality, Poincar´e inequality and Sobolev inequality,

Z ν eu · ∇(e2 · ∂t u e L∞ ke θe ≤ kθk e e k ≤ k∂t ω e k2 + Cke uk2 , ) dx uk k∇∂t u 8 Ω Z ν ≤ ke ve(e2 · ∂t u e k ≤ C ke e k ≤ k∂t ω e v k k∂t u v k k∇∂t u ) dx e k2 + Cke v k2 8 Ω

and

Z (∂t u e · ∇e e ) dx ≤ k∇e e k2L4 ≤ C kωk k∂t u e kk∇∂t u ek u) · (∂t u uk k∂t u Ω ν e k2 . ≤ k∂t ω e k2 + Ckωk2 k∂t u 4

Inserting the estimates above in (3.18) yields (3.19)

d e (t)k2 + νk∂t ω e k2 + C ke k∂t u e k2 ≤ C kωk2 k∂t u uk2 . dt

Gronwall’s inequality, together with (3.9) and (3.13), implies

Poincar´e’s inequality then implies (3.20)

Z

0

Writing (3.19) as

e (t)k2 + ν k∂t u

∞

Z

t

0

e (t)k2 dt ≤ C k∂t u

k∂τ ω e (τ )k2 dτ ≤ C. Z

0

∞

k∂t ω e (t)k2 dt ≤ C.

 Rt 2 d  −C R t kωk2 dτ 0 e k2 ≤ C e−C 0 kωk dτ ke e k∂t u uk2 , dt

integrating in time and using the fact the right-hand side is bounded uniformly in time, we obtain (3.21)

|A(t) − A(s)| ≤ C |t − s|,

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 17

where A is given by A(t) = e−C

In addition, (3.20) implies

Rt 0

kω(τ )k2 dτ

e (t)k2 . k∂t u

A(t) ∈ L1 (0, ∞).

Lemma 3.1 then asserts A(t) → 0 as t → ∞. Therefore, (3.22)

e (t)k2 = eC k∂t u

Rt 0

kω(τ )k2 dτ

A(t) ≤ C A(t) → 0.

e in (1.6) that It then follows from the equation of u

e 2 kH −1 ≤k∂t u e k + ke e k + k∇e u⊗u uk k∇Pe − θe

(3.11), (3.12) and (3.22) then lead to

e k + C ke ≤k∂t u uk k∇e uk + k∇e uk.

e 2 kH −1 → 0 k∇Pe − θe

as t → ∞.

e in (1.6) can be written in the component form as Finally we prove (1.9). We note that the equation for u e · ∇e ν∂y ω e + ∂x Pe = −∂t u e−u u,

e e · ∇e −ν∂x ω e + ∂y Pe = −∂t ve − u v + θ.

Due to the boundary condition ω e |∂Ω = 0, it holds that

kν∂y ω e + ∂x Pek2 + k − ν∂x ω e + ∂y Pek2 = ν 2 k∇e ω k2 + k∇Pek2 + 2ν

= ν 2 k∇e ω k2 + k∇Pek2 + 2ν = ν 2 k∇e ω k2 + k∇Pek2 ,

Z

ZΩ

Ω

(∂y ω e ∂x Pe − ∂x ω e ∂y Pe)dx ∇ · (−e ω ∂y Pe, ω e ∂x Pe)dx

which, together with the preceding equations, implies (3.23)

e 2. e k2 + ke k∇e ω k2 . k∂t u u · ∇e uk2 + kθk

According to the Ladyzhenskaya inequality, Lemma 2.1 and Poincar´e inequality, we have ke u · ∇e uk2 . ke uk2L4 k∇e uk2L4

. ke ukk∇e uk2 k∇e ωk .

1 k∇e ω k2 + ke uk2 k∇e uk4 , 2

which, together with (3.23), implies (3.24)

e 2. e k2 + ke k∇e ω k2 . k∂t u uk2 k∇e uk4 + kθk

e 2 ≤ 2(kθk2 + kαy + θk ¯ 2 ), the uniform boundedness of k∇ω(t)k then follows from (3.24), (3.11), Since kθk (3.12), (3.22) and (2.5). As a consequence of Lemma 2.1 and the Poincar´e inequality, we conclude that ku(t)kH 2 is uniformly bounded in time. This completes the proof of Theorem 1.2. 

4. P ROOF OF T HEOREM 1.4 This section proves Theorem 1.4. For the sake of clarity, it is divided into three subsections. The first subsection derives the system of equations (1.13) while the second proves the linear stability when T 0 (y) = α > 0 and the third proves the linear instability of periodic perturbations when T 0 (y) = α < 0.

18

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

e ∈ H 2 (Ω) and θe ∈ H 1 (Ω), we present a 4.1. Derivation of (1.13). Although our functional setting is u derivation that is valid for more general Sobolev spaces e ∈ W 2,p (Ω), u

θe ∈ W 1,p (Ω),

Let PΩ denote the projection operator from L (Ω) to Applying PΩ to the velocity equation in (1.6) yields p

Lpσ (Ω)

1 < p < ∞.

with 1 < p < ∞ as stated in Lemma 1.3.

e 2 ). e + PΩ (e ∂t u u · ∇e u) = ν PΩ ∆e u + PΩ (θe

(4.1)

We first show that (4.2)

PΩ ∆e u = ∆e u.

In general the operators PΩ and ∆ do not commutate and (4.2) may not be true, but the stress-free e ∈ W 2,p (Ω) with ∇ · u e = 0, we have boundary conditions in (1.6) allow us to prove (4.2). For u ∆e u ∈ Lp (Ω)

(4.3)

In addition, we show that (4.4)

and ∇ · ∆e u = 0 in Ω.

∆e u · n = 0 on ∂Ω.

In fact, ω e = 0 on ∂Ω implies the directional derivative of ω e along the tangential direction τ on ∂Ω is zero, namely de ω = 0 on ∂Ω, dτ where τ denotes the unit tangent vector to ∂Ω. That is,

(4.5)

τ · ∇e ω = 0 or τ1 ∂1 ω e + τ2 ∂2 ω e = 0 on ∂Ω.

Since τ · n = 0, we have τ = n⊥ or (τ1 , τ2 ) = (−n2 , n1 ). Thus, (4.5) becomes Invoking ω e = ∂1 ve − ∂2 u e yields

−n2 ∂1 ω e + n1 ∂2 ω e = 0 on ∂Ω.

−n2 ∂1 (∂1 ve − ∂2 u e) + n1 ∂2 (∂1 ve − ∂2 u e) = 0

or

which is exactly (4.4). Now (4.3) and (4.4) imply that

− n2 ∆e v − n1 ∆e u = 0 on ∂Ω,

∆e u ∈ Lpσ (Ω).

(1.10) in Lemma 1.3 allows us to conclude (4.2),

PΩ ∆e u = ∆e u. e 2 ). According to Lemma 1.3, we can write, for some Next we derive an explicit formula for PΩ (θe 1,p h∈W , e 2 = PΩ (θe e 2 ) + ∇h, θe

where the operator ∆−1 is defined in (1.11). Therefore, (4.6)

e 2 ). h = ∆−1 ∇ · (θe

e 2 ) = θe e 2 − ∇∆−1 ∇ · (θe e 2 ). PΩ (θe

Inserting (4.2) and (4.6) in (4.1) yields

e 2 − ∇∆−1 ∇ · (θe e 2 ). e + PΩ (e ∂t u u · ∇e u) = ν ∆e u + θe

We thus have obtained the following equivalent form of (1.6) (4.7)

 e 2 − ∇∆−1 ∇ · (θe e 2 ),  e + PΩ (e ∂u u · ∇e u) = ν ∆e u + θe   t e · ∇θe + ve T 0 (y) = 0, ∂t θe + u    e = 0. ∇·u

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 19

Ignoring the nonlinear terms in (4.7), we obtain the associated linearized system  −1  ∂t U = ν ∆U + Θ e2 − ∇∆ ∇ · (Θ e2 ),     ∂ Θ + V T 0 (y) = 0,   t ∇ · U = 0,     (U, Θ)(x, 0) = (U0 (x), Θ0 (x)),     U · n|∂Ω = 0, ω|∂Ω = 0,

(4.8)

where we have used capital letters to distinguish the solutions of the linearized system from those of the full system (4.7), and U = (U, V ) with V being the vertical component and ω = ∇ × U. 4.2. Linear stability for T 0 (y) = α > 0. This subsection proves the second part of Theorem 1.4. Assume T 0 (y) = α > 0 and the conditions in the second part of Theorem 1.4. The existence and uniqueness of the corresponding solution can be proven similarly as in the proof of Theorem 1.1. It suffices to prove (1.14). As in the proof of (3.6) in Theorem 1.2, we have
√ 1 d αkU(t)k2 + kΘ(t)k2 + ναk∇U(t)k2 + ναk κ (U · τ )(t)k2L2 (∂Ω) = 0, 2 dt where κ denotes the curvature of ∂Ω and κ = 0 in the case of a rectangular domain, and we have invoked

the fact that

Z

Ω

U · ∇∆

−1

∇ · (Θe2 ) dx =

Integrating in time yields 2

2

αkU(t)k + kΘ(t)k + να ≤ αkU0 k2 + kΘ0 k2 ,

(4.9)

Z

Z

∂Ω

U · n ∆−1 ∇ · (Θe2 ) dS(x) = 0.

t 2

0

k∇U(s)k ds + να

which yields, as in the proof of Theorem 1.2, kU(t)k → 0,

In addition, (4.9) also implies

kΘ(t)k → C0 <

k∇U(t)k2 ∈ L1 (0, ∞) and

(4.10)

It remains to prove the limits

k∇U(t)k → 0,

Z

k∇Θ(t)k → C1 ,

Ω

∆U · ∂t U dx = ν k∆Uk2 +

Z

Ω

0

t

√ k κ (U · τ )(s)k2L2 (∂Ω) ds

p αkU0 k2 + kΘ0 k2 .

√ k κ (U · τ )(t)k2L2 (∂Ω) ∈ L1 (0, ∞).

Dotting the equation of U in (1.13) with ∆U yields (4.11)

Z

k∆U(t)k → 0

Θ∆V dx −

Z

as t → ∞.

∆U · ∇∆−1 ∇ · (Θ e2 )) dx.

Ω

We integrate by parts in the first term. To be more precise, we write the term in the form of components and apply the divergence theorem, Z

Ω

∆U · ∂t U dx

=

Z

ZΩ

∂k ∂k U · ∂t U dx

Ω

∂k (∂k U · ∂t U) −

Z

∂k U · ∂t ∂k U dx Z Z 1 d = nk ∂k U · ∂t U dS(x) − |∇U|2 dx, 2 dt Ω ∂Ω =

Ω

where n = (n1 , n2 ) denotes the unit outward normal to ∂Ω. To further simplify it, we write, due to U · n = 0 on ∂Ω, U = (U · τ )τ + (U · n)n = (U · τ )τ on ∂Ω

20

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

and invoke Lemma 3.2 to obtain Z

Ω

∆U · ∂t U dx

(4.12)

Z

Z 1 d |∇U|2 dx 2 dt Ω ∂Ω Z Z 1 d |∇U|2 dx = − κ(U · τ ) ∂t (U · τ ) dS(x) − 2 dt Ω ∂Ω Z  Z 1 d = − κ (U · τ )2 dS(x) + |∇U|2 dx . 2 dt ∂Ω Ω

=

(n · ∇)U · τ ∂t (U · τ ) dS(x) −

The last term in (4.11) is actually zero. In fact, by ∇ · U = 0 and the divergence theorem Z

Ω

∆U · ∇∆−1 ∇ · (Θ e2 )) dx

Z

=

Ω

∆Uk ∂k ∆−1 ∇ · (Θ e2 )) dx

Ω

∂k (∆Uk ∆−1 ∇ · (Θ e2 )) dx

Z

=

Z

=

∂Ω

∆U · n ∆−1 ∇ · (Θ e2 ) dx.

Going through a similar process as in the derivation of (4.4), we have ∆U · n = 0 on ∂Ω

and thus (4.13)

Z

∆U · ∇∆−1 ∇ · (Θ e2 )) dx = 0.

Ω

By the divergence theorem, Z

Θ∆V dx

=

Ω

=

Z

ZΩ

∇ · (Θ∇V ) dx −

∂Ω

Z

∇Θ · ∇V dx Z Θ n · ∇V dS(x) − ∇Θ · ∇V dx. Ω

Ω

In the case when Ω is a rectangle, n · ∇V = n1 ∂1 V + n2 ∂2 V = 0 on the two sides of Ω, and on the top and the bottom of Ω, Z t

Θ(x, t) = Θ0 (x) +

V (x, τ ) dτ = 0

0

due to V (x, τ ) = 0 and Θ0 (x) = 0 on the top and the bottom. When Ω is a general domain, we need the condition that Θ(x, t) = 0 for x ∈ ∂Ω. Therefore, Z

(4.14)

Ω

Θ∆V dx = −

Z

Ω

∇Θ · ∇V dx.

Inserting (4.12), (4.13) and (4.14) in (4.11), we obtain (4.15)

1 d 2 dt

Z

∂Ω

κ (U · τ )2 dS(x) +

Z

Ω

 Z |∇U|2 dx + νk∆Uk2 = ∇Θ · ∇V dx. Ω

Taking the gradient of the equation of Θ in (1.13) and then dotting the resulting equation with ∇Θ, we obtain 1 d k∇Θk2 + α 2 dt

(4.16)

Z

Ω

∇Θ · ∇V dx = 0.

Multiplying (4.15) by α and then adding to (4.16) lead to (4.17)

 √ 1 d  αk∇Uk2 + αk κ (U · τ )k2L2 (∂Ω) + k∇Θk2L2 (Ω) + νk∆Uk2 = 0. 2 dt

Integrating in time yields

(4.18)

Z t √ αk∇Uk2 + αk κ (U · τ )k2L2 (∂Ω) + k∇Θk2L2 (Ω) + 2ν k∆Uk2 dτ 0 √ = αk∇U0 k2 + αk κ (U0 · τ )kL2 (∂Ω) + k∇Θ0 k2L2 (Ω) ≡ Z02 .

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 21

As a consequence, (4.16) implies √ d k∇Θk2 ≤ α k∇Θkk∇Uk ≤ αZ02 dt

and thus, for any 0 ≤ s < t < ∞,

√ k∇Θ(t)k2 − k∇Θ(s)k2 ≤ αZ02 |t − s|.

In addition, integrating (4.17) in time from s to t yields

√ αk∇U(t)k2 + αk κ (U · τ )(t)k2L2 (∂Ω) + k∇Θ(t)k2L2 (Ω) √ ≤ αk∇U(s)k2 + αk κ (U · τ )(s)k2L2 (∂Ω) + k∇Θ(s)k2L2 (Ω) ,

which especially implies

(4.19)

    √ √ α k∇U(t)k2 + k κ (U · τ )(t)k2L2 (∂Ω) − k∇U(s)k2 + k κ (U · τ )(s)k2L2 (∂Ω) √ ≤ k∇Θ(t)k2 − k∇Θ(s)k2 ≤ αZ02 |t − s|.

(4.10) and (4.19), together with Lemma 3.1, allow us to conclude that and

k∇U(t)k → 0

√ k κ (U · τ )(t)k2L2 (∂Ω) → 0

as t → ∞

and k∇Θ(t)k → C1 < Z0 .

Finally we show (4.20)

as t → ∞.

k∇ω(t)k → 0

The vorticity ω = ∇ × U satisfies (4.21)

∂t ω = ν∆ω + ∂1 Θ.

Taking the inner product of (4.21) with ∆ω yields Z

∆ω ∂t ω dx = νk∆ωk2 +

Ω

Due to ω = 0 on ∂Ω,

Z

∂1 Θ∆ω dx.

Ω

1 d k∇ωk2 + νk∆ωk2 ≤ k∇Θkk∆ωk. 2 dt

By Young’s inequality, 1 d ν 1 k∇ωk2 + νk∆ωk2 ≤ k∇Θkk∆ωk ≤ k∆ωk2 + Z0 . 2 dt 2 2ν

Integrating in time yields, for any 0 ≤ s < t < ∞, (4.22)

k∇ω(t)k2 − k∇ω(s)k2 ≤ 1 |t − s|Z0 . ν

Due to the divergence-free condition ∇ · U = 0, ∂1 ω = ∆U2 and ∂2 ω = −∆U1 . Thus, according to (4.18), Z

0

or

∞

k∇ω(t)k2 dt =

Z

0

∞

k∆Uk2 dt < ∞

k∇ω(t)k2 ∈ L1 (0, ∞).

Lemma 3.1 then implies (4.20). This completes the proof of linear stability.

22

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

4.3. Linear instability for T 0 (y) = α < 0. Finally we assume T 0 (y) = α < 0 and prove the instability of any periodic perturbation. If the perturbation is periodic, namely (U, Θ) is spatially periodic, then the instability can be easily seen from the Fourier transform of (1.13). In fact, taking the Fourier transform of (1.13), we obtain 

  b (k, t) −ν|k|2 U  b    V (k, t)  =  0 b t) Θ(k, 0

0 −ν|k|2 −α

1 k2 − k|k| 2

k12 |k|2

0



 b (k, t) U  b    V (k, t)  b t) Θ(k,

When α < 0, the coefficient matrix has a positive eigenvalue. In fact, the characteristic polynomial is given by p(λ) = λ + ν|k|2

and the eigenvalues of the coefficient matrix






λ2 + ν|k|2 λ + α

λ1 = −ν|k|2 , 1 1 λ2 = − ν|k|2 − 2 2 1 1 λ3 = − ν|k|2 + 2 2

s s

k12 |k|2



ν 2 |k|4 − 4α

k12 , |k|2

ν 2 |k|4 − 4α

k12 . |k|2

Clearly, for α < 0, λ1 < 0, λ2 < 0 and λ3 > 0. This implies any periodic perturbation is linearly unstable. This completes the proof of Theorem 1.4. 5. N UMERICAL S IMULATIONS This section is devoted to the numerical illustration of the analytical results recorded in Theorems 1.1 through Theorem 1.4, and further investigation of some unresolved problems associated with (1.4). For numerical simplicity, we simulate the partially dissipative system subject to the stress free boundary conditions on Ω = [−1, 1] × [−1, 1]. The main purpose of this section is try to understand the explicit decay rate of the velocity field and long-time dynamics of the temperature function associated with (1.4). 5.1. Numerical Method. We use a projection finite element method for u and P (see [10] for details). The projection method in the strong form is Du∗,k+1 − ν∆u∗,k+1 ∆t u∗,k+1 P k+1

= −∇P Ex,k+1 + (0, θ)Ex,k+1 − (u · ∇u)Ex,k+1 , = uk+1 + ∆t∇φ, = P Ex,k+1 + 1.5φ − ν∇ · u∗,k+1 ,

where Du∗,k+1 1.5u∗,k+1 − 2uk + 0.5uk−1 = ∆t ∆t k+1 is the BDF2 temporal discretization of ∂u at time t , and ∂t (u · ∇u)Ex,k+1 = 2(u · ∇u)k − (u · ∇u)k−1

is a second order extrapolation of (u · ∇u) at time tk+1 . Similarly, P Ex,k+1 = 2P k − P k−1 and θEx,k+1 = 2θk − θk−1 . We use the P2/P1 finite element method to solve for u and P . The weak form is as follows. At time step tk+1 , we seek uk+1 ∈ (H 1 (Ω))2 with uk+1 · n|∂Ω = 0 and P k+1 ∈ L2 (Ω), such that for any

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 23

v ∈ (H 1 (Ω))2 with v · n|∂Ω = 0, ψ ∈ H 1 (Ω) and Q ∈ L2 (Ω), it holds that   Du∗,k+1 , v + ν(∇u∗,k+1 , ∇v) = (P Ex,k+1 , ∇ · v) + ((0, θ)Ex,k+1 , v) − ((u · ∇u)Ex,k+1 , v), ∆t   ∇ · u∗,k+1 (∇φ, ∇ψ) = − ,ψ , ∆t (P k+1 , Q)

=

(uk+1 , v)

=

(P Ex,k+1 + 1.5φ − ν∇ · u∗,k+1 , Q),

(u∗,k+1 − ∆t∇φ, v).

Finally, the transport equation of θ is solved by a third order accurate WENO scheme [51, 53]. The numerical tests show this scheme is second order accurate for velocity, pressure, and temperature in L∞ norm (data not shown). For the long-time numerical simulation for system (1.4), we have used the resolutions 50 × 50, 100 × 100, and 200 × 200. When the time exceeds 100730, all these simulations exhibit some instabilities whose results might be not reliable anymore. But before this time, we observe convergence when the mesh is refined. What we present here are the results from 200 × 200 resolution before time = 100730. 5.2. Decay of Velocity Field. In this subsection, we numerically illustrate the analytical results obtained in Theorem 1.2, regarding the decay of the velocity field associated with (1.4). Numerical results are presented in Figure 1. ||(U,V)|| H1 ||(U,V)||

0

10

||(U,V)|| ||(U,V)t||

L∞ L2 L2

−1.0

2.0xt

1.1xt−1.1

−5

10

1.0xt−1.2

−10

10

0

5

10

10 t

F IGURE 1. Decay of velocity field of numerical solution of (1.4) on [−1, 1] × [−1, 1] with boundary conditions u · n|∂Ω = 0 and ω|∂Ω = 0, initial data u0 (x, y) = sin(2πx) cos(2πy), v0 (x, y) = − cos(2πx) sin(2πy), θ0 (x, y) = sin(3x) cos(2y) + (x − 0.5)3 + 1/(y + 10), and ν = 1. This figure plots the evolution of various norms of the velocity for time up to t = 100730. The dashed lines are the least-squares fittings to the corresponding curves of the same color. In the simulation of (1.4), we choose u0 (x, y) = sin(2πx) cos(2πy), v0 (x, y) = − cos(2πx) sin(2πy), θ0 = sin(3x) cos(2y) + (x − 0.5)3 + 1/(y + 10), and ν = 1. By Theorem 1.2, the kinetic energy and the first order derivatives of the velocity field converge to zero. Fig. 1 plots the evolution of various norms of the solution for time up to t = 100730, which illustrates that ku(t)kH 1 , k∂t u(t)k and ku(t)kL∞ asymptotically approach zero. More interestingly, from Fig. 1 we see that although the decay of the velocity field exhibits a slightly oscillatory fashion, the bulk of the curves presents algebraic decay rates by comparing with polynomial functions, which are slower than exponential decay rates.

24

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

5.3. Evolution of Pressure and Temperature. In this subsection, we numerically illustrate the time evolution of the pressure and temperature functions. Numerical results are presented in Figures 2 and 3.

F IGURE 2. Numerical solutions of the pressure function of (1.4) subject to the stress free boundary conditions on [−1, 1] × [−1, 1] at several time steps.

F IGURE 3. Numerical solutions of the temperature function of (1.4) subject to the stress free boundary conditions on [−1, 1] × [−1, 1] at several time steps. Fig. 2 and Fig. 3 show that the pressure and temperature become horizontally homogeneous and stratify in the vertical direction as time evolves – a scenario that is consistent with the analytical results recorded in Theorem 1.2, especially the third point of (1.8). The numerical results in Fig. 3 show three stages for the heat to get to the stratification. The first is the interior bulk convection as shown from time

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 25

= 0 to time = 20, where the bulk red region of high temperature and blue region of low temperature are swept to the upper and lower parts of the domain, respectively. The second is boundary layer eruption. One relatively low temperature layer is held on the top wall (Fig. 3, time = 20) and the slip boundary condition (u · n = 0) is unable to move this layer down to the corresponding positions in the desired stratification. However, driven by gravity, the heat erupts from some separate boundary spots into the interior domain whose pattern exhibits the Rayleigh-Taylor instability (Fig. 3, from time = 2230 to 10730). This stage actually proceeds from time = 1600 to time = 40000 in the numerical simulation. Finally, after all the eruptions are emitted, the whole system enters a relaxation stage. On account of the complication produced by this instability, we remark that the identification of the exact ansatz of the pressure and temperature functions is still elusive. R To look into the stratification in more detail, we define T (y, t) = Ω θ(x, y, t)dx/2, the average value of temperature at a y -layer. We plotted the convergence behavior of θ(x, y, t) to T (y, t) in Fig. 4 [a, b, c]. The L2 norm shows almost monotone decreasing but the L∞ norm shows overall very slow decreasing with oscillations. The profile of T (y, t) shows monotonic increasing in y -direction at time = 105 (Fig. 4 [b, c]). To study the convergence of the mechanical energy, we take the field T (y, t) at t = 120830 as the reference state. Fig. 4[d] shows the difference of the total mechanical energy (kinetic + potential) at time ≤ 105 and that of the reference state. This indicates a power-rule decreasing of the energy. ||θ(x,y,t)−T(y,t)|| in L∞ and L2 norms L∞ L2

0

10

64.8t−0.8 −1

10

−2

10

−3

10

3

10

4

10

time

5

10

[a]

[b] ∫Ω((u2+v2)/2−yθ)dxdy −2

mechanical energy 79029t−2.0

10

−3

10

−4

10

−5

10

[c]

−6

10

3

10

4

10

5

10

[d]

F IGURE 4. Numerical results. [a] Convergence of θ(x, y, t) to T (y, t) over time. [b,c] Plots of the function T (y, t) from two viewpoints. [d] Difference of total mechanical energy before time ≤ 105 with that at time = 120830. 6. C ONCLUSION AND L OOKING A HEAD We have studied the global well-posedness and long-time asymptotic behavior of large-amplitude classical solutions to the 2D Boussinesq equations without thermal diffusion on bounded domains with non-smooth boundaries and subject to the stress-free velocity boundary conditions. Utilizing energy methods we showed that for initial data with low regularity, there exist unique global-in-time solutions

26

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

to the initial-boundary value problems (IBVPs) of the model, and the L2 norms of the velocity field and its first order spatial and temporal derivatives converge to zero as t → ∞. Consequently we found that the pressure and temperature functions stratify in the vertical direction in a weak topology. Moreover, we established the linear stability of the hydrostatic equilibrium T (y) = αy + θ¯ when α > 0, and the instability of periodic perturbations of the ansatz when α < 0. Our results indicate how buoyancy plays a dominant role when thermal diffusion is negligible, and partially demonstrate the effectiveness of the 2D Boussinesq equations for modeling stratification-dominated situations in the applied sciences. In addition, the analytical approach developed for proving the long-time behavior results is of independent interest and may be adopted to study other PDE systems with similar partially dissipative structures. Finally we mention that several fundamental questions regarding the long-time behavior of largeamplitude classical solutions to the partially dissipative system (1.4) remain open. We list three of them: • Explicit Decay Rates. An important piece of information that is missing from our analytical results is the explicit decay rate of the velocity field or the total mechanical energy. This is mainly due to the lack of dissipation in the temperature equation. Unlike the fully dissipative free energy formulation (1.2) which generates exponential decay of the perturbation, the Lyapunov functional associated with the partially dissipative system (1.15) is not capable of bringing about any explicit decay rate of the perturbation. However, our numerical simulation indicate that the velocity field might converge to zero as a power law which provides a strong hint for a rigorous proof of the phenomenon. We leave this investigation for the future. • Thermal Structure and Stability of the Final State. Another piece of information missing from Theorem 1.2 is a precise description of the final buoyancy distribution in case of general initial conditions. In fact the buoyancy field is just rearranged by the flow so the area of the domain where θ(x, t) ranges between any two values is conserved. This suggests that the final state of the relaxation problem studied here should generically be the unique stably stratified distribution compatible with all these constraints. ˆ : in a rectangular domain Ω = [0, W ] × We provide a formula of the stably stratified state θ(y) [0, H], for any initial value θ0 (x, y), define a height function of temperature value a as (6.1)

h(a) =

1 area{(x, y) ∈ Ω : θ0 (x, y) ≤ a}. W

For a general domain Ω whose horizontal width is given by w(y) at a height y , this formula R h(a) becomes 0 w(y)dy = area{(x, y) ∈ Ω : θ0 (x, y) ≤ a}. The final stably stratified state stratifiˆ cation function, call it θ(y), is then the inverse function of h(a): ˆ = S(y) , h−1 (y). θ(y)

The reasoning behind these formulas is as follows. First, because of mass conservation, h(a) must be a constant over time for any a. Second, we speculate that the ultimate final temperature profile must be nondecreasing in the final state. Note the formula (6.1) satisfies this condition. It remains to be shown that any other stratified states compatible with the initial data are linearly unstable, which does not rule out their emergence from properly prepared initial data but suggests a unique final state evolving from generic conditions. For the numerical test presented in the last section, we computed the final state S(y) numerically and tested the convergence of θ(x, y, t) to it. The results are shown in Figure 5. Note the L2 norm of θ(x, y, t) − S(y) decays as a power rule again and the mean value T (y, t) agrees well ˆ = S(y) with S(y), which is a nonlinear function of y , when time is large. This demonstrates θ(y) is indeed the thermal structure of the final state.

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 27

Comparison between T(y,t) and S(y)

||θ(x,y,t)−S(y)|| in L∞ and L2 norms

2

0

10

1

T(y,t) at t=100730 S(y)

0

−1

10

−1

∞

L L2

−2

10

−3

10

−2 −3

5.7t−0.5 3

10

4

10

time

5

10

[a]

−4 0

0.5

1

1.5

height from domain bottom

2

[b]

F IGURE 5. [a] Convergence of θ(x, y, t) to S(y). The fitting of L2 error uses a power function of 5.7t−0.5 with least squares method. [b] Comparison between T (y, t) at t = 100730 and S(y). R Here, T (y, t) = Ω θ(x, y, t)dx/2 and Ω = [0, 1]2 .

• Infinite Prandtl Number Convection. The ratio of a fluid’s kinematic viscosity ν to its thermal diffusion coefficient κ is the Prandtl number P r = ν/κ so the κ = 0 model studied here corresponds to an infinite Prandtl number limit of the Boussinesq system. Traditionally (see, e.g., [58]) the time t is rescaled to the dimensionless variable κt/`2 where ` is an appropriate domain length scale (e.g., the height of a rectangular domain or the thickness of a spherical shell) and velocities are measured in units of κ/`. Then, suitably rescaling buoyancy and pressure, the formal P r → ∞ limit in a stress-free thermally insulating container is  ∇P = ∆u + θe2 ,      ∂t θ + u · ∇θ = ∆θ,    ∇ · u = 0,

 θ(x, 0) = θ0 (x),      u · n|∂Ω = 0, ω|∂Ω = 0,     n · ∇θ|∂Ω = 0.

This is a very different system than (1.4), ostensibly describing dynamics on even longer O(`2 /κ) times scales than those captured by the strictly non-diffusive model: diffusion is still

effective for the temperature on these time scales while the flow is totally viscous dominated. Boundary conditions are still required for the temperature field, but no initial data is required for the velocity field. The dynamics modeled by the system (1.4) studied in this paper, corresponding to measuring time on viscous (O(`2 /ν)) time scales, would be considered an “initial layer” problem for these infinite Prandtl number equations of motion. √ Alternatively we could consider dynamics on an intermediate, say O(`2 / νκ), time scale in which case the formal P r → ∞ limit for appropriately rescaled dependent variables is  ∇P = ∆u + θe2 ,       ∂ θ + u · ∇θ = 0,   t ∇ · u = 0,     θ(x, 0) = θ0 (x),    u · n| = 0, ω| = 0 ∂Ω ∂Ω

28

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

requiring neither boundary conditions for the temperature nor initial data for the flow. It is natural to wonder to what extend the long time behavior of the relaxation on viscous time scales might be quantitatively captured by an “over damped” model like this—which keeps the conserved quantities in play for the temperature—keeping in mind that the initial data for the temperature immediately above is not simply proportional to the initial temperature distribution in (1.4). Rather, θ0 (x) for this system would correspond to the sort of rapidly-rearranged temperature distribution seen in the simulations (like that seen around time 20 in Figure 3). For our purposes the real value of such a reduced model might be to facilitate accurate estimation of the final state and observed kinetic and/or total energy relaxation kinetics. Acknowledgements. Support for this work came in part from NSF Awards DMS-1515161 (CRD) and DMS-1614246 (JW), a John Simon Guggenheim Foundation Fellowship (CRD), a Simons Foundation Collaboration Grant for Mathematicians 413028 (KZ), and from Louisiana Board of Regents Research Competitiveness Subprogram Award LEQSF(2015-18)-RD-A-24 (KZ). KZ also gratefully acknowledges start up funding from the Department of Mathematics at Tulane University. We thank Professor Chi-Wang Shu for helpful discussions regarding numerical methods. R EFERENCES [1] H. Abidi and T. Hmidi, On the global well-posedness for Boussinesq system, J. Diff. Equ. 233 (2007): 199–220. [2] D. Adhikari, C. Cao, H. Shang, J. Wu, X. Xu and Z. Ye, Global regularity results for the 2D Boussinesq equations with partial dissipation, J. Diff. Equ. 260 (2016): 1893–1917. [3] D. Adhikari, C. Cao and J. Wu, The 2D Boussinesq equations with vertical viscosity and vertical diffusivity, J. Diff. Equ. 249 (2010): 1078–1088. [4] D. Adhikari, C. Cao and J. Wu, Global regularity results for the 2D Boussinesq equations with vertical dissipation, J. Diff. Equ. 251 (2011): 1637–1655. [5] D. Adhikari, C. Cao, J. Wu and X. Xu, Small global solutions to the damped two-dimensional Boussinesq equations, J. Diff. Equ. 256 (2014): 3594–3613. [6] C. Bardos, F. Di Plinio and R. Temam, The Euler equations in planar nonsmooth convex domains, J. Math. Anal. Appl. 407 (2013): 69-89. [7] A. Biswas, C. Foias and A. Larios, On the attractor for the semi-dissipative Boussinesq equations, Ann. Inst. H. Poincare Anal. Non Lineaire 34 (2017): 381–405. [8] J.P. Bourguignon and H. Brezis, Remarks on the Euler equation, J. Funct. Anal. 15 (1975): 341–363. [9] L. Brandolese and M. Schonbek, Large time decay and growth for solutions of a viscous Boussinesq system, Transactions AMS 364 (2012): 5057–5090. [10] D.L. Brown, R. Cortez and M.L. Minion, Accurate projection methods for the incompressible Navier-Stokes equations, J. Comput. Phys. 168 (2001): 464–499. [11] J.R. Cannon, and E. DiBenedetto, The initial value problem for the Boussinesq equations with data in Lp , Approximation methods for Navier-Stokes problems (Proc. Sympos., Univ. Paderborn, Paderborn, 1979), Lecture Notes in Math. 771, Springer, Berlin, 1980, pp. 129–144, [12] C. Cao and J. Wu, Global regularity for the 2D anisotropic Boussinesq equations with vertical dissipation, Arch. Rational Mech. Anal. 208 (2013): 985–1004. [13] D. Chae, Global regularity for the 2D Boussinesq equations with partial viscosity terms, Adv. Math. 203 (2006): 497–513. [14] D. Chae, P. Constantin and J. Wu, An incompressible 2D didactic model with singularity and explicit solutions of the 2D Boussinesq equations, J. Math. Fluid Mech. 16 (2014): 473–480. [15] D. Chae and H. Nam, Local existence and blow-up criterion for the Boussinesq equations, Proc. Roy. Soc. Edinburgh Sect. A 127 (1997):935–946. [16] D. Chae, S. Kim, and H. Nam, Local existence and blow-up criterion of H¨older continuous solutions of the Boussinesq equations, Nagoya Math. J. 155 (1999): 55–80. [17] D. Chae and O.Y. Imanuvilov, Generic solvability of the axisymmetric 3-D Euler equations and the 2-D Boussinesq equations, J. Diff. Equ. 156 (1999): 1–17.

LONG TIME BEHAVIOR OF THE TWO-DIMENSIONAL BOUSSINESQ EQUATIONS WITHOUT BUOYANCY DIFFUSION 29

[18] D. Chae and J. Wu, The 2D Boussinesq equations with logarithmically supercritical velocities, Adv. Math. 230 (2012): 1618–1645. [19] D. C´ordoba, C. Fefferman and R. De La LLave, On squirt singularities in hydrodynamics, SIAM J. Math. Anal. 36 (2004): 204–213. [20] R. Danchin and M. Paicu, Existence and uniqueness results for the Boussinesq system with data in Lorentz spaces, Phys. D 237 (2008): 1444–1460. [21] R. Danchin and M. Paicu, Global well-posedness issues for the inviscid Boussinesq system with Yudovich’s type data, Comm. Math. Phys. 290 (2009): 1–14. [22] R. Danchin and M. Paicu, Global existence results for the anisotropic Boussinesq system in dimension two, Math. Models Methods Appl. Sci. 21 (2011): 421–457. [23] W. E and C.-W. Shu, Small-scale structures in Boussinesq convection, Phys. Fluids 6 (1994): 49–58. [24] J. Geng and Z. Shen, The Neumann problem and Helmholtz decomposition in convex domains, J. Funct. Anal. 259 (2010): 2147-2164. [25] P. Grisvard, Elliptic problems in nonsmooth domains, Monographs and Studies in Mathematics, 24, Pitman (Advanced Publishing Program), Boston, MA, 1985. [26] L. He, Smoothing estimates of 2d incompressible NavierStokes equations in bounded domains with applications, J. Funct. Anal. 262 (2012): 3430-3464. [27] T. Hmidi and S. Keraani, On the global well-posedness of the 2D Boussinesq system with a zero diffusivity, Adv. Diff. Equ. 12 (2007): 461–480. [28] T. Hmidi and S. Keraani, On the global well-posedness of the Boussinesq system with zero viscosity, Indiana Univ. Math. J. 58 (2009): 1591–1618. [29] T. Hmidi and S. Keraani and F. Rousset, Global well-posedness for a Boussinesq-Navier-Stokes system with critical dissipation, J. Differ. Equ. 249 (2010): 2147–2174. [30] T. Hmidi and S. Keraani and F. Rousset, Global well-posedness for Euler-Boussinesq system with critical dissipation, Comm. Partial Differ. Equ. 36 (2011): 420–445. [31] J. Holton, An Introduction to Dynamic Meteorology, International Geophysics Series, 4th edition, Elsevier Academic Press, 2004. [32] T. Hou and C. Li, Global well-posedness of the viscous Boussinesq equations, Disc. Cont. Dyn. Sys. 12 (2005): 1–12. [33] L. Hu and H. Jian, Blow-up criterion for 2-D Boussinesq equations in bounded domain, Front. Math. China 2 (2007): 559–581. [34] W. Hu, I. Kukavica, M. Ziane, Persistence of regularity for a viscous Boussinesq equations with zero diffusivity, Asymptot. Anal. 91(2) (2015): 111–124. [35] W. Hu, I. Kukavica and M. Ziane, On the regularity for the Boussinesq equations in a bounded domain, J. Math. Phys. 54 (2013), 081507, 10 pp. [36] D. Jerison and C.E. Kenig, The inhomogeneous Dirichlet problem in Lipschitz domains, J. Funct. Anal. 130 (1995): 161-219. [37] Q. Jiu, C. Miao, J. Wu, and Z. Zhang, The 2D incompressible Boussinesq equations with general critical dissipation, SIAM J. Math. Anal. 46 (2014): 3426–3454. [38] Q. Jiu, J. Wu, and W. Yang, Eventual regularity of the two-dimensional Boussinesq equations with supercritical dissipation, J. Nonlinear Science 25 (2015): 37–58. [39] C. Lacave, E. Miot, and C. Wang, Uniqueness for the 2-D Euler equations on domains with corners, Indiana Univ. Math. J. 63 (2014): 1725-1756. [40] M. Lai, R. Pan and K. Zhao, Initial boundary value problem for 2D viscous Boussinesq equations, Arch. Rational Mech. Analysis 199 (2011): 739–760. [41] A. Larios, E. Lunasin and E.S. Titi, Global well-posedness for the 2D Boussinesq system with anisotropic viscosity and without heat diffusion, J. Diff. Equ. 255 (2013): 2636–2654. [42] D. Li and X. Xu, Global wellposedness of an inviscid 2D Boussinesq system with nonlinear thermal diffusivity, Dyn. Par. Diff. Equ. 10 (2013): 255–265. [43] J. Li, H. Shang, J. Wu, X. Xu and Z. Ye, Regularity criteria for the 2D Boussinesq equations with supercritical dissipation, Comm. Math. Sci. 14 (2016): 1999–2022. [44] J. Li and E.S. Titi, Global well-posedness of the 2D Boussinesq equations with vertical dissipation, Arch. Ration. Mech. Anal. 220 (2016): 983-1001. [45] S.A. Lorca and J.L. Boldrini, The initial value problem for a generalized Boussinesq model, Nonlinear Analysis 36 (1999): 457–480.

30

CHARLES R. DOERING, JIAHONG WU, KUN ZHAO, AND XIAOMING ZHENG

[46] A. Majda, Introduction to PDEs and Waves for the Atmosphere and Ocean, Courant Lecture Notes in Mathematics, no. 9, AMS/CIMS, 2003. [47] A. Majda and A. Bertozzi, Vorticity and Incompressible Flow, Cambridge University Press, Cambridge, UK, 2002. [48] J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1987. [49] R. Salmon, Lectures on Geophysical Fluid Dynamics, Oxford Univ. Press, 1998. [50] A. Sarria and J. Wu, Blowup in stagnation-point form solutions of the inviscid 2d Boussinesq equations, J. Diff. Equ. 259 (2015): 3559–3576. [51] C.-W. Shu, High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD, Inter. J. Comput. Fluid Dyn., 17 (2003): 107–118. [52] A. Stefanov and J. Wu, A global regularity result for the 2D Boussinesq equations with critical dissipation, J. d’Analyse Mathematique, to appear. [53] S. Tan and C.-W. Shu, Inverse Lax-Wendroff Procedure for Numerical Boundary Conditions of Conservation Laws, J. Comput. Phys. 229 (2010): 814–8166. [54] Y. Taniuchi, A note on the blow-up criterion for the inviscid 2-D Boussinesq equations, in: R. Salvi (Ed.), The NavierStokes Equations: Theory and Numerical Methods, Lec. Notes Pure Appl. Math. 223 (2002): 131–140. [55] R. Temam, Navier-Stokes Equations, North Holland, 1977. [56] G.K. Vallis, Atmospheric and Oceanic Fluid Dynamics, Cambridge University Press, Cambridge, U.K., 2006. [57] W. von Wahl, Estimating ∇u by div u and curl u, Math. Methods Appl. Sci. 15 (1992): 123–143. [58] X. Wang and J.P. Whitehead, A bound on the vertical transport of heat in the ‘ultimate’ state of slippery convection at large Prandtl numbers, J. Fluid Mech. 729 (2013): 103–122. [59] J.P. Whitehead and C.R. Doering, Ultimate State of Two-Dimensional Rayleigh-B´enard Convection between Free-Slip Fixed-Temperature Boundaries, Phys. Rev. Lett. 106 (2011): 244501. [60] J. Wu and X. Xu, Well-posedness and inviscid limits of the Boussinesq equations with fractional Laplacian dissipation, Nonlinearity 27 (2014): 2215–2232. [61] J. Wu, X. Xu, L. Xue and Z. Ye, Regularity results for the 2D Boussinesq equations with critical and supercritical dissipation, Comm. Math. Sci. 14 (2016): 1963–1997. [62] W. Yang, Q. Jiu and J. Wu, Global well-posedness for a class of 2D Boussinesq systems with fractional dissipation, J. Differential Equations 257 (2014): 4188–4213. [63] Z. Yoshida and Y. Giga, Remarks on spectra of operator rot, Math. Z. 204 (1990): 235–245. [64] K. Zhao, 2D inviscid heat conductive Boussinesq system in a bounded domain, Michigan Math. J. 59 (2010): 329–352. (C. R. Doering) D EPARTMENT OF M ATHEMATICS , U NIVERSITY OF M ICHIGAN , A NN A RBOR , MI 48109 E-mail address: [email protected] (J. Wu) D EPARTMENT OF M ATHEMATICS , O KLAHOMA S TATE U NIVERSITY, S TILLWATER , OK 74078 E-mail address: [email protected] (K. Zhao) D EPARTMENT OF M ATHEMATICS , T ULANE U NIVERSITY, N EW O RLEANS , LA 70118 E-mail address: [email protected] (X. Zheng) D EPARTMENT OF M ATHEMATICS , C ENTRAL M ICHIGAN U NIVERSITY, M OUNT P LEASANT, MI 48859 E-mail address: [email protected]

HIGHLIGHTS • We study the 2D Boussinesq equations without buoyancy diffusion. • We establish the global-in-time existence and uniqueness of classical solutions. • We study linear stability and instability of some stratified hydrostatic equilibria. • Direct numerical simulations corroborate analytical results.