Nonuniqueness of generalised weak solutions to the primitive and Prandtl equations

We develop a convex integration scheme for constructing nonunique weak solutions to the hydrostatic Euler equations (also known as the inviscid primitive equations of oceanic and atmospheric dynamics) in both two and three dimensions. We also develop such a scheme for the construction of nonunique weak solutions to the three-dimensional viscous primitive equations, as well as the two-dimensional Prandtl equations. While in [D.W. Boutros, S. Markfelder and E.S. Titi, Calc. Var. Partial Differential Equations, 62 (2023), 219] the classical notion of weak solution to the hydrostatic Euler equations was generalised, we introduce here a further generalisation. For such generalised weak solutions we show the existence and nonuniqueness for a large class of initial data. Moreover, we construct infinitely many examples of generalised weak solutions which do not conserve energy. The barotropic and baroclinic modes of solutions to the hydrostatic Euler equations (which are the average and the fluctuation of the horizontal velocity in the $z$-coordinate, respectively) that are constructed have different regularities.


Problems considered in this paper and context
In this work, we consider the following general equation (with (d + 1)-dimensional spatial domain for d = 1, 2) where the horizontal velocity field u : T d+1 × (0, T ) → R d , the vertical velocity field w : T d+1 × (0, T ) → R, and the pressure p : T d+1 → R are unknown, and the horizontal and vertical viscosity parameters ν * h , ν * v ≥ 0 are given constants. The d-dimensional horizontal gradient is denoted by ∇ h and the d-dimensional horizontal Laplacian by ∆ h . Since the pressure p is only determined up to an additive constant, we may require that p is mean-free.
In this paper we are interested in the following cases: • Taking d = 2 and ν * h = ν * v = 0 gives the three-dimensional hydrostatic Euler equations of an incompressible fluid (also known as the inviscid primitive equations of oceanic and atmospheric dynamics). In this paper, the terms inviscid primitive equations and hydrostatic Euler equations will be used interchangeably.
In this paper, we will develop a convex integration scheme for system (1.1)-(1.3) for the cases mentioned above. In particular, we will work with a generalised notion of weak solution. While classical weak solutions have sufficient Lebesgue integrability for the nonlinearity to make sense as an L 1 (T d+1 × (0, T )) function, another notion of weak solution was introduced in [7] where the nonlinearity is interpreted as a paraproduct. The generalised weak solutions introduced in this paper treat the nonlinearity in an even more general way, see section 1.3.3 below.
In all the cases of (1.1)-(1.3) that we are interested in, we will show the existence of such generalised weak solutions (for a dense set of initial data in the relevant spaces). In addition, we will show that such weak solutions are nonunique.
If ν * h = ν * v = 0, we recall that classical spatially analytic solutions of (1.1)-(1.3) (see [40,50,51]) conserve the energy, i.e. the spatial L 2 (T d+1 ) norm of u. In [7] an analogue of Onsager's conjecture was studied for the three-dimensional hydrostatic Euler equations and it was found that there exist several sufficient regularity criteria for weak solutions which guarantee the conservation of energy. In particular, there exist several notions of weak solutions for these equations, each of which have their own version of the analogue of the Onsager conjecture.
In this work, we will construct generalised weak solutions to these equations, which do not conserve energy and do not satisfy the regularity criteria mentioned above. In other words, in this paper we prove a first result towards the aim of resolving the dissipation part of the analogue of the Onsager conjecture for the inviscid primitive equations (hydrostatic Euler), while the conservation part of the analogue of the Onsager conjecture has been studied in [7], as was mentioned before.

Literature overview
In this section we will provide an overview of some of the literature that is related to this work. As both the primitive and Prandtl equations as well as the Onsager conjecture have been the subject matter of many works in recent years, this overview is by no means comprehensive and is by necessity incomplete in reviewing all the relevant work.
Onsager's conjecture was originally posed in [75] for the incompressible Euler equations. The conjecture states that if a weak solution lies in L 3 ((0, T ); C 0,α (T 3 )) for α > 1 3 it must conserve energy. If α < 1 3 energy might not be conserved. In [37] a proof of a slightly weaker result than the first half of the conjecture was given. A full proof of the first half was then given in [30]. In [35] a different proof was presented, well-posedness for the case without rotation and Dirichlet boundary conditions. Finally, [59] considered the case of impermeable and stress-free boundary conditions. The linear and nonlinear ill-posedness of the inviscid primitive equations in all Sobolev spaces was proven in [44,80]. The ill-posedness results in Sobolev spaces suggest that the natural space for showing local well-posedness of the inviscid primitive equations is the space of analytic functions, which was proved in [40,50,51]. In [40] the role of fast rotation in prolonging the life-space of solutions was investigated.
In [16] it was shown that smooth solutions of the inviscid primitive equations can form a singularity in finite time, see also [90]. In [29] the existence and nonuniqueness of weak solutions with L ∞ data was proven. In [7] several sufficient criteria for energy conservation were proven. In the inviscid setting there have also been works studying the case of initial data with a monotonicity assumption, see [9,49,66].
The Prandtl equations for the boundary layer were derived by Prandtl in [78]. In [73,74] the local well-posedness of the equations was shown under a monotonicity assumption. In [82] the local well-posedness for analytic data was proven, while in [36] the blow-up of solutions for certain classes of C ∞ data was proven. Further local well-posedness results were proved in [46,49,52,76,77]. In [42] it was shown that the equations are nonlinearly unstable.
The linear ill-posedness of the Prandtl equations in all Sobolev spaces was shown in [38] (for further work see [63] and references therein). In the three-dimensional case a convex integration scheme was developed in [65]. The analytic local well-posedness has been improved to Gevrey function spaces, see [58] and references therein.

Baroclinic and barotropic modes
Now we introduce the notion of barotropic and baroclinic modes, which is an important decomposition of the solutions which has been explored extensively in the investigation of the primitive equations. In the construction of the convex integration scheme for the primitive equations we will not use this decomposition explicitly. However, it is an important idea underlying the scheme.
We will illustrate this concept for the equations in the inviscid case, the viscous case is similar and can be found in [23]. The 3D inviscid primitive equations are given by (1.6) where u : T 3 × (0, T ) → R 2 is the horizontal velocity field, w : T 3 × (0, T ) → R the vertical velocity field and p : T 3 × (0, T ) → R the pressure.
The barotropic mode u of a velocity field u is defined as follows u(x 1 , x 2 , t) :=ˆT u(x 1 , x 2 , z, t) dz. (1.7) The baroclinic mode u is defined as the fluctuation u := u − u. (1.8) The primitive equations (1.4)-(1.6) can then be written formally as a coupled system of evolution equations for the barotropic and baroclinic modes u and u, which are (1.10) Moreover, we have the following incompressibility conditions (1.11) which formally follow from equation (1.6) and the periodicity of the functions.
In the convex integration scheme, we will add separate barotropic and baroclinic perturbations. This leads to different regularities of the barotropic and baroclinic modes of the solution and allows us to control different parts of the error.
The following estimates on the baroclinic and barotropic modes are standard

Notation
Throughout the paper we will use the following notation.
• The components of the spatial variable are given by x = (x 1 , z) if d = 1, and x = (x 1 , x 2 , z) if d = 2. For d = 1, x 1 represents the horizontal direction, for d = 2 the horizontal position is given by (x 1 , x 2 ). In both cases, z is the vertical direction.
• The horizontal velocity field is called u, the vertical velocity is denoted by w and the full velocity by u = (u, w). They are d-, 1-and (d + 1)-dimensional, respectively.
• We use the symbol ∇ h for the horizontal gradient (which equals ∂ x 1 if d = 1), and ∇ for the full ((d + 1)-dimensional) gradient.
1 p + 1 p ′ = 1. • Let 1 < p ≤ ∞. In section 1, p− denotes any parameter 1 ≤ p− < p. In the other sections we have to be a bit more precise. In particular there is a need to quantify the '−' in p−. More precisely there will be a δ > 0 and we set p− := 1 1 p +δ . Here we tacitly assume that δ is sufficiently small, such that p− ≥ 1.
• For 1 ≤ p, q ≤ ∞ and s ∈ R the Besov space B s p,q (T 3 ) is defined in appendix A.1. Let us emphasise here that B s 2,2 (T 3 ) = H s (T 3 ), see Remark A.2.
• Throughout this paper, we will omit the domain of a space-time norm if it is T d+1 ×[0, T ], e.g. we write · L p (H s ) = · L p ((0,T );H s (T d+1 )) .
• In view of section 1.3.1 we define the barotropic and baroclinic part of any quantity a = a(x) by a =ˆT a(x) dz, a = a − a.

Generalised weak solutions
In [7] two new types of weak solutions to the hydrostatic Euler equations (1.4)-(1.6) were introduced. In the present paper we will consider a slightly different notion of weak solution, which we will refer to as a generalised weak solution. This notion of solution is inspired by the notion of a type III weak solution, as introduced in [7].
Before we state the results for the different cases of the system (1.1)-(1.3), we will be more specific regarding the notion of weak solution used in this paper. The weak solutions of (1.1)-(1.3) we consider are defined as follows: We assume that u ∈ L 2 (T 3 × (0, T )), w ∈ D ′ (T 3 × (0, T )) and uw ∈ L 1 ((0, T ); B −s 1,∞ (T 3 )) for some suitably large s ∈ R. System (1.1)-(1.3) must then be satisfied in the sense of distributions, where the vertical advection termˆT is interpreted as a duality bracket between the term uw and the test function φ ∈ D(T 3 × (0, T )).
If u and w happen to have sufficient regularity, for example when u ∈ L 2 ((0, T ); H s+δ (T 3 )) and w ∈ L 2 ((0, T ); H −s (T 3 )) (for some small δ > 0), then by applying the paradifferential calculus (see appendix A) we know that uw ∈ L 1 ((0, T ); B −s 1,∞ (T 3 )) . This is a stronger notion of solution compared to the notion of a generalised weak solution that we introduced above, as u is required to have (positive) Sobolev regularity and w has to possess some regularity (i.e., it is more than just a distribution).
In the next few subsections, we will give precise definitions of the notion of weak solution we will use, and we will state the theorems we will prove for the different cases of the system (1.1)-(1.3). But generally speaking, we will split the nonlinearity uw into the barotropicvertical and baroclinic-vertical interactions, i.e., the terms uw and uw.
The baroclinic mode u of the constructed solutions will have sufficient regularity such that uw can be interpreted as a paraproduct. The terms u and w do not have sufficient regularity to apply the paradifferential calculus. However, as part of the convex integration scheme we will obtain separate estimates on uw in order to show that it lies in L 1 ((0, T ); B −s 1,∞ (T 3 )) for some suitable s. Therefore the weak solutions we obtain are partly 'generalised' (as for the baroclinic-vertical part of the nonlinearity) and partly 'paradifferential' (for the barotropicvertical part of the nonlinearity).

Results for the 3D inviscid primitive equations
We first introduce the notion of weak solution for the 3D inviscid primitive equations (1.4)-(1.6).
where s > 0 is referred to as the regularity parameter) and the equations are satisfied in the sense of distributions, i.e.
Remark 1.2. We emphasise that this definition of weak solutions to (1.4)-(1.6) is more general than the notion of weak solution introduced in [7]. While in [7] the velocity field of a weak solution has sufficient regularity to automatically guarantee that uw ∈ L 1 ((0, T ); B −s 1,∞ (T 3 )) (by using the paradifferential calculus), in Definition 1.1 we do not have sufficient (separate) regularity requirements on u and w such that the product uw is well-defined. Hence uw ∈ L 1 ((0, T ); B −s 1,∞ (T 3 )) is a separate independent requirement of Definition 1.1.
In this paper we will prove the following result. . Moreover, let 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 1 0 < s 1 , s 3 be parameters satisfying Then there exists a weak solution (u, w, p) in the sense of Definition 1.1 with regularity parameter s = 1 and with the following properties: 1. The solution satisfies that where u and u denote the barotropic and baroclinic modes of u respectively. Remark 1.4. Alternatively one can construct a weak solution with the properties stated in Theorem 1.3 where the only difference is that the endpoint time integrability is attained for u rather than w. In other words see Remarks 2.6 and 5.6 below. To this end however, we have to require that q 3 < q 1 (strictly) in (1.15).
Remark 1.5. By proceeding as in section 7 below, we can achieve in addition that u, u ∈ L 1 ((0, T ); W 1,1 (T 3 )). To this end, however, we have to require the constraints (1.20) rather than (1.15), see also Theorem 1.10 and Remark 1.13, below.
Remark 1.6. Again we would like to remark that the solutions constructed in Theorem 1.3 are partially 'generalised' (see section 1.3.3) and partially 'paradifferential' as in [7]. In particular, they have been inspired by the type III weak solutions that were introduced in [7]. More precisely, from the regularities of u and w stated in Theorem 1.3 it follows that uw ∈ L 1 ((0, T ); B −1 1,∞ (T 3 )) (see the proof of Theorem 1.3 in sections 4-6 for details). The term uw is estimated directly in L 1 ((0, T ); B −1 1,∞ (T 3 )) as part of the convex integration scheme, as one cannot obtain the regularity of the product uw simply from the regularities of u and w (as they are insufficient to apply the paradifferential calculus directly). The specific form of the perturbations allows for a direct estimate, as was done for example in [28]. Therefore the interpretation of the term uw can be seen as 'paradifferential', while the interpretation of the term uw is in the sense of a 'generalised weak solution' (as in Definition 1.1).
Remark 1.7. In addition, we would like to emphasise that in the presence of physical boundaries the primitive equations are often studied with no-normal flow boundary conditions on the top and bottom of the channel, i.e. w| z=0,1 = 0. However, in the convex integration scheme developed in this paper we will work on the three-dimensional torus rather than the channel. Note that solutions in the torus can be understood as solutions in the channel with an in-flow out-flow boundary condition, i.e.
for a flow w B . In our case w B will be constructed as part of the convex integration scheme. In other words we will not solve the boundary value problem for given w B and in particular, not for the case of the impermeability boundary condition w B = 0. We also remark that the constructed flow w B belongs to the space L q ′ 2 ((0, T ); L 2 (T 2 )) ∩ L q ′ 3 ((0, T ); H −s 3 (T 2 )), where the parameters q ′ 2 , q ′ 3 and s 3 are the same as in Theorem 1.3. Theorem 1.3 allows to show the nonuniqueness and existence of solutions which do not conserve energy: Corollary 1.8. For any analytic initial data there exist infinitely many global-in-time weak solutions (u, w, p) of the hydrostatic Euler equations (1.4)-(1.6) (in the sense of Definition 1.1 which satisfy the regularity properties of Theorem 1.3) and they do not conserve energy.
Proof. We take the smooth local-in-time solution for the given choice of analytic data (whose existence can be proven using the methods from [40,51]) as the first solution (u 1 , w 1 , p 1 ) on [0, T /2], and the zero solution on [T /2, T ] as (u 2 , w 2 , p 2 ). Then Theorem 1.3 yields a weak solution, which we may extend by zero for t > T . If the initial data are non-zero, we can conclude that the energy is not conserved as it is positive on [0, T /4) and zero on (3T /4, ∞).
Another global-in-time weak solution can be constructed similarly with replacing T by T /4. This solution has positive energy on [0, T /16) while the energy is zero on (3T /16, ∞). Consequently the two solutions cannot coincide. Repeating this argument leads to infinitely many global-in-time weak solutions with the same initial data, which are smooth and unique for a small initial interval of time, but which do not conserve energy.
For zero initial data, we observe that Theorem 1.3 allows one to 'connect' any analytic initial data with any analytic data in finite time. Hence we may connect the zero initial data to arbitrary analytic data with positive energy at t = T . On the time interval [ T , ∞) we then proceed as above.

Results for the 3D viscous primitive equations
We now consider the viscous primitive equations, which are given by where ν * h and ν * v are the horizontal and vertical viscosities. As before, u : T 3 × (0, T ) → R 2 is the horizontal velocity field, w : T 3 × (0, T ) → R the vertical velocity and p : T 3 × (0, T ) → R the pressure. We have the following notion of weak solution for these equations. Definition 1.9. A triple (u, w, p) is called a weak solution of the viscous primitive equations where s > 0 is referred to as the regularity parameter) and the equations are satisfied in the sense of distributions, i.e.
In this paper we will prove the following result. . Moreover, let 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 0 < s 1 , s 3 be parameters satisfying the following relations (1.20) Then there exists a weak solution (u, w, p) in the sense of Definition 1.9 with regularity parameter s = 1 and with the following properties: 1. The solution satisfies that 2. We have that Remark 1.11. Similar to Theorem 1.3 one can even obtain endpoint time integrability for w.
With the modification described in Remarks 1.4, 2.6 and 5.6 one can alternatively establish endpoint time integrability for u. Remark 1.12. The reader should notice that there exist parameters 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 0 < s 1 , s 3 satisfying (1.20). Indeed for every q 3 < 3/2, we have Hence there exists q 1 with 1 Thus q 3 < q 1 and for q 2 > 2 sufficiently large, the estimate holds since the right-hand side converges to 1 2q 3 for q 2 → ∞. This allows to choose s 1 and s 3 such that 2 so all constraints in (1.20) are satisfied.
Remark 1.13. We would like to emphasise that Theorem 1.10 holds for any choice of viscosities ν * h , ν * v ∈ R, in particular even for the inviscid case ν * h = ν * v = 0. Remark 1.14. Finally, we emphasise that the solutions constructed in Theorem 1.10 are not of Leray-Hopf type, as they do not have a finite rate of mean energy dissipation (i.e. the horizontal velocity field does not belong to the space L 2 ((0, T ); H 1 (T 3 ))).
We now obtain the global existence of weak solutions as a corollary. Proof. The proof works exactly as the proof of Corollary 1.8 where the corresponding local (even global) well-posedness result can be achieved by using the methods from [23].
Remark 1.16. The proof of nonuniqueness of global weak solutions works equally well in the three cases of full, horizontal or vertical viscosity, which were studied in the works [19][20][21]. Moreover, in the case of full viscosity the result can also be adapted to classes of initial data belonging to different function spaces, by relying on the well-posedness results from [41].

Results for the 2D hydrostatic Euler equations
It is also possible to develop a convex integration scheme for the two-dimensional hydrostatic Euler equations. They are given by where u : T 2 × (0, T ) → R is the horizontal velocity, w : T 2 × (0, T ) → R is the vertical velocity and p : T 2 × (0, T ) → R is the pressure. We first state the definition of weak solution to these equations.
Definition 1.17. A triple (u, w, p) is called a weak solution of the two-dimensional hydrostatic Euler equations (1.21)-(1.23) if u ∈ L 2 (T 2 × (0, T )), w ∈ D ′ (T 2 × (0, T )) and p ∈ L 1 (T 2 × (0, T )) such that uw ∈ L 1 ((0, T ); B −s 1,∞ (T 2 )) (where s > 0 is the regularity parameter) and the equations are satisfied in the sense of distributions, i.e., In particular, we will prove the following theorem. Moreover, let 1 ≤ q 2 , q 3 ≤ ∞ and 0 < s 3 be parameters satisfying 2 Then there exists a weak solution (u, w, p) in the sense of Definition 1.17 with regularity parameter s = 1 and with the following properties: 1. The solution satisfies that

We have that
Remark 1.19. It might seem slightly odd to label the parameters by q 2 , q 3 and s 3 (rather than q 1 etc.). The reason we chose to do so is because it will allow for easy comparisons with the three-dimensional scheme from Theorem 1.3. We emphasise that there are no equivalent parameters to q 1 and s 1 in the two-dimensional version of the scheme. Remark 1.21. By proceeding as in section 9 we can achieve in addition that u, u ∈ L 1 ((0, T ); W 1,1 (T 3 )).
In contrast to the three-dimensional case (cf. Remark 1.5) in two dimensions there is no need to require stronger constraints for the parameters, see also Theorem 1.23 and Remark 1.25.
We observe that it is possible to establish a two-dimensional analogue of Corollary 1.8 using the local well-posedness result from [50,51] for analytic data in the channel. This yields existence of infinitely many global weak solutions for suitable initial data.

Results for the two-dimensional Prandtl equations
Now we turn to studying the two-dimensional Prandtl equations, which are given by where u : T 2 × (0, T ) → R is the horizontal velocity, w : T 2 × (0, T ) → R is the vertical velocity and p : T 2 × (0, T ) → R the pressure.
We observe that these equations differ from the two-dimensional hydrostatic Euler equations (1.21)-(1.23) by the vertical viscosity term ν * v ∂ zz u. We introduce the following notion of weak solution to the Prandtl equations (1.25)-(1.27).
We will prove the following result. . Moreover, let 1 ≤ q 2 , q 3 ≤ ∞ and 0 < s 3 be parameters satisfying Then there exists a weak solution (u, w, p) in the sense of Definition 1.22 with regularity parameter s = 1 and with the following properties: 1. The solution satisfies that Remark 1.25. Similar to Remark 1.13, Theorem 1.23 holds for any ν * v ∈ R, in particular even for the inviscid case ν * v = 0. We note that it is possible to establish an analogue of Corollary 1.8 where one has to use the local well-posedness result from [77, p. 6] (see also [89, p. 7186]) for analytic data in the strip/channel. A straightforward adaption of the proof of Corollary 1.8 yields the existence of infinitely many global weak solutions for suitable initial data.

Further remarks and outline of the paper
Now that we have presented the results for the four cases of system (1.1)-(1.3) that we consider in this paper, we would like to make some further remarks on these results. Some conclusions that can be drawn are: 1. There exist weak solutions of the inviscid primitive equations (1.4)-(1.6) that do not conserve energy. Compared to the solutions constructed in [29], the solutions that we construct in this paper have Sobolev regularity. Moreover, they are related to the notion of type III weak solutions, as introduced in [7], as the barotropic-vertical part of the nonlinearity is interpreted as a paraproduct.
2. In addition, the scheme is able to construct solutions where the baroclinic and barotropic modes have different regularities. This is expected, as the loss of derivative in the advective term only occurs in the baroclinic equation. The barotropic mode must have higher Sobolev regularity than the baroclinic mode in the scheme as otherwise the paraproduct between the vertical velocity and the barotropic mode will not make sense.
3. As far as we can tell, this is the first proof of nonuniqueness of weak solutions for the viscous primitive equations. It shows that although the system is globally well-posed (as shown in [23]), at low regularity the system has nonunique weak solutions. This is true even if one has sufficiently regular Sobolev data for which global well-posedness holds in the class of strong solutions.
4. To the best of our knowledge, this is the first convex integration scheme for the twodimensional Prandtl equations (in the three-dimensional case there is the work [65]), as well as the two-dimensional hydrostatic Euler equations.
There are a few new features of the scheme that we wish to highlight: • We have introduced a splitting of the Reynolds stress tensor into a barotropic and baroclinic part. We add perturbations to separately deal with both these parts of the error. We then ensure that the interactions between the two perturbations are controlled.
• The splitting of the Reynolds stress tensor requires us to construct and use horizontal and vertical inverse divergence operators, as the barotropic part depends only on the horizontal variables, while the baroclinic part is mean-free with respect to the z-variable.
• Having two parts of the perturbation allows us to use different scalings of the temporal intermittency functions for the barotropic and baroclinic parts of the perturbation. This makes it possible to ensure that the different perturbations have different regularities, such that the interactions between the different parts can be controlled. In particular, this is crucial to control the terms u p ⊗ u p and w p u p (the barotropic-baroclinic and vertical-baroclinic parts of the nonlinearity).
Now we present an outline of the paper. In sections 2-6 we will develop the convex integration scheme for the 3D inviscid primitive equations, in order to prove Theorem 1.3. In section 2 we state the core inductive proposition of the convex integration scheme and prove Theorem 1.3 using this proposition. In section 3 we discuss several preliminaries. In particular, we introduce the inverse divergence operators, the spatial building blocks for the convex integration, as well as the temporal intermittency functions. In addition, we will discuss the choice of the frequency parameters.
In section 4 we introduce the perturbation that will be used in each iteration of the convex integration scheme, and compute the new Reynolds stress tensor after adding the perturbation. We will prove the estimates on the perturbation required for Proposition 2.4 in section 5. The estimates on the Reynolds stress tensor will be proven in section 6.
In sections 7-9 we will develop convex integration schemes to study the other cases of equations (1.1)-(1.3) that we are interested in this paper. These schemes differ from the scheme presented in sections 3-6 in some aspects, while other parts are similar. Therefore for the sake of brevity, in sections 7-9 we will focus on the parts that differ from the convex integration scheme for the 3D inviscid primitive equations.
In section 7 we provide an extension of the convex integration scheme to the viscous primitive equations with full viscosity and prove Theorem 1.10. The cases with anisotropic viscosities can be studied in a similar manner. In section 8 we investigate the two-dimensional hydrostatic Euler equations and prove Theorem 1.18. Finally, in section 9 we consider the (two-dimensional) Prandtl equations and provide the proof for Theorem 1.23.
In appendix A we give a short introduction to Littlewood-Paley theory, Besov spaces and paradifferential calculus, in order to make the paper self-contained. In appendix B we state the improved Hölder inequality, which was introduced in [68, Lemma 2.1], and we prove an oscillatory paraproduct estimate based on this inequality. Moreover, we provide another proof of Lemma 5.3 as an alternative to the proof given in section 5.1.2. This Lemma states a new inequality needed to control the interaction between the vertical velocity and the baroclinic mode, which turns out to be a critical part of the scheme.

The inductive proposition
The following underdetermined system of equations is called the hydrostatic Euler-Reynolds system where u, w, p, R h and R v are the unknowns. Here the horizontal Reynolds stress tensor 3 R h : T 2 × [0, T ] → S 2×2 is a function of (x 1 , x 2 , t), while the vertical Reynolds stress tensor , and which is mean-free with respect to z, i.e.´1 0 R v dz = 0. We will only work with smooth solutions to this system.
Moreover by definition R v is mean-free with respect to z and thus The following definition is inspired by [27, Definition 2.1].
Definition 2.2. We say that a smooth solution (u, w, p, R h , R v ) of the hydrostatic Euler- The core of the proof of Theorem 1.3 will revolve around proving the following inductive proposition.
is a smooth solution of the hydrostatic Euler-Reynolds system (2.1)-(2.3) which is well-prepared with associated time interval I and parameter τ > 0. Moreover consider parameters 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 0 < s 1 , s 3 which satisfy the following constraints 4 Finally let δ, ǫ > 0 be arbitrary. Then there exists another smooth solution (u + u p + u p , w + w p , p + P, R h,1 , R v,1 ) of the hydrostatic Euler-Reynolds system (2.1)-(2.3) which is well-prepared with respect to the same time interval I and parameter τ /2, and has the following properties: 3 We denote the set of all symmetric 2 × 2 matrices by S 2×2 . 4 Note that the constraints (2.4) are weaker than (1.15). Indeed from (1.15) we deduce 2. The following estimates are satisfied 56 3. Moreover, we have the following bounds Remark 2.5. Note that when writing u p , u p , we implicitly require u p to be independent of z and u p mean-free with respect to z, see section 1.3.1.
Remark 2.6. Another version of Proposition 2.4 where (2.12) and (2.13) are replaced by is true as well, see Remark 5.6. This way we obtain the endpoint time integrability for u rather than w in Theorem 1.3, see Remark 1.4.
Next we prove Theorem 1.3 using Proposition 2.4.
Proof of Theorem 1.3. We first take χ 1 and χ 2 to be a C ∞ partition of unity of [0, T ] such that χ 1 ≡ 1 on [0, 3T /8] and χ 2 ≡ 1 on [5T /8, T ]. Then we define (u 0 , w 0 , p 0 ) as follows For a suitable choice of χ 1 and χ 2 , (u 0 , w 0 , p 0 ) is no longer a solution of the hydrostatic Euler equations, but with a proper definition 7 of R h,0 , R v,0 it solves the hydrostatic Euler- Taking the sequence ǫ n = 2 −n , δ n = δ with a suitable choice of δ > 0 (see below) and applying Proposition 2.4 inductively, we find a sequence of well-prepared solutions of the hydrostatic Euler-Reynolds system with a sequence of well-preparedness parameters {τ n } (and the same time interval I). Note that τ n → 0.
Estimates (2.5), (2.7)-(2.14) imply that the sequence u 0 + n k=1 u k is a Cauchy sequence in the space ). In particular, by choosing δ appropriately we can identify the limits u, u and w, as n → ∞, lying in the spaces stated in Theorem 1.3. Now we define the pressure p by where u = u + u.
Next, we check that the triple (u + u, w, p) is a weak solution in the sense of Def- Here we have also used that 1 ≥ s 1 > s 3 (which follows from (1.15)) as well as Lemma A.3, and the fact that 1 ≤ 1 (which follows from q 3 ≤ q 1 , see (1.15)) in order to obtain L 1 integrability in time. In addition, we have that ) by estimates (2.6) and (2.15).
Furthermore, we observe that (1.13) immediately follows from the definition of p. Moreover since for any n ∈ N the quintuple (2.19) satisfies (2.3), we find that (u, w, p) complies with (1.14). In order to show (1.12) we define the abbreviations for any n ∈ N and any test function ϕ ∈ D(T 3 × (0, T )). Note that (2.5) and (2.6) imply ), respectively. Hence by taking the limit we deduce from (2.21) that for any test function ϕ ∈ D(T 3 × (0, T )) which is either mean-free with respect to z, or independent of z with ∇ h · ϕ = 0. Here we have used that for any n ∈ N and ϕ mean-free with respect to z, according to (2.2), and, furthermore, that for any n ∈ N and ϕ independent of z with Now we are ready to prove (1.12). We may split the test function φ 1 = φ 1 + φ 1 into the barotropic and baroclinic parts, and use the Helmholtz decomposition to find test functions ϕ, ψ, which are independent of z, and such that φ 1 = ϕ + ∇ h ψ and ∇ h · ϕ = 0. Then by where the latter equality follows from the fact that p is independent of z, (1.14) and the definition of p in equation (2.20).
Finally, we observe that equation (1.16) follows from Proposition 2.4 because the time interval I of well-preparedness stays the same for the sequence (2.19). In particular, all the perturbations have support in the time interval I. Therefore since (u 0 , w 0 , p 0 ) agrees with (u 1 , w 1 , p 1 ) on [0, T /4) and with (u 2 , w 2 , p 2 ) on (3T /4, T ], the constructed solution (u, w, p) will have the same properties, as no perturbations with support in [0, T /4) ∪ (3T /4, T ] have been added.
After recalling some preliminaries in section 3, we will prove Proposition 2.4 in sections 4-6.

Outline
In this paper we are going to use Mikado flows as building blocks. These have been introduced in [31] and are built upon a geometric lemma which goes back to [71], see also [87,Lemma 3.3]. Later on concentrated Mikado flows have been introduced in [68] in order to construct solutions with Sobolev regularity. In this paper the term Mikado flows will always refer to such Mikado flows with concentration.
In the proof of Proposition 2.4 we will handle the two error terms R h and R v separately.
To treat R v we use two-dimensional Mikado flows and Mikado densities in two directions, whereas for R h we use two-dimensional Mikado flows in several directions which are given by the above mentioned geometric lemma. We use the version of the Mikado flows which was introduced in [27]. We recall these flows in section 3.4. We call the perturbation which reduces the error R v vertical and the one which reduces R h horizontal.
The above mentioned concentration is represented by the spatial concentration parameters µ h and µ v which are used in the horizontal and vertical perturbation, respectively. Moreover, the perturbations will be highly oscillating flows, and this oscillation is represented by the spatial oscillation parameters σ h and σ v , which are again used in the horizontal and vertical perturbation, respectively.
Finally we will use intermittent flows. To this end we introduce temporal intermittency functions in section 3.5. These are time-dependent functions which contain the temporal concentration parameters κ h , κ v and temporal oscillation parameters ν h , ν v .

The parameters
As mentioned in section 3.1 we have two sets of four parameters, namely {µ h , σ h , κ h , ν h } and {µ v , σ v , κ v , ν v }, so eight parameters in total. In addition to that, we work with two "master parameters" λ h and λ v , which the other parameters depend on via These exponents are determined in the following Lemma. We will later fix λ h , λ v . These parameters will be very large and such that σ h , σ v ∈ N, as well as κ h , κ v > 1.
Lemma 3.1. Let 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 0 < s 1 , s 3 satisfy the following conditions 8 2 and in addition Proof. We choose 0 < a h , a v < 1 and set By taking logarithms, inequalities (3.3)-(3.7) are equivalent to Using the definition of b v and c v we immediately conclude that (3.10) is valid.
Since these upper bounds for γ i are positive, it remains to show that the upper bounds given by the left-hand sides of (3.8), (3.9), (3.11) and (3.12) are positive as well.
It is obvious that δc v > 0. Furthermore from our choice of can be made small (see section 5.1.2). Moreover (3.7) will be used at several points during the proof. Finally, inequalities (3.4) and (3.6) make sure that the temporal parts of the linear error are controlled, see section 6.3.1.

Inverse divergence operators
Like in most of the convex integration schemes in the context of fluid dynamics in the literature, we will need inverse divergence operators in order to define the new Reynolds stress tensors R h,1 and R v,1 . In this context the first inverse divergence operator goes back to [33]. In this paper we will work with three inverse divergence operators. The horizontal inverse divergence R h and its bilinear version B will be used to define the new horizontal Reynolds stress tensor R h,1 . Those operators are treated in sections 3.3.1 and 3.3.2, respectively. In order to determine the new vertical Reynolds stress tensor R v,1 we need a "vertical inverse divergence" which is just an integral in z. It is introduced in section 3.3.3.

Horizontal inverse divergence
Our horizontal inverse divergence coincides with the two-dimensional inverse divergence from [27]. It is based upon the inverse divergence introduced in [33] and is defined as follows.
The following Lemma, which can also be found in [27,Appendix B], summarizes some properties of the map R h .
1. The following identities hold (3.17) 9 We denote the set of all symmetric 2 × 2 matrices with zero trace by S 2×2 0 . 10 We are using the Einstein summation convention here, in particular we sum over k = 1, 2. Moreover, we recall that δ ij is the Kronecker delta and in the definition of the inverse horizontal Laplacian ∆ −1 h we assume the spatial horizontal average to be zero (in order to ensure uniqueness).
If f is mean-free, i.e.´T 2 f dx = 0, then ) is a Calderón-Zygmund operator, in particular for any 1 < p < ∞ and all A ∈ C ∞ (T 2 ; R 2×2 ) we have For the proof we refer to [27, Appendix B].

Horizontal bilinear inverse divergence
Next we recall the bilinear inverse divergence operator from [27]. For our purposes we call it horizontal bilinear inverse divergence. 20) or written without components (where we have abused notation) We will also use the following Lemma from [27].
Moreover, we have the following estimate The proof of Lemma 3.6 can be found in [27, Appendix B].

Vertical inverse divergence
Finally we introduce the vertical inverse divergence as follows.
Definition 3.7. We define the map 11 R (3.24) 11 We denote the space of all functions in C ∞ (T 3 ; R 2 ) which have zero-mean with respect to z by The vertical inverse divergence operator has the following properties.

The following identities hold for any
(3.27) 2. For 1 ≤ p, q ≤ ∞ and s ∈ R, the operator R v satisfies the following estimates Proof. The identities (3.25), (3.26) are just a simple consequence of the definition of R v . We also observe that i.e. (3.27).
Estimate (3.28) is established simply by moving the L p norm inside the integral. In order to prove estimate (3.29), we first observe that which implies (3.29).
Finally observe that which immediately yields (3.31).

Building blocks for the perturbation
Next we recall the building blocks. We begin with the Mikado flows and Mikado densities which we use to handle R v . We state their existence together with their most important properties in the following proposition. The construction of the Mikado flows and densities is nowadays standard and goes back to [31]. For the proof of the following proposition we refer to [27,Section 4.1].
Proposition 3.9. For each k ∈ {1, 2} there exist functions W k ∈ C ∞ (T 2 ; R 2 ) and φ k ∈ C ∞ (T 2 ; R) (referred to as the Mikado flows and Mikado densities respectively) depending on a parameter µ v , with the following properties: 1. The functions W k , φ k have zero mean for all k ∈ {1, 2}. Moreover where e k denotes the k-th standard basis vector in R 2 , and by construction W k = φ k e k .

For any
3. For all s ≥ 0, 1 ≤ p ≤ ∞ and k, k ′ ∈ {1, 2} with k = k ′ the following estimates hold: Here the implicit constant may depend on s, p, but it does not depend on µ v .
Let us now recall the Mikado flows which we will use to treat R h . In the following proposition  1. The flows W k have zero mean, i.e.´T 2 W k dx = 0, for all k ∈ Λ. Moreover 3. For all s ≥ 0, 1 ≤ p ≤ ∞ and k, k ′ ∈ Λ with k = k ′ the following estimates hold: (3.40) Here the implicit constant may depend on s, p, but it does not depend on µ h . 12 We denote the set of all skew-symmetric 2 × 2 matrices by A 2×2 .

Horizontal temporal intermittency functions
We set g h (t) := κ 1/2 h G(κ h t), more precisely, g h is the 1-periodic extension of the right-hand side (where we require that κ h > 1). Note that from (3.47) we obtain the normalisation identitŷ and furthermore it is straightforward to verify that

Vertical temporal intermittency functions
The vertical temporal oscillation functions are given by the 1-periodic extension (assuming that κ v > 1) of The corresponding temporal correction function is defined by In addition to that we need temporal oscillation functions where the argument of G is shifted. Those are defined as the 1-periodic extension of with corresponding correction function Since G has compact support in (0, 1/2) and κ v > 1, the functions g ± v,1 and g ± v,2 have disjoint supports.
Note that due to the fact that q 2 > 2, we have 1/q 2 < 1 − 1/q 2 which justifies the notation g − v,k , g + v,k for k = 1, 2. Similarly to the horizontal temporal functions which we introduced in section 3.5.1, we have the following estimates for any p ∈ [1, ∞] and k ∈ {1, 2}

Velocity perturbation and new Reynolds stress tensor
In sections 4, 5 and 6 we prove Proposition 2.4 hence we suppose that the assumptions of Proposition 2.4 hold.
The perturbation will be written as where u p,h and u p,v are referred to as the horizontal and vertical principal parts of the perturbation, while u c,h , u c,v , u t,h and u t,v are referred to as the horizontal and vertical spatial and temporal correctors.
Remark 4.1. We would like to remark that u p,h , u c,h and u t,h do no depend on z, while u p,v , u c,v and u t,v are mean-free with respect to z. Therefore, the first three are indeed a barotropic perturbation, while the latter three are a baroclinic perturbation. This is already hidden in (4.1) and (4.2).
In sections 4.1 and 4.3 we define the horizontal perturbation u p and the vertical perturbation u p , respectively. The pressure perturbation P is determined in section 4.2. Finally we define the new Reynolds stress tensors R h,1 and R v,1 in section 4.4.

The horizontal perturbation
We begin by constructing the horizontal perturbation which consists (see above) of a principal part u p,h , a spatial corrector u c,h and a temporal corrector u t,h .
In order to construct u p,h , we introduce a cutoff function χ. First we choose χ ∈ C ∞ ([0, ∞)) to be increasing and satisfying Next we define the function χ(x, t) := χ |R h (x, t)| .
It is straightforward to check that I − R h χ ∈ B 1/2 (I) for all (x, t) ∈ T 2 × [0, T ]. This means in particular that we can evaluate the functions Γ k (see Proposition 3.10) at I − R h χ . We now introduce a temporal smooth cutoff function θ ∈ C ∞ ([0, T ]; [0, 1]) which satisfies in order to achieve the desired property of the supports of the perturbations u p , u p , w p . The horizontal principal perturbation is then defined by where the W k are given by Proposition 3.10 and the amplitude functions are Notice that u p,h does not need to be divergence free. To overcome this we define the corrector u c,h as u c,h : Hence as Ω k is skew-symmetric. Moreover, using the definition of θ in (4.4), we have u p,h = u c,h = 0 whenever dist(t, I c ) ≤ τ /2.
Next, we define the horizontal temporal corrector to be It is straightforward to check that (∇ h ⊗ ∇ h ) : R h is mean-free (so that the inverse Laplacian ∆ −1 h can be applied to this expression), ∇ h · u t,h = 0, and that u t,h = 0 whenever dist(t, I c ) ≤ τ /2.
Finally notice, that u p,h , u c,h and u t,h are indeed independent of z.

The pressure perturbation
The pressure perturbation is defined as follows Note that ∂ z P = 0, since R h , and hence also χ, are independent of z.

The vertical perturbation
We denote the components of the vertical Reynolds stress tensor as R v,k , k = 1, 2, i.e., This allows us to define the vertical principal perturbation as where W k and φ k are given by Proposition 3.9.
We now introduce the vertical spatial corrector in order to make the perturbation divergence free. We set (4.14) Observe that and hence since Ω k is skew-symmetric. Notice that w p,v is independent of z. Again according to the definition of θ in (4.4), Moreover, we introduce the vertical temporal corrector to be Since u t,v does not need to be divergence free, we introduce the corrector It is then simple to check that ∇ h · u t,v + ∂ z w t,v = 0. Similar to u t,h , see above, we have u t,v = 0 and w t,v = 0 whenever dist(t, I c ) ≤ τ /2.
Finally notice, that u p,v , u c,v and u t,v are mean-free with respect to z because R v is mean-free with respect to z.

New Reynolds stress tensors
The goal of this section is to define the new Reynolds stress tensors R h and R v . These will consist of several pieces.

Horizontal oscillation error
Let us first define Proof. Let us first look at the term ∇ h · (u p,h ⊗ u p,h + R h ). We may write Using the definition of the a k , items 1 and 2 of Proposition 3.10 and Lemma 3.6 we find Next we compute Hence we have shown If dist(t, I c ) ≤ τ , then R h = 0 by well-preparedness. If dist(t, I c ) ≥ τ then θ(t) = 1 and hence 1 − θ 2 = 0. This completes the proof of (4.22).

Vertical oscillation error
We define and set Proof. First we observe that since g ± v,1 g ± v,2 = 0, see section 3.5.2. Hence we obtain using Proposition 3.9 Here we used that (1−θ 2 )∂ z R v,k = 0, see the proof of Lemma 4.2. Moreover a straightforward computation shows which finishes the proof of Lemma 4.3.

Linear errors
Next we define the horizontal and vertical linear errors by Note that the arguments of the operators R h and R v satisfy the required properties, i.e. they are independent of z and mean-free with respect to z, respectively.
For convenience let us write

Corrector errors
Finally, we define the horizontal and vertical corrector errors by As in section 4.4.3 we remark that the arguments of the operators R h and R v satisfy the required properties, i.e. they are independent of z and mean-free with respect to z, respectively.

Conclusion
The new Reynolds stress tensors R h,1 , R v,1 are then given by First, we note that by definition R osc,h , R lin,h and R cor,h (and consequently also R h,1 ) are independent of z. Moreover, by definition R osc,v is mean-free with respect to z, and so are R lin,v and R cor,v according to Lemma 3.8 (consequently R v,1 has the same property).
Next we remark that R h,1 (x, t) = R v,1 (x, t) = 0 whenever dist(t, I c ) ≤ τ 2 . Indeed, for the oscillation errors R osc,h , R osc,v this follows from the fact that R h = R v = 0 whenever dist(t, I c ) ≤ τ , and the definition of θ, see (4.4). The fact that R lin,h = 0, R lin,v = 0, R cor,h = 0, and R cor,v = 0 whenever dist(t, I c ) ≤ τ 2 immediately follows from item 1 of Proposition 2.4, which we already proved in section 4.1 and 4.3 Finally, with the help of • the fact that ∇ h · (u + u p + u p ) + ∂ z (w + w p ) = 0, • Lemmas 4.2 and 4.3, • Lemmas 3.4 and 3.8, • the fact that´T 2 u p,h + u c,h dx = 0, see (4.8), • and ∂ z u p,h = ∂ z u c,h = ∂ z u t,h = ∂ z w p,v = 0, see section 4.1 and section 4.3, a long but straightforward computation shows Hence (u + u p + u p , w + w p , p + P, R h,1 , R v,1 ) solves (2.1).

Estimates on the perturbation
In the remaining sections we will use the following convention: for quantities Q 1 and Q 2 we write Q 1 Q 2 if there exists a constant C such that Q 1 ≤ CQ 2 . In general we require that C does not depend on (u, w, p, R h , R v ). However if the right-hand side Q 2 only contains powers of the parameters µ i , σ i , κ i , ν i , λ i (i = h, v), then the implicit constant C may depend on (u, w, p, R h , R v ).

Horizontal principal perturbation
First, we estimate the horizontal part of the principal perturbation. We recall the following Lemma from [27].
Lemma 5.1. We have the following estimates for all n, m ∈ N 0 , p ∈ [1, ∞] The implicit constant in (5.1) might depend on u or R h , whereas the implicit constant in (5.2) neither depends on t nor on u or R h .
For the proof we refer to [27,Lemma 5.2].
With Lemma 5.1 at hand we are ready to prove the estimates on the horizontal principal perturbation.
Lemma 5.2. If λ h is chosen sufficiently large (depending on R h ), then the horizontal principal perturbation satisfies the following estimates Proof. By applying the improved Hölder inequality in Lemma B.1 and Lemmas 3.11 and 5.1, we find that with a constant C u,R h depending on u, R h . Since t → ´T 2 χ(·, x) dx 1/2 is smooth, we can apply Lemma B.1 once again to obtain where we made use of (3.48). Since So have shown that which implies (5.3) by taking λ h sufficiently large (depending on R h ).
Furthermore by applying Lemmas 3.1, 3.11 and 5.1 we find that which proves (5.4).

Vertical principal perturbation
Let us first show the following Lemma.
Proof. Using equation (3.15) we find Next, we observe that according to Lemma 3.11 Then by Lemmas 3.4 and A.3, and inequality (5.6) we obtain Remark 5.4. We present an alternative proof of Lemma 5.3 in the appendix, see section A.3.
Next we estimate the vertical principal perturbation.
Lemma 5.5. If λ v is chosen sufficiently large (depending on R v ), then the vertical principal perturbation satisfies the following inequalities Proof. According to (3.51) and Lemmas 3.1 and 3.11, we obtain So we have shown (5.7), (5.8).
Next, notice that in accordance with Proposition 3.9 and Lemma 3.11 (keeping in mind that s 3 ≤ 1 and hence −s 3 + 1 ≥ 0)
Finally, we derive estimate (5.13) for the product u p,v w p,v . Because g ± v,1 g ± v,2 = 0, see section 3.5.2, and the improved Hölder inequality (Lemma B.1) we have .
First, observe that according to (3.51), (3.52). Next, we estimate where C Rv is a constant depending on R v , and where we used Lemmas 3.11 and A.3. Moreover, we obtain by Lemma 5.3 Hence we have shown which implies (5.13) by choosing λ v sufficiently large (depending on R v ).

Spatial correctors
Lemma 5.7. The spatial correctors satisfy the following estimates Proof. By using Lemmas 3.1, 3.11 and 5.1 as well as estimate (3.51), one gets that

Temporal correctors
Lemma 5.8. The temporal correctors satisfy the estimates where n ∈ N is arbitrary, and the implicit constant may depend on n.

Conclusion
We have already shown in section 4 that u p , u p and w p satisfy as well as item 1 of Proposition 2.4. Hence (u + u p + u p , w + w p ) fulfills (2.3). Additionally, we have shown in section 4 that ∂ z P = 0 and hence p + P satisfies (2.2). Moreover, we proved that (2.1) holds. Consequently (u + u p + u p , w + w p , p + P, R h,1 , R v,1 ) is indeed a solution of the Euler-Reynolds system (2.1)-(2.3). We also showed in section 4 that (u + u p + u p , w + w p , p + P, R h,1 , R v,1 ) is well-prepared for the time interval I and parameter τ /2.
Furthermore, estimates (2.7)-(2.14) of Proposition 2.4 are a simple consequence of Lemmas 5.2, 5.5, 5.7 and 5.8, where one has to choose λ h , λ v sufficiently large, depending on R h and R v , respectively.

Estimates on the stress tensor
In order to finish the proof of Proposition 2.4 it remains to show estimates (2.5), (2.6). These two estimates simply follow from Lemmas 6.1, 6.2, 6.3 and 6.4, which we prove in this section, below.
6.1 Oscillation error 6 Here we have used that (similar to (5.6)) according to Lemma 3.11.
Next, we obtain from (3.50) Finally, using Lemma 5.1 and Proposition 3.10 we get By choosing λ h large enough (depending on R h ), we conclude with (6.1).

Vertical part
Lemma 6.2. If λ v is chosen sufficiently large (depending on R v ), then the vertical oscillation error satisfies Proof. Using (5.17) and Lemma 5.3 we find For the temporal part of the error we obtain by (3.53) Consequently (6.3) follows by choosing λ v sufficiently large, depending on R v .
Finally, according to Lemmas 3.8, 5.2, 5.5, 5.7, 5.8 and A.3 In these estimates we have used the fact that the time interval [0, T ] is finite. Then (6.5) follows by choosing λ h and λ v large enough (again depending on R h and R v , respectively).

Linear error
Lemma 6.4. If λ h and λ v are chosen sufficiently large (depending on R h and R v , respectively), then the linear errors satisfy the estimates In order to prove this Lemma, we consider the time derivative (see section 6.3.1) and advective terms (see section 6.3.2) separately.
Proof of Lemma 6.4. We simply conclude using Lemmas 6.5 and 6.6 below by choosing λ h and λ v large enough.

Advective terms
Lemma 6.6. For the advective part of the linear error, the following bounds hold Proof. Lemmas 3.4, 5.2, 5.5, 5.7 and 5.8 yield For the vertical advective terms we have according to Lemmas 3.8, 5.2, 5.5, 5.7, 5.8 and A.3

The viscous primitive equations
In this section we consider the viscous primitive equations (1.17)- (1.19). We begin by stating the viscous primitive-Reynolds system We prove the following version of Proposition 2.4. Theorem 1.10 can be proven in exactly the same fashion as Theorem 1.3.

The perturbation and Reynolds stress tensors satisfy the following estimates
3. Moreover, we have the following bounds 14) In order to prove Proposition 7.1, we need the following version of Lemma 3.1.
Lemma 7.2. Let 1 ≤ q 1 , q 2 , q 3 ≤ ∞ and 0 < s 1 , s 3 satisfy the conditions (7.4). Then we can Proof. Similar to the proof of Lemma 3.1, it suffices to show that there is a choice of Let us first fix 0 < a h < 1/2, and then b h > 0 such that Note that such a choice is possible since s 1 + 1 − 2 q 1 < 0 according to (7.4). Because (7.4) implies q 1 < 2, (7.24) is equivalent to Moreover as a h < 1/2, we have which are equivalent to (7.18), (7.19) and (7.22) respectively.
Next we choose a v , b v , c v > 0. We simply deduce from (7.4) that Thus we can choose 0 Then we fix which is positive as 1 q 3 − 1 q 2 > 0 which in turn follows from (7.4). The choice of c v immediately implies (7.20), while (7.27) is equivalent to (7.23). Finally (7.4) ensures This is equivalent to which in turn allows for the choice of a small a v > 0 such that (7.21) holds.
Now we can prove Proposition 7.1.
Proof of Proposition 7.1. We make the same choice of perturbations u p , u p , w p , P as we did in the inviscid case, see section 4. The only errors that change compared to the inviscid case are the linear errors, which now contain the additional terms respectively. Thus the validity of (7.8), (7.10)-(7.15) follows immediately from sections 4-6.
As in the proof of (5.4) one can deduce from (7.16) that Analogously we find Similarly we obtain from (7.17) Since u t,h L 1 (W 1,1 ) λ −γ h h and u t,v L 1 (W 1,1 ) λ −γv v , which follow immediately from Lemma 5.8, we deduce (7.7) and (7.9) by choosing λ h , λ v sufficiently large (depending on R h and R v , respectively).
In order to show (7.5) and (7.6), it remains to estimate the terms in (7.28). Using Lemmas 3.4, 3.8 and A.3, as well as u p L 1 (W 1,1 ) λ −γ h h and u p L 1 (W 1,1 ) λ −γv v which we established above, we find

Two-dimensional hydrostatic Euler equations
In this section we will develop a convex integration scheme for the two-dimensional hydrostatic Euler equations (1.21)-(1.23), which is somewhat different in nature to the scheme in the three-dimensional case. A similar scheme will be established in section 9 for the (two-dimensional) Prandtl equations (1.25)-(1.27).
In two dimensions the hydrostatic Euler-Reynolds system (2.1)-(2.3) reduces to We observe that in contrast to the three-dimensional case u, R h and R v are now just scalar quantities, where R h does not depend on z and R v is mean-free with respect to z. The former allows to include R h as part of the pressure. In other words by setting we may assume without loss of generality that R h = 0 (up to a redefinition of the pressure). So all in all the two-dimensional hydrostatic Euler-Reynolds system we will work with, reads with unknowns u, w, p and R v . As R h is no longer there, the only task is to minimize R v . For this reason we will only have a baroclinic perturbation u p in Proposition 8.1 below.
In this section, we will prove the following version of the inductive proposition (cf. Proposition 2.4). Theorem 1.18 then follows exactly as Theorem 1.3. satisfy the constraints in (1.24). Finally, let δ, ǫ > 0 be arbitrary. Then there exists another smooth solution (u+ u p , w +w p , p+P, R v,1 ) of the two-dimensional hydrostatic Euler-Reynolds system (8.1)-(8.3) which is well-prepared with respect to the same time interval I and parameter τ /2, and has the following properties:

It satisfies the following estimates
3. Moreover, we have that

Preliminaries
Due to the fact that there is no longer an error R h in the two-dimensional hydrostatic Euler  (8.11) in addition to (3.5), (3.7) and Proof. Similar to the proof of Lemma 3.1, it suffices to show that there is a choice of We know from (1.24) that By setting the latter implies that which in turn allows for the choice of a small a v > 0 such that (8.13) holds. The definition of c v immediately implies (8.12). Finally (1.24) guarantees that which is equivalent to (8.14).
Remark 8.3. Notice that compared to Lemma 3.1 we have replaced (3.6) by (8.11), where the latter is a weaker restriction than the former. Indeed it is simple to see that (3.6) implies (8.11) provided q 2 > 2. Note furthermore that (8.11) suffices to prove Lemma 6.5. Moreover, the additional inequality (8.10) is needed to deal with an additional spatial corrector for the vertical velocity.
Moreover, we need a two-dimensional version of the vertical inverse divergence operator R v , cf. Definition 3.7.
Note that the vertical inverse divergence defined in Definition 8.4 has the same properties as stated in Lemma 3.8.
Regarding the building blocks we will use the following version of Proposition 3.9. 15 Again we denote the space of all functions in C ∞ (T 2 ; R) which have zero-mean with respect to z by C ∞ 0,z (T 2 ; R).
Then we introduce a corrector for the vertical velocity in order to get a divergence-free perturbation: where R v is now given by Definition 8.4. Note that R v is mean-free with respect to z and hence the operator R v can be applied to the expression in (8.20). With Lemma 3.8 it is simple to see that ∂ x 1 u p,v + ∂ z w p,v = 0.
Finally, we introduce a temporal corrector of the form To keep the whole velocity field divergence-free, we now must introduce a temporal corrector for the vertical velocity We observe that ∂ x 1 u p + ∂ z w p = 0, item 1 of Proposition 8.1 holds and u p,v and u t,v are indeed mean-free with respect to z.

The new Reynolds stress tensor
As in the three-dimensional scheme, the new Reynolds stress tensor R v,1 will be written as Exactly as in Lemma 4.3 one can show that Next, we define the linear and corrector errors by Note that the arguments of the operator R v are indeed mean-free with respect to z.
In the three-dimensional convex integration scheme we introduced a pressure perturbation, see section 4.2. An analogue of this perturbation is not needed in two dimensions as there is no barotropic perturbation. However, the pressure has to absorb some terms that were covered by the horizontal Reynolds stress tensor in the three-dimensional scheme. To this end we define We observe that ∂ z P = 0.
Let us finally remark that it is straightforward to see that R v,1 is mean-free with respect to z, that R v,1 (x, t) = 0 whenever dist(t, I c ) ≤ τ 2 , and that (u + u p , w + w p , p + P, R v,1 ) solves (8.1).

Estimates on the perturbation
Now we claim the following estimates on the perturbation. Lemma 8.7. If λ v is chosen sufficiently large (depending on R v ), then we have that Proof. In fact the proof of (8.25)-(8.28) can be taken verbatim from Lemma 5.5. Estimate (8.29) can also be proven in a similar way as in Lemma 5.5. To this end we need a version of Lemma 5.3, namely the estimate Consequently R h has the properties stated in Lemma 3.8. Using property (3.30) we are able to prove (8.30) and thus (8.29) follows.
The spatial and temporal correctors can be estimated as follows.
Lemma 8.8. The spatial and temporal correctors satisfy the following estimates where n ∈ N is arbitrary, and the implicit constant may depend on n.
Proof. We only prove (8.31), as the proof of the estimates for the temporal correctors (8.32), (8.33) is similar to the proof of Lemma 5.8. We obtain using Lemmas 3.8 and 8.2, and equations (3.51), (8.16) Analogously Here we have used that R v is bounded in H −s 3 , see (3.29) and keeping in mind that H −s 3 = B −s 3 2,2 , and s 3 ≤ 1.

Estimates on the Reynolds stress tensor
In order to finish the proof of Proposition 8.1, it remains to show (8.4).
Lemma 8.9. If λ v is chosen sufficiently large (depending on R v ), then the errors satisfy the following estimates Proof. Lemma 8.9 can be proven similarly to Lemmas 6.2-6.6, with only small modifications: To obtain the required estimate for R osc,x,v one has to use (8.30) rather than Lemma 5.3. Moreover, since in two dimensions there is no spatial corrector u c,v , the time derivative part of the linear error R lin,t,v must be estimated slightly differently compared to Lemma 6.5. To this end we write Hence we find (by using inequality (8.11)) This finishes the proof of Proposition 8.1.

The two-dimensional Prandtl equations
In this section we will study the following Prandtl-Reynolds system with unknowns (u, w, p, R v ). As in section 8 there is no horizontal Reynolds stress tensor R h . In this setting we have the following version of the inductive proposition (cf. Proposition 2.4). As before, the proof of Theorem 1.23 then works in the same way as the proof of Theorem 1.3.

The perturbation and the new Reynolds stress tensor satisfy the following estimates
3. Finally, the products of the vertical and horizontal perturbations satisfy that Proof. In order to prove Proposition 9.1 we modify the proof of Proposition 8.1 in the same way as we did in the three-dimensional case, cf. proof of Proposition 7.1. Specifically we choose u p , w p and P as in the proof of Proposition 8.1, while the linear error now contains the additional term Then it follows from section 8 that estimates (9.6)-(9.10) hold. In order to show (9.5) we proceed as in the proof of Proposition 7.1. To this end we need the additional parameter estimate We can achieve this as soon as we see that (9.12) holds according to (1.28).
It remains to estimate the additional term in (9.11). This works exactly as in the proof of Proposition 7.1.

A.1 Littlewood-Paley theory and Besov spaces
We first introduce a dyadic partition of unity {ρ j } ∞ j=−1 as follows ρ 0 (ξ) = ρ(ξ), ρ j (ξ) = ρ(2 −j ξ) for j = 1, 2, . . . , ρ −1 (ξ) = 1 − It is also possible to define Besov spaces using difference quotients. We first define the forward finite difference operator . The higher order finite differences are defined inductively. For m ≥ 2, we define that Let s > 0 and 1 ≤ p, q ≤ ∞ and let ⌊s⌋ denote the integer part of s, then Besov norm may be defined as follows (if q < ∞) If q = ∞, the Besov norm is given by These two different definitions of the Besov norm are equivalent. Remark A.1. Note that the index h in the notation for the finite differences ∆ m h represents the shift. This is in contrast to the main body of this paper where the index h always means "horizontal". Remark A.2. The reader should note that B s p,p (T 3 ) = W s,p (T 3 ) when s ∈ R\Z and 1 ≤ p ≤ ∞. This is stated in [1,Equation 3.5] for example. The equivalence of function spaces in the case 0 < s < 1 can also be found in [88,Definition 32.2] (when the interpolation space definition of Besov spaces is used) and in [86,Proposition 2] and [8, Page 1686] (where the Sobolev spaces are defined using the Sobolev-Slobodeckij seminorm). We note that the case 0 < s < 1 can easily be extended to all s > 0 with s / ∈ N by using Theorem 2.3 in [83]. Finally, we recall that B s 2,2 (T 3 ) = H s (T 3 ) for all s ∈ R (even when s is an integer), see [3,Page 99]. The latter is used several times in this paper.

A.2 Paradifferential calculus
We recall Bony's product decomposition f g = T f g + T g f + R(f, g). (A.5) The terms T f g and T g f are called paraproducts and are given by The term R(f, g) is referred to as the resonance term and is defined as • If α < 0 then for any f ∈ B α A.3 Alternative proof of Lemma 5.3 In this section we present an alternative proof of Lemma 5.3 using a paraproduct decompostion and Lemmas A.4 and A.5.
Alternative Proof of Lemma 5.3. We start by decomposing the term under consideration in terms of Bony's decomposition, i.e. .
Using Lemma A.4 we find Another application of Lemma A.4 yields

Finally Lemmas A.3 and A.5 lead to
where 0 < s ≪ 1 can be chosen arbitrary.
To conclude we proceed similar to the original proof. Due Lemmas 3.4 and A.3, and estimate (5.6) This finishes the proof of Lemma 5.3.

B Other estimates B.1 Improved Hölder inequality
Let us recall the following estimate from [68].

B.2 Oscillatory paraproduct estimate
We are going to prove a version of Lemma B.1 in the case of Besov spaces.
Proof. We first recall the low-frequency cut-off operator S j from [3]: Hence we may write T (f σ , g) = ∞ j=−1 S j−1 f σ ∆ j g = ∞ j=1 S j−1 f σ ∆ j g, where the latter equation follows from the fact that S −2 f σ = S −1 f σ = 0. In order to estimate the Besov norm of T (f σ , g), we will use [3,Lemma 2.69]. To be able to use this lemma we need to show that supp F (S j−1 f σ ∆ j g) lies in 2 j C for any j ∈ N 0 , where C is a fixed annulus.
In order to show this, let j ∈ N and i ∈ {−1, ..., j − 2}. By construction of the dyadic partition of unity, there exist radii 0 < r < r 0 < R such that the support of ρ −1 is contained in the ball with radius r 0 , the support of ρ 0 is contained in the annulus with inner radius r and outer radius R, and r 0 < 2r < R < 4r. Next we observe supp F (∆ i f σ ∆ j g) = supp (ρ i f σ ) * (ρ j g) .
For i = −1 this yields that supp F (∆ i f σ ∆ j g) is contained in the annulus with inner radius 2 j (r − 2 −j r 0 ) and outer radius 2 j (R + 2 −j r 0 ), which is in turn a subset of the annulus with inner and outer radii 2 j (r − 1 2 r 0 ) and 2 j (R + 1 2 r 0 ). Note that the inner radius is positive due to 2r > r 0 . Similarly for i ≥ 0 we obtain that supp F (∆ i f σ ∆ j g) is contained in the annulus with inner radius 2 j (r − 2 i−j R) and outer radius 2 j (R + 2 i−j R), which is in turn subset of the annulus with inner and outer radii 2 j (r − 1 4 R) and 2 j (R + 1 4 R). Again note that 4r > R implies that the inner radius is positive. Hence there exists an annulus C such that supp F (S j−1 f σ ∆ j g) ⊂ 2 j C.
Hence we may apply Lemma 2.69 from [3] to conclude that T (f σ , g) B s p,q 2 js S j−1 f σ ∆ j g L p j∈N l q (N) .
In order to estimate S j−1 f σ ∆ j g L p we define φ i := F −1 ρ i . Hence we have ∆ i f = φ i * f . A direct computation shows φ i L 1 = φ 0 L 1 for any i ∈ N. Hence we obtain for 1 ≤ p < ∞ and for any i, j by Minkowski's integral inequality and Lemma B.1 ˆR n φ i (y)f σ (x − y) dy p dx 1/p ≤ˆR n ˆR n |f σ (x − y)φ i (y)∆ j g(x)| p dx 1/p dy =ˆR n |φ i (y)| f (·σ − σy)∆ j g(·) L p dy ˆR n |φ i (y)| f (· − σy) L p ∆ j g L p + σ −1/p ∆ j g C 1 f (· − σy) L p dy For the case p = ∞ one obtains the same result, the details are left to the reader. Hence for any j ∈ N. Combining (B.3) and (B.4) we get for any 1 ≤ q < ∞ where we have used that ∆ j g C 1 ≤ ∇∆ j g L ∞ + ∆ j g L ∞ , ∇∆ j g = ∆ j ∇g and Lemma A.3. For the case q = ∞ we proceed analogously.
Lemma B.2 can be used to prove the following slightly weaker version of Lemma 5.3.
Lemma B.3. For all 1 ≤ p ≤ ∞, δ > 0 and k ∈ {1, 2} we have Proof. Compared to the proof presented in section A.3 the only difference is how we handle the paraproduct in (A.9). We use Lemma B.2 and estimate (5.6) to find

With Lemma B.3 one can show
which is slightly weaker than (5.13). This finally allows to prove a version of Theorem 1.3 where the regularity parameter is s = 1+.