Well-posedness of strong solutions to the anelastic equations of stratified viscous flows

We establish the local and global well-posedness of strong solutions to the two- and three-dimensional anelastic equations of stratified viscous flows. In this model, the interaction of the density profile with the velocity field is taken into account, and the density background profile is permitted to have physical vacuum singularity. The existing time of the solutions is infinite in two dimensions, with general initial data, and in three dimensions with small initial data.


Introduction
The anelastic Navier-Stokes system for stratified flows, ρ(∂ t u + u · ∇u) + ρ∇p = ∆u in Ω, div(ρu) = 0 in Ω, is derived as the limiting system of the compressible Navier-Stokes system after filtering out the acoustic waves for strong stratified flows. Here the velocity field u and the pressure p are the unknowns while the background density ρ is given as a time-independent, negative function. The rigorous derivation of (1) can be found in [21]. Comparing to the incompressible Navier-Stokes system (see, e.g., [5,26]), the main difference is the incompressible condition divu = 0 is replaced by the anelastic relation div(ρu) = 0 with the background density profile ρ, which represents the strong stratification owing to the balance of the gravity and the pressure (see, e.g., [10]). Such an approximation preserves slight compressibility while filtering out the acoustic waves, which significantly simplifies the original compressible Navier-Stokes system, and enables more efficient computation applications to relevant model flows in physical reality. In particular, the anelastic approximation is used to describe the semi-compressible ocean dynamics (see, e.g., [7,8]), as well as the tornado-hurricane dynamics (see, e.g., [22,24]). We refer the readers to [1, 2, 9, 14-16, 20, 23] for related topics and comparisons of various models of the atmospheric and oceanic dynamics. We remark that the background density profile ρ in the anelastic relation div(ρu) = 0 is given by the resting state ∇P (ρ) = ρg e z , where P (ρ) denotes the pressure potential and g is the gravity acceleration. For the sake of simplifying the presentation, we have choosed the gravity to point upwards, which can be done after performing a vertical reflection of the coordinates. In the case when the flow connects to vacuum continuously, the resting state yields a degenerate density profile. For an isentropic flow with P (ρ) = ρ γ , γ > 1, this implies ρ γ−1 ≃ z, referred to as the physical vacuum in the study of compressible flows (see, e.g., [13,17]). The main characteristics of the physical vacuum is the Hölder continuity of the background density profile, whose derivatives are singular at z = 0. While there are some recent developments in the global stability of background solutions to compressible Euler or Navier-Stokes equations for one-dimensional or radial-symmetric flows (see, e.g., [11,12,18,19]), the corresponding multi-dimensional problem is mostly open. On the other hand, after formally filtering out the acoustic waves by sending the Mach number and the Froude number to zero at the same rate in the compressible Navier-Stokes equations with physical vacuum, the resulting equations appear to be the aforementioned anelastic system with ρ = z α , α = 1/(γ − 1) > 0.
Compatibility conditions for u in are given by, ∂ z v in z=0,1 , w in z=0,1 = 0, div(ρu in ) = 0, ρu t t=0 = ρu in,1 := ∆u in − ρu in · ∇u in − ρ∇p in ∈ L 2 (Ω), (6) where p in is the solution to the following elliptic problem div(ρ∇p in ) = div∆u in − div(ρu in · ∇u in ), Hereafter, for any subscript s, we will alway use the letter v s to denote the horizontal velocity, w s to denote the vertical velocity, and u s to denote the velocity field, i.e., u s = (v s , w s ) ⊤ .
Remark 1. The meaning of (2) is that the non-degenerate profile ρ can be extended as a non-degenerate and smooth function in 2T n . On the other hand, one can replace ρ pv in (3) with any density profile which satisfies the aforementioned physical vacuum near z = 0 (i.e., ρ ≃ z α near z = 0), and is smooth and non-degenerate at z = 1. The requirement of smoothness and the non-degeneracy of the density profile at z = 1 is owing to technical reasons.

Remark 2.
We require the supports of v in and w in to be away from the boundary {z = 0, 1}. This can be modified as follows: let η, ζ be even, odd in the z-variable, respectively, and (η, ζ) ∈ H 2 (2T n ; R n−1 ) × H 2 (2T n ; R); Equation (7) should be considered as a compatibility condition of u in .
In comparison to the Navier-Stokes system (see, e.g., [5]), the density profile interacts with the velocity field. To explain this statement, let us forget about the boundary and consider (1) in 2T n for a moment. Also we assume the density profile ρ is non-degenerate and smooth. Let u be any smooth vector field, and suppose that it can be decomposed, in analogy with the Helmholtz decomposition, as where div(ρu ρ,σ ) = 0, and φ is a smooth function. One can see that u ρ,σ and ∇φ are orthogonal in L 2 (ρ d x), i.e., the square-integrable space with respect to the measure ρ d x, where d x is the Lebesgue measure. Then one can easily see that, the H s -regularity of u ρ,σ , s ≥ 1, depends not only on the regularity of u, but also on that of ρ. Such an interaction of ρ in the H s -regularity estimates causes the main difficulty in the study of strong solutions. As one will see later, one will need to consider the interaction of the pressure term ρ∇p with the nonlinearity term ρu · ∇u as well as the viscosity term ∆u. We take advantage of the regular density profile in (2), and consider an associated problem in 2T n , whose solutions satisfies (1) in Ω. Then we employ an elementary approach in the Galerkin's approximation which takes into account the aforementioned interactions. After studying the regularity, we restrict our solutions back to the original domain Ω and obtain a unique strong solution to (1) with (4) in the non-degenerate case, i.e., (2).
To deal with the physical vacuum profile in (3), we approximate the problem with a sequence of non-degenerate profiles in the class of (2). The existence theorem in the non-degenerate case yields a sequence of approximating solutions. Then we derive the necessary uniform weighted estimates. To handle the physical vacuum density profile, the desired strong solutions to (1) in the physical vacuum case, i.e., (3), are constructed as the limit of the approximating sequence. However, the solutions that we obtain lack regularity on the boundary {z = 0}, due to the weighted estimates. In particular, the solutions are not regular enough to have trace of ∇u on {z = 0}, which causes troubles when one try to show the uniqueness of solutions. We employ the arguments originated in [25] for the Navier-Stokes system to establish the uniqueness of strong solutions.
Next, we sum up the main theorems. The first theorem concerns the local well-posedness of strong solutions to (1): Theorem 1. Let ρ satisfy either (2) or (3). Consider initial data u in ∈ H 2 satisfying (5) and (6). There exists a unique strong solution (u, p) to (1) with (4) in [0, T ], for some T ∈ (0, ∞). In the case of (2), the strong solution satisfies the regularity: In the case of (3), the strong solution satisfies the regularity: See section 2 for the notations ∂ zz , ∇ h etc..
We refer the detailed description of local well-posedness to Theorem 3 and Theorem 4, below. At the same time, we also have the following theorem concerning global well-posedness of strong solutions: Theorem 2. Under either one of the following conditions, the existing time of the local strong solutions constructed in Theorem 1 becomes infinite: 1. n = 2; 2. n = 3, provided initial velocity u in satisfies, for ρ satisfying either (2) or (3), with some µ ∈ (0, 1), small enough.
The rest of this work is organized as follows. In section 2, we summarize the notations, definitions and inequalities which will be used in this paper. In section 3, we construct the local strong solutions to (1) with the non-degenerate density profile. In section 4, we consider (1) with the physical vacuum profile, when the uniform estimates and approximation arguments are presented. These two sections finish the proof of Theorem 1. In the last two sections, i.e., sections 5 and 6, we employ some global a priori estimates, which lead to the proof of Theorem 2.

Preliminaries
Throughout this paper, we use the following definition of strong solutions: Definition 1 (Strong solutions). (u, p) is called a strong solution to (1) if system (1) holds almost everywhere in Ω.
We use the notation ∂ x to denote the spatial derivative in the horizontal direction, i.e. derivative with respect to x ∈ 2T, when n = 2, and x 1 , x 2 for x = (x 1 , x 2 ) ⊤ ∈ 2T 2 , when n = 3; the notation ∂ z to denote the spatial derivative in the vertical direction; the notation ∂ t to denote the temporal derivative; u s for s ∈ {t, x, z} is short for ∂ s u; also u s1s2 and ∂ s1s2 u for s 1 , s 2 ∈ {t, x, z} are short for ∂ s1 ∂ s2 u. div h , ∇ h , ∆ h are used to denote the divergence, the gradient, the Laplace, respectively, in horizontal direction, i.e. div h = ∂ x , when n = 2, and div h = ∇ h ·, when n = 3, In addition, we abuse the notation: · d x, depending on the context. L p , H k are used to denote L p (Ω), H k (Ω) or L p (2T n ), H k (2T n ), depending on the context. Also, we summarize the symmetric-periodic extensions in the following: Definition 2 (Symmetric-periodic extensions). For any smooth function f defined in Ω := 2T n−1 × (0, 1), one can extend it to an even function in Ω ± := 2T n−1 × (−1, 1), using the even-symmetric extension E + s , defined by In addition, if lim z→0 + f (x, z) = 0, one can also extend f to a odd function in Ω ± , using the odd-symmetric extension E − s , defined by Also, for any smooth function g defined in Ω ± , one can extend it to a function in 2T n using the periodic extension E p , defined by k is the integer such that z − 2k ∈ (−1, 1). (10) Then the even-symmetric-periodic extension operator is defined by The odd-symmetric-periodic extension operator is defined by To study (1) in the case of physical vacuum, i.e., (3), we will need to apply the following Hardy-type inequalities: Lemma 1 (Hardy-type inequalities). Let k = −1 be a real number. Suppose that a function f ∈ C 1 ([0, 1]) satisfies Then for some positive constant C k ∈ (0, ∞), independent of ε ∈ (0, 1), one has 1. for k > −1, In particular, after taking ε = 0 in (13) and (14), one will arrive at the standard Hardy's inequalities.
Our strategy of constructing solutions is: first, we introduce a problem in 2T n , which is associated with (15); then we construct the solutions with enough regularity to the associated problem; by restricting such solutions to the associated problem back in Ω, we obtain the required solutions to (15).
The following theorem is the main part of this section: Theorem 3. Let ρ be a strict positive scalar function in Ω that satisfies (16).
Consider initial data u in ∈ H 2 satisfying (5) and (6). Then there exists a unique strong solution (u, p) to (15), Moreover, the strong solution satisfies the following regularity: Furthermore, the following estimates hold: where C in,ρ ∈ (0, ∞) depends only on the initial data u in and Also, one can choose p to satisfy In addition, let u 1 , u 2 be strong solutions with initial data u 1,in , u 2,in , respectively. Then the following estimate holds, where T * ∈ (0, ∞) is the co-existence time of the solutions, and C in,1,2,T * ∈ (0, ∞) depends on the initial data and T * .
In fact, we will only show the proof of Theorem 3 when n = 2. The case when n = 3 is similar and we omit it for the sake of clarity of our presentation. The proof is done in the following steps: introducing the associated problem via the symmetric-periodic extension; introducing the Galerkin approximating problem; establishing existence of strong solutions; improving the regularity; establishing uniqueness and continuous dependency on the initial data.
Step 0: the associated problem. We observe that system (15) is invariant with respect to the following symmetry: ρ, v, w, p are even, even, odd, even, respectively, with respect to the z-variable. (SYM) Recalling Ω = 2T × (0, 1), that is to say, by extending any solution (ρ, v, w, p) to system (15) to (ρ ± , v ± , w ± , p ± ) satisfies the same equations as in system (15) in domain Ω ± := 2T × (−1, 1). Then the new system for the extended functions (ρ ± , v ± , w ± , p ± ) is invariant with respect to translation z z ± 2. Thus, we further extend (ρ ± , v ± , w ± , p ± ) periodically in the z-variable by applying E p to (ρ ± , v ± , w ± , p ± ). Combining these two extensions together, we obtain, and (ρ * , v * , w * , p * ) satisfies the same equations as in system (15) in domain 2T 2 . Therefore, we end up with the same set of equations as in system (15) in periodic domain 2T 2 , and for simplicity, the same notations ρ, v, w, p are used to denote ρ * , v * , w * , p * . Such a convention will be adopted in the following. Then we have got the following system, with symmetry (SYM). We adopt initial data (E + sp v in , E − sp w in ) ⊤ for (15'), which we will denote by the same notation as the original initial data, i.e., (v in , w in ) ⊤ . Notice, the boundary conditions in (4) are automatically implied by symmetry (SYM). In the next step, a Galerkin approximating procedure will be used to construct solutions to (15').
Step 1: the Galerkin approximating problem. Given any non-negative integer m, we consider the finite dimensional space, denoted by X m and defined as follows: with a v k , a w k , b k being complex-valued scalar functions of t only and satisfying the reality condition, i.e., Notice that, the dimension of X m over R is 3(2m + 1)(m + 1) − 1. Also, we define the lower-m th -frequency projection operator P m , m ≥ 0, as follows.
Then P m projects (v, w, p) with symmetry (SYM) into X m via P m (v, w, p) = (P m v, P m w, P m p), where we have taken Ω p d x = 0.
Consider any non-negative integer m and (v m , w m , p m ) ∈ X m with a v k , a w k , b k given as in (20). To solve the problem (15'), we consider the following system of ODE: To (22), we will need to reformulate (22) into a system of dimension 3(2m+1)(m+1)−1. In fact, we claim that {b k } can be represented as functions of {(a v k (t), a w k (t))} by inverting a linear algebraic system of dimension (2m + 1)(m + 1) − 1, and one can derive a first- Taking ∂ x and ∂ z to (22) 1 and (22) 2 , respectively, and summing the results together yield, using (22) 3 , which is, due to the even symmetry and the strict positivity of ρ, a non-singular linear system with the unknowns {b k } of dimension (2m + 1)(m + 1) − 1. Thus after solving for {b k }, (22) can be written as the following 3(2m + 1)(m + 1) − 1 dimensional system, In particular, (23) 1 and (23) 2 form the 2(2m We remark that, (22) 3 is preserved by the solutions to (23) with compatible initial data, since (23) (23) 3 and substituting the solutions to (23) 1 and (23) 2 , we will have an ODE system of the form with {F k } k∈Zm being locally Lipschitz continuous, in fact quadratic functions, with respect to the arguments. Here we need again that ρ is strictly positive. Then the existence theorem of ODE systems yields that given initial data where Q m = k∈Zm/{(0,0)} q k e πik1x cos(πk 2 z), with q (k1,k2) = q (−k1,k2) , is determined by solving, as above, the non-singular linear algebraic system there exists a solution (v m (t), w m (t), p m (t))| t∈(0,T * m ) ∈ X m to system (23), or equivalently system (22), for some positive constant T * m ∈ (0, ∞). We remark that, as m → ∞, Q m → 0 in H 3 . In fact, owing to the fact div Hence (v m , w m ) t=0 is an approximation of (v in , w in ).
Step 2: existence of strong solutions. In order to pass the limit m → ∞ in (22) to obtain a solution to (15'), we need to establish that the existence time T * m , obtained above, is independent of m. This is done via some uniform-inm estimates. Let T * * m ≥ T * m be the maximal existing time of the solutions (v m , w m ). All the estimates below in this step are done in the time interval [0, T * * m ).
After taking the L 2 -inner product of (22) 1 and (22) 2 with v m and w m , respectively, summing up the resulting equations and applying integration by parts yield, where we have used (22) 3 . Next, we take the L 2 -inner product of (22) 1 and (22) 2 with ∂ t v m and ∂ t w m , respectively. Similarly, after summing up the resulting equations and applying integration by parts, one will have, since ρ has uniform upper bound and strictly positive lower bound, where we have applied the two-dimensional Sobolev embedding inequality. Thus, we have, after applying Young's inequality and (25), In order to estimate ∇ 2 v m , ∇ 2 w m , we rewrite (22) 1 and (22) 2 in the following pressure-viscosity form: which yield Meanwhile, direct calculations show that Since we have, after applying integration by parts, where we need ρ ∈ C 2 (2T 2 ). Therefore, (28) and (29) imply, On the other hand, taking the L 2 -inner product of (23) 3 with −p m yields Therefore, after applying the two-dimensional Sobolev embedding inequality, together with the fact that ρ is strictly positive, we arrive at Then, (30) and (31) imply Consequently, (26) and (32) yield Thus, (33) implies that, there exists T * ∈ (0, T * * m ), independent of m, such that where C in ∈ (0, ∞) depends only on the initial data and inf x∈2T ρ( x), ρ C 2 (2T 2 ) , and T * is independent of m. Then, after passing m → ∞ with a suitable subsequence according to the weak compactness theorem of Sobolev spaces and Aubin's compactness theorem (see, e.g., [26]), we have obtianed such that p m ⇀ p, weakly in L 2 (0, T * ; H 1 ).
Thus it is easy to verify that (u = (v, w), p) is a strong solution to (15') with (35), which satisfies, according to (34), Step 3: improving the regularity. In this step, we establish the regularity of solution (u, p) to (15') via some a priori estimates. In the following, we use C in,ρ ∈ (0, ∞) to denote a generic constant depending only on the initial data and on inf Here we focus with our estimates over [0, T * ]. We emphasize that all the estimates in this step are formal and can be proved rigorously via the Galerkin method.
First, we obtain the time-derivative estimate. After applying a time derivative to (15') 1 , the resulting equation is ρ(∂ t u t + u · ∇u t + u t · ∇u) + ρ∇p t = ∆u t .
Then after taking the L 2 -inner product of (37) with u t , one has The right-hand side of (38) can be estimated as follows: Together with Young's inequality, (38) and (39) imply Thus applying Grönwall's inequality to the above yields, together with (36), where T * is given in step 2. Then, following similar arguments as in (32), one can obtain, sup Next, we will sketch the H 3 estimate of u. First, applying ∂ x to (15') 1 yields, After taking the L 2 -inner product of (42) with u xxx and applying integration by parts, we obtain where we have applied (40), (41), Hölder's and the Sobolev embedding inequalities. Then, after noticing ∆∂ x u = ∂ xxx u + ∂ xzz u and using (42), (43) implies Moreover, since one can also derive What is left is to obtain the estimate of ∇ 2 p L 2 , or equivalently, the estimate of ∆p L 2 . We rewrite (15') 1 , after multiplying it with ρ and applying div to the resulting, as, div(ρ 2 ∇p) = div ρ∆u − ρ 2 (∂ t u + u · ∇u) where we have used (15') 2 and the identity Thus we have which yields, together with (40) and (41), Then (40), (45) and (47) yield Now we can restrict the solution in Ω. It is easy to verify (u| Ω , p| Ω ) is the strong solution to (15) with the boundary condition (4), and it satisfies (17) owing to (40), (41), and (48). In particular, the regularity of p allows us to take the trace of ∇p on the boundary, and it follows from the construction that p satisfies (18).
In particular, for u in,1 = u in,2 , we have u 1 ≡ u 2 and T * 1 = T * 2 . This finishes the proof of Theorem 3 in the case when n = 2. The case when n = 3 follows by similar arguments, employing the three-dimensional Sobolev embedding inequalities.
Remark 5. It is worth stressing that thanks to the uniqueness of solutions in Step 4, all strong solutions to (15) with (4) with the same initial data as described in Theorem 3 should be equal to the one constructed by our extensionrestriction techniques through Step 0 to Step 3.
In addition, we have the estimates where C in depends only on the initial data. Suppose that there are two solutions u 1 , u 2 with initial data u in,1 , u in,2 in time interval [0, T ]. Then for some constant C in,T ∈ (0, ∞) depending on initial data and T .
Remark 6. Property 4 in Lemma 2 implies that Hardy's inequalities in Lemma 1 can be applied with z + ε replaced by q ε .
Lemma 3. For any fixed ε ∈ (0, 1), assume that α > 1 and (u ε , p ε ) is the solution to (53) as mentioned above. There exists a constant T ∈ (0, ∞) independent of ε, such that the following estimates hold: where C in is a constant depending only on initial data.
Proof. Taking the L 2 -inner product of (53) 1 with u ε implies, after substituting (53) 2 and (4), 1 2 In the meantime, the L 2 -inner product of (53) 1 with u ε,t implies, similarly, The right-hand side of (57) can be estimated as follows, where we have used Property 4 in Lemma 2. Notice, after applying the Sobolev embedding inequality and the Hardy-type inequality in Lemma 1, one can derive provided α/2 − 1 > −1/2, or equivalently, α > 1, where we have used Property 4 in Lemma 2.
Then after taking the L 2 -inner product of (60) with u t , the result is Similarly, the right-hand side of (60) can be estimated as follows, Therefore, combining (56), (57), (58), (59), (61) and (62) gives us where we have used Property 4 in Lemma 2 and applied Young's inequality. In particular, the above yields (55) for a short time.
Next, to obtain the estimates of the spatial derivatives of u ε requires a little work. In fact, we shall proceed with the following steps: 1. obtain estimate for the horizontal derivative; 2. obtain estimate for the pressure; 3. obtain estimate for the L 2 -norm of ∂ zz u ε . In conclusion, we will obtain the following: Lemma 4. In addition to the assumptions in Lemma 3, assume that α > 3/2.
In particular, (63) together with (55) yields, where T is the same as in (55), and C in is some constant depending only on initial data and is independent of ε.
Proof. As mentioned above, we establish the proof in three steps.
Step 1: Obtain estimate for the horizontal derivative. Taking the L 2 -inner product of (53) 1 with ∆ h u ε implies Then, applying Hölder's and the Sobolev embedding inequalities to the righthand side of (65) yields that, together with Property 4 in Lemma 2,

Therefore (65) implies
Step 2: Obtain estimate for the pressure. Notice Therefore, after multiplying (53) 1 with q 3α ε and applying div to the resulting equation, we end up with Recall that p ε satisfies (54), and the integration by parts in the following is allowed. Thus, after taking the L 2 -inner product of (68) with −p ε and applying integration by parts in the resultant using the boundary conditions (4) and (54), we arrive at where Now we need to evaluate the right-hand side of (69). Indeed, applying the Hölder and the Sobolev embedding inequalities in (I) yields To estimate (II), notice that from (53) 2 and (4), we have Then after substituting (70) in (II) and applying integration by parts, it follows, Then applying Properties 4 and 5 in Lemma 2 and the Hölder inequality yields, where in the last inequality, we have applied Hardy-type inequality in the vertical direction (see Lemma 1), with α − 2 > −1/2, i.e., α > 3/2. Therefore, (69) implies, for α > 3/2, q 2α ε ∇p ε L 2 q α/2 ε u ε,t L 2 + u ε H 1 + (z + ε) α ∇u ε 1/2 Step 3: Obtain estimate for ∂ zz u. We rewrite (1) 1 as, Then directly, we have where the last term on the right-hand side can be estimated as Notice, Consequently, (66), (72) and (74) yield (63). Now we collect (55) and (63) to finish the proof. Indeed, after applying the Hardy-type inequality in Lemma 1, we have the following inequalities Therefore, together with Property 4 in Lemma 2, (55), (63) and (76) imply the estimates in (64).
5 Global-in-time a priori estimates when n = 2 In this and the following sections, we present some global-in-time a priori estimates of solutions to (51). These arguments can be rigorously justified following the arguments in sections 3 and 4. The estimates of solutions to (15) are similar, and will be omitted.
Notice that, the regularity of u in (78) allows us to take the following actions. Taking the L 2 -inner product of (51) 1 with u implies, as in (56), where E(t) := e + |∇v(t)| 2 d x + ρ pv |u t (t)| 2 d x.
for some constant C in depending only on the initial data. (92) and (93) imply the global well-posedness.
6 Small data global-in-time a priori estimates when n = 3 Similarly, the estimates in (86), (87) and (88) hold. That is, We estimate the nonlinearities on the right-hand side of (87) and (88) as follows, Then, after denoting E(t) := ρ pv |u| 2 d x + |∇u| 2 d x + ρ pv |u t | 2 d x, Thus we have shown the global well-posedness with small initial data.