Global Well-posedness for the Primitive Equations Coupled to Nonlinear Moisture Dynamics with Phase Changes

In this work we study the global solvability of the primitive equations for the atmosphere coupled to moisture dynamics with phase changes for warm clouds, where water is present in the form of water vapor and in the liquid state as cloud water and rain water. This moisture model contains closures for the phase changes condensation and evaporation, as well as the processes of autoconversion of cloud water into rainwater and the collection of cloud water by the falling rain droplets. It has been used by Klein and Majda in \cite{KM} and corresponds to a basic form of the bulk microphysics closure in the spirit of Kessler \cite{Ke} and Grabowski and Smolarkiewicz \cite{GS}. The moisture balances are strongly coupled to the thermodynamic equation via the latent heat associated to the phase changes. In \cite{HKLT} we assumed the velocity field to be given and proved rigorously the global existence and uniqueness of uniformly bounded solutions of the moisture balances coupled to the thermodynamic equation. In this paper we present the solvability of a full moist atmospheric flow model, where the moisture model is coupled to the primitive equations of atmospherical dynamics governing the velocity field. For the derivation of a priori estimates for the velocity field we thereby use the ideas of Cao and Titi \cite{CT}, who succeeded in proving the global solvability of the primitive equations.


Introduction
Moisture and precipitation still cause major uncertainties in numerical weather prediction models and it is our aim here to develop further the rigorous analysis of atmospheric flow models. In a preceding paper [16] we assumed the velocity field to be given and studied the moisture model for water vapor, cloud water and rain water coupled to the thermodynamic equation through the latent heat in the setting of Klein and Majda [19] corresponding to a basic form of a bulk microphysics model in the spirit of Kessler [18] and Grabowski and Smolarkiewicz [14]. In this work we couple the moisture dynamics to the primitive equations of the atmospheric dynamics by taking over the ideas of Cao and Titi [8] for their recent breakthrough on the global solvability of the latter system. Moreover, cases of partial viscosities and diffusions, arising from the asymptotical analysis in [19], will be analyzed in a future work capitalizing on the results by Cao et al. [4,5,6].
A study of a moisture model coupled to the primitive equations has already been carried out by Coti Zelati et al. in [11]. The moisture model there consists of one moisture quantity coupled to temperature and contains only the process of condensation during upward motion, see e.g. [17]. Since the source term there is modeled via a Heavy side function as a switching term between saturated and undersaturated regions, the analysis, however, requires elaborate techniques. Coti Zelati et al. in [3,10,12] therefore used an approach based on differential inclusions and variational techniques, which have then been coupled to the primitive equations in [11].
The moisture model we are analyzing here is physically more refined and consists of three moisture quantities for water vapor, cloud water and rain water and contains besides all the phase changes due to condensation and evaporation also the autoconversion of cloud water to rain water after a certain threshold is reached, as well as the collection of cloud water by the falling rain droplets.
In the remainder of the introduction we first state the moisture model in pressure coordinates, which have the advantage that under the assumption of hydrostatic balance the continuity equation takes the form of the incompressibility condition.
1.1. Governing equations. Solvability of the full geophysical governing equations (without moisture) is a long standing problem. Assuming hydrostatic balance where g denotes the graviational acceleration, the equations reduce to the well-known primitive equations and only recently the global well-posedness of strong solutions could be proven by Cao and Titi [8] for the incompressible ocean dynamics. The density of air ρ in the atmosphere in comparison to the ocean however varies strongly with height, and the incompressibility assumption is only justified when describing shallow phenomena. Thus for the atmosphere in general the full compressible governing equations need to be considered. However, under the assumption of hydrostatic balance (1), which in particular guarantees the pressure to decrease monotonically with height, the pressure p can be used as the vertical coordinate. Switching to the pressure coordinates (x, y, p) has the main advantage that the continuity equation takes the form of the incompressibility condition see, e.g., Lions et al. [21] and Petcu et al. [22]. Thus the ideas of Cao and Titi [8] can be taken over for the atmospheric primitive equations in pressure coordinates, as we shall also see below and we therefore work in the following with the governing equations in the pressure coordinates and use hereafter the notation We note that the vertical velocity ω in pressure coordinates takes the opposite sign as the vertical velocity in the cartesian coordinates, i.e., ω < 0 for upward motion and ω > 0 for downward motion. Moreover the derivatives as well as the velocity components have different units in vertical and horizontal directions. Nevertheless the total derivative in pressure coordinates reads For the eddy viscosity closure of turbulence and molecular transport we use whereT =T (p) corresponds to some background distribution of the temperature, being uniformly bounded from above and below, and R is the individual gas constant. The operator D * thereby provides a close approximation to the full Laplacian in cartesian coordinates, see also [21,22]. The governing equations in pressure coordinates (x, y, p) with corresponding velocities (u, ω) become where Φ denotes the geopotential ∂ z Φ = g, and the second equation combines the ideal gas law (see (8) below) and the hydrostatic balance equation (1). Here v := (u, ω), k = (0, 0, 1), and (k × v) h are the first two components of k × v. Now the density ρ does not appear in the system anymore, which is connected to the other thermodynamic quanitities via the ideal gas law A typical measure for quantification of moisture are mixing ratios, which compare the density of the moisture component to the density of dry air (denoted as ρ d ). We assume to be in the warm cloud regime, where no ice and snow phases occur and water is therefore present in the form of water vapor (with density ρ v ) and in the liquid state as cloud water and rain water (with corresponding densities ρ c and ρ r ), such that we have the mixing ratios For these mixing ratios we then have the moisture balances where the total derivative is given according to (3) and the diffusion terms are as in (4). The source terms S ev , S cd , S ac , S cr are, respectively, the rates of evaporation of rain water, the condensation of water vapor to cloud water and the inverse evaporation process, the auto-conversion of cloud water into rainwater by accumulation of microscopic droplets, and the collection of cloud water by falling rain. Moreover V denotes the terminal velocity of falling rain and is assumed to be constant. The thermodynamic equation accounts for the diabatic source and sink terms, such as latent heating, radiation effects, but we will in the following only focus on the effect of latent heat in association with phase changes (see, e.g., [10,11,16,19]). The temperature equation in pressure coordinates then reads, see, e.g., [13,15], where the heat capacity c p and the latent heat L are assumed to be constant.
To describe the state of the atmosphere a common thermodynamic quantity used instead of the temperature is the potential temperature The potential temperature has the main advantage that the left-hand side of (13) simply reduces to T θ d dt θ. This property was essential in the preceding works [12] and [16] to derive a priori nonnegativity and boundedness of the temperature and moisture components.
Remark 1. In the present model the difference of the gas constants for dry air and water vapor as well as the dependence of the internal energy on the moisture components via different heat capacities is neglected. These additional terms that would arise in a more precise thermodynamical setting are small in principle and therefore often not taken into account. It has been, however, revealed in [15] that, e.g., in the presence of deep convective clouds the refined thermodynamical setting can be essential and the according moisture model leading in particular to a much stronger coupling of the thermodynamic equation to the moisture components will be investigated in a forthcoming paper.

1.2.
Explicit expressions for the source terms. The threshold for phase changes is saturation, which is defined via the saturation mixing ratio q vs . Saturation thereby is reached when q v = q vs , whereas the air is undersaturated if q v < q vs and oversaturated if q v > q vs , respectively. For the given function q vs we pose the natural assumption to depend continuously on p and T and to vanish below and above some critical temperatures (given in Kelvin), i.e.
for some 0 ≤ T A ≤ T B , which is helpful for proving nonnegativity of the solution. Moreover we assume q vs to be nonnegative, uniformly bounded and to be Lipschitz continuous with respect to T , i.e., we assume for a positive constant C. This constant actually depends on the pressure as it grows approximately inversely proportional to p. We are however only interested in the lower part of the atmosphere, where weather related phenomena are taking place. There the pressure is uniformly bounded from below (e.g. by 100 hPa), such that the constant C in (16) can be assumed to be positive and uniformly bounded. For more details we refer also to [16].
For the source terms of the mixing ratios we take over the setting of Klein and Majda [19] corresponding to a basic form of the bulk microphysics closure in the spirit of Kessler [20] and Grabowski and Smolarkiewicz [14], which has also been used in the preceding work [16]: where C ev , C cr , C ac , C cd , C cn are dimensionless rate constants. Moreover, (g) + = max{0, g} and q * ac denotes the threshold for cloud water mixing ratio beyond which autoconversion of cloud water into precipitation becomes active.
Exponents β ∈ (0, 1) cause difficulties in the analysis, in particular for the uniqueness of the solutions. This problem however was overcome in [16] by introducing new unknowns, which allow for certain cancellation properties of the source terms and reveal advantageous monotonicity properties.
The rest of this paper is arranged as follows. In section 2, we formulate the full problem with boundary conditions and state the main results on the global existence and uniqueness of solutions. In section 3, we prove the existence and uniqueness, and the uniform a priori estimates of solutions to an approximate system, which is nothing but the original one by replacing S ev (which may be only Hölder continuous in q r ) by S ev,ε which is Lipschitz in q r . In section 4, based on the results in section 3, we give the proof of the global existence result, and the uniqueness is also established by using the idea in our previous work [16].
Throughout this paper, unless explicitly specified, we use C to denote a generic positive constant depending only on the given functions in the boundary conditions, the initial data, and the physical parameters appearing in the original system (but not on the parameter ε arising in the approximate system introduced in the next section).

Formulation of the problem and main results
Recall the momentum equation where v = (u, ω), k = (0, 0, 1), and (k × v) h are the first two components of k × v.
According to the incompressibility condition (7) and the boundary condition (24), we have a diagnostic equation for ω x, y, s)ds .
By the hydrostatic balance (6), we have Similar to [10,16] we consider M to be a cylindrical domain of the form where M ′ is a smooth bounded domain in R 2 , and 0 < p 1 < p 0 . The boundary is given by The boundary conditions are: where α 0j , α ℓj , α 0T , α ℓT are given nonnegative constants, and T b0 , T bℓ , q b0j , q bℓj , which can depend on time, are given nonnegative and sufficiently smooth functions. Here n denotes the outward normal direction on ∂M ′ . Note that the boundary conditions (24)-(26) reduce to those in [11] if α ℓT , α ℓv , α ℓc , and α ℓr are set tozero. The initial condition is Throughout this paper, we use the abbreviation According to the weight in the vertical diffusion terms, we introduce the weighted norm which, since the weight gp R dT is uniformly bounded from above and below by positive constants, is equivalent to the L 2 -norm. Moreover, we shall often use for convenience the notation We state our main result on the global existence and uniqueness of solutions to the fully coupled system in the following: p 0 ∇ h · u 0 dp = 0 on M ′ . Then, system (10)-(13), (21)-(23), subject to (24)-(27), has a unique global in time solution (u, T, q v , q c , q r ), satisfying for any T ∈ (0, ∞).
3. An approximated system: existence and a priori estimates This section is devoted to proving the global existence, uniqueness, and uniform a priori estimates of the following ε-approximate system to (10)- (12), (13), and (21)-(23), where ε ∈ (0, 1) is fixed and Unlike the original S ev , the corresponding approximation S + ev,ε is Lipschitz with respect to q r , and it approximates S ev as ε tends to zero.
Since all the nonlinear terms S + ev,ε , S + cr , S + ac , and S + cd are Lipschitz with respect to q v , q c , q r , and T , the local, in time, existence and uniqueness of strong solutions to the initial boundary value problem of the ε-approximate system (28)-(34) follows the standard contraction mapping fixed point principle and we obtain the following proposition on the local, in time, existence and uniqueness result.
Then, there is a positive time T 0 depending only on the upper bound of (u 0 , T 0 , q v0 , q c0 , q r0 ) H 1 (M) , which is independent of ε, such that system (28)-(34), subject to (24)-(27), has a unique strong solution u, T, q v , q c , q r on M × (0, T 0 ), satisfying In the following we let (u, T, q v , q c , q r ) be the solution obtained in Proposition 1, and extend it to the maximal interval of existence (0, T max ), where T max is characterized as The following assumption will be made in the subsequent propositions throughout this section.
Assumption 1. Let all the assumptions in Proposition 1 hold, and let the unique solution (u, T, q v , q c , q r ) obtained in Proposition 1 be extended in the above way to the maximal time of existence T max .
The aim of the rest of this section is to show that T max = ∞, and to establish the a priori estimates that are independent of ε ∈ (0, 1).
As in [16] we have the nonnegativity and uniform boundedness of (T, q v , q c , q r ) stated in the next proposition.
with q * vs = max q vs . Moreover q * c , q * r , T * are continuous in T and depend on the following quantities: Proof. We only state the key steps of the proof and refer to the proof of Proposition 3.2 in [16] for more details. Due to the vanishing of the antidissipative term in the potential temperature θ, the latter is used instead of the temperature for the derivation of the maximum principle. As a first step then the nonnegativity of q c , q v , q r and θ are derived in exactly this order by employing the Stampacchia method. Therefore the according equations are multiplied with the corresponding negative parts of the solution components and after integration and integration by parts the Gronwall inequality is applied to show that if the solution was nonnegative initially it has to remain nonnegative for all times. The same method can then be used to show uniform boundedness of q v by deriving an The boundedness of q c , q r and θ again in this order follow by employing iterative estimations for cutoff functions in L m and then passing to the limit as m → ∞ (see Proposition 3.2 in [16]).
Due to the nonnegativity of q v , q c , q r , and T we obtain also for the source terms S + cr = S cr , S + cd = S cd and S + ev,ε = S ev,ε = C ev T q r (q r + ε) β−1 (q vs − q v ) + , such that in the following we can drop the positive signs in the according notations. Proposition 2 in particular also implies the boundedness of the source terms. To see the uniform boundedness of S ev,ε we notice that since β ∈ (0, 1), we have q r (q r + ε) β−1 ≤ q r q β−1 r = q β r , and recall that q vs (p, T ) = 0, for T ≥ T B (see (15)). Therefore, one has where the constant C, as the upper bounds in Proposition 2, depends on T , the initial data and the given functions in the boundary conditions (24)-(26).
Proposition 3. Let Assumption 1 hold, then there exists a function K 0 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that the estimate holds for any T ∈ [0, T max ).
Proof. The conclusion follows from multiplying (32), (33), and (34), respectively, with q v , q c , and q r , performing integration, integrating by parts and using the boundary conditions as well as Young's inequality. More precisely by the incompressibility condition (30) and the boundary conditions (24)-(26), the integrals involving the convection terms Moreover for the diffusion terms we obtain from integration by parts the following coercive estimate Finally all the integrals involving the source terms S ev,ε , S cr , S ac , S cd are uniformly bounded, due to the a priori bounds obtained in Proposition 2 concluding the proof.
Furthermore we have the basic energy inequality contained in the following proposition.
Proposition 4 (Basic energy estimate). Let Assumption 1 hold, then there exists a function K 1 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that for any T ∈ [0, T max ).
Proof. Multiplying equation (28) by u, integrating the resultant over M, then it follows from integration by parts that here we have used the facts that (k×v) h ·u = u ⊥ ·u = 0, and that the integral involving the convection terms vanishes after integration by parts, due to the incompressibility condition (30) and the boundary conditions (24)-(26). For the integral of the dissipation term, it follows from integration by parts and the boundary conditions (24)-(26) for u that For the term − M ∇ h Φ · udM, it follows from integrating by parts, using the boundary conditions (24)-(26) for u, and equations (29)-(30), that Substituting (43) and (44) into (42) yields Multiplying equation (31) by c p , integrating the resultant over M, and recalling the nonnegativity of T , it follows from integration by parts that By the boundary conditions (24)-(26) for T , it follows from integration by parts that where we note that C depends on the boundary data. Recalling (41) we have moreover Thanks to (47) and (48), it follows from (46) that with the constant C depending again on the boundary data. Adding (45) and (49) yields d dt from which by integration the conclusion of the proposition follows.
Thanks to the basic energy estimate, we can now derive the L ∞ (L 2 ) ∩ L 2 (H 1 ) estimate on T contained in the following proposition.
Proposition 5. Let Assumption 1 hold, then there exists a function K 2 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that for any T ∈ [0, T max ).

Proof. Performing integration by parts and using the boundary conditions for T we obtain
By the trace inequality, the boundary integrals in the above inequality can be bounded as Multiplying (31) with T and using (50) we get from which, recalling (41), we obtain By Lemma 2 from the Appendix and Young's inequality, we deduce which, substituted into (51), leads to from which, by the Gronwall inequality and Proposition 4, the conclusion follows.
We have the following L ∞ (L 2 ) estimate on the vertical derivative of the velocity, which was first established by Cao and Titi in [8] for the primitive equations without coupling to the moisture equations. Proposition 6. Let Assumption 1 hold, then there exists a function K 3 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that for any T ∈ [0, T max ).
Proof. The proof is adapted from that in [8]. Decompose the velocity into the barotropic and baroclinic modesū and u as follows p 0 p 1 f dp. For the L 6 -norm of u, following the arguments in [8], with tiny modifications on dealing with the integral involving Φ, we have (see the inequality above (66) in page 257 of [8]) Noticing f + f ≤ C f , ∇f + ∇f ≤ C ∇f , it follows from Proposition 4, Proposition 5 and the Gronwall inequality that sup 0≤t≤T u 6 for a continuous function K ′ 3 on [0, ∞), which depends on the initial and boundary data. For the barotropic modeū it holds that (see the first inequality in page 259 of [8]) from which, by Proposition 4 and (52), it follows from the Gronwall inequality that for a continuous function K ′′ 3 on [0, ∞). Note that by the Sobolev embedding inequality and Proposition 4 one can easily obtain from (52) and (53) that for a continuous function K ′′′ 3 on [0, ∞) depending on the initial and boundary data. For the vertical gradient we have the following bound (see the inequality above (75) in page 260 of [8]) from which, by (52)-(53), Proposition 5 and the Gronwall inequality it follows that for a continuous function K ′′′′ 3 on [0, ∞) depending on the initial and boundary data, which completes the proof.
We are now ready to establish the estimates for the vertical derivative of the moisture components.
Proposition 7. Let Assumption 1 hold, then there exists a function K 4 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that for any T ∈ [0, T max ).
Proof. We only prove the estimate for q r , since the derivation for the other moisture quantities follows the same steps. Multiplying equation (34) by −∂ 2 p q r and integrating over M yield Recalling (41), it follows from the Hölder and Young inequalities that For the fist term on the left-hand side of (54), using the boundary conditions (24) and (25), it follows from integration by parts that where in the last step we have used Lemma 1 (see Appendix), Proposition 2 and the Hölder inequality to estimate We next bound the integrals involving the diffusion terms. It should be noticed that the following derivations are formal, since the regularity of q r does not guarantee the validity of the integrals on the boundary Γ 0 = M ′ × {p 0 } directly. However, due to the C 2 (M × [0, T ]) functions, enjoying the boundary conditions (24)-(26) for q r , being dense in the space {f |f ∈ L 2 (0, T ; H 2 (M)), ∂ t f ∈ L 2 (0, T ; L 2 (M))}, which q r belongs to, one can rigorously justify the desired estimates (62) below for q r in the standard way ( i.e., choosing a sequence {q rn } ∞ n=1 in C 2 (M × [0, T ]) satisfying the corresponding boundary conditions (24)-(26), proving that (62) holds for each q rn , and passing the limit n → ∞).
Using the boundary conditions (25), it follows from integration by parts that Integrating by parts and using the boundary conditions (24) and (26), we deduce where we have used Using the boundary condition (26), one has Substituting (58) and (61) into (57) yields For the term involving the convection, it follows from integration by parts and the boundary condition (25) that By Proposition 2, Lemma 1 (see Appendix), and Young's inequality, we obtain α 0r where in the last step we have used Proposition 4 and Proposition 6. Note that the boundary condition (26) for u implies ∂ p u · n = 0 on Γ ℓ . Thanks to this, by Proposition 2 and Proposition 6, it follows from integration by parts and the Young inequality that for any p ∈ [p 1 , p 0 ]. Based upon this fact, it follows from Lemma 2 (see Appendix) that Substituting (64), (65), and (67) into (63), one obtains the estimate Substituting (55), (56), (62), and (68) into (54) leads to d dt from which, noticing that the conclusion follows, by Proposition 3, Proposition 6, and the Gronwall inequality.
Based on the estimates on the vertical derivatives of the velocity and moisture components, we can now bound the corresponding horizontal derivatives. Proposition 8. Let Assumption 1 hold, then there exists a function K 5 (t), which depends on the initial and boundary data, and is continuous for all t ≥ 0, such that for any T ∈ [0, T max ).
Proof. Following the arguments in Section 3.3.3 of [8], we have for the estimate of the horizontal gradient, from which, by Proposition 6 and Young's inequality, one obtains leading to Next, we estimate the horizontal gradient of q r . Multiplying equation (34) by −∆ h q r and integrating over M yield Recalling (41), it follows from the Young inequality, Proposition 2 and Proposition 7 that Integrating by parts and using the boundary condition (26), we obtain where in the last step we have used which is guaranteed by the trace inequality and Proposition 2.
Recalling the expression of D qr q r , we have where, we have used Proposition 2, Proposition 6, and Proposition 7. By the elliptic estimates and Proposition 2, we have for i ∈ {v, c, r}. Then it follows from (79) and the Young inequality that Thanks to (72), (73), (78), and (81), it follows from (71) that d dt The same estimate as above also holds for q v , q c , such that d dt from which, noticing that the conclusion follows by the Gronwall inequality and Propositions 4-7.
Finally, we have the control of the gradient of the temperature.
Proof. We need to prove T max = ∞. Assume, by contradiction, that T max < ∞. By Propositions 2-9, we have the estimate sup 0≤t≤T (u, T, q v , q c , q r ) H 1 (M) ≤ C 0 , for any T ∈ (0, T max ), and C 0 is a positive constant depending on the initial and boundary data, but which is independent of T ∈ (0, T max ). This contradicts (35) and, thus, T max = ∞. The a priori estimates except those involving the time derivatives follow directly from Propositions 2-9; while the desired estimates for the time derivative follow from those in Propositions 2-9, using equations (28), (31)-(34). The proof is lengthy but standard, and therefore omitted.