{'fraction_non_alphanumeric': 0.05390060463058546, 'fraction_numerical': 0.021250553015779383, 'mean_word_length': 4.047714753610243, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vortices in stably-stratified rapidly rotating Boussinesq convection

We study the Boussinesq approximation for rapidly rotating stably-stratified fluids in a three dimensional infinite layer with either stress-free or periodic boundary conditions in the vertical direction. For initial conditions satisfying a certain quasi-geostrophic smallness condition, we use dispersive estimates and the large rotation limit to prove global-in-time existence of solutions. We then use self-similar variable techniques to show that the barotropic vorticity converges to an Oseen vortex, while other components decay to zero. We finally use algebraically weighted spaces to determine leading order asymptotics. In particular we show that the barotropic vorticity approaches the Oseen vortex with algebraic rate while the barotropic vertical velocity and thermal fluctuations go to zero as Gaussians whose amplitudes oscillate in opposite phase of each other while decaying with an algebraic rate.


Background
The rotating Boussinesq equations have been widely used to study the effects of rotation and density stratification on flow dynamics in the Earth's atmosphere and oceans, as well as in many other geophysical settings. These equations are an approximation of the compressible, rotating Navier-Stokes equations where one assumes weak density variations in only the equation of state and the buoyancy term. Furthermore, one assumes that the density varies linearly with temperature, and, if one is studying ocean dynamics, salinity.
After such approximations and nondimensionalization, these equations take the form of the incompressible Navier-Stokes equation posed in a non-inertial rotating frame coupled with an equation for temperature fluctuations in the stratified fluid, Here u describes the velocity field of the fluid, g > 0 is a gravitational constant, ν > 0 represents the kinematic viscosity, and ν ′ > 0 the thermal viscosity. The term Ωe 3 × u arises from the effect of the Coriolis force due to the rotating frame. The angular frequency Ω is also proportional to the inverse of the Rossby number which measures the rate of rotation relative to the characteristic length scale of the fluid. While Ω is in general x-dependent (say for spherical geometry in the case of the atmosphere of the Earth), it suffices in many situations to use the "f-plane" approximation where Ω is constant. We will assume this for the rest of the work. Next θ and p denote thermal and pressure fluctuations about a horizontally homogeneous mean state. We assume that the mean velocity is zero, the mean density varies linearity with x 3 , and that mean density and mean pressure are in hydrostatic balance, where the pressure gradient is balanced by the force of gravity on the fluid, This approximation arises from the fact than in typical applications of this model, such as the Earth's oceans, the horizontal scale is much larger than than the vertical scale; see [41] for more explanation. If dρ dx3 > 0, so that the fluid is convective with less dense fluid lying below more dense fluid (say the fluid is heated from below for example), the system is said to be unstably-stratified. See more Section 1.4 for more discussion.
In this work we focus on the case dρ dx3 < 0, where less dense fluid is above more dense fluid, known as stable stratification. The Brunt-Väisälä buoyancy frequency N := −dρ/dx 3 denotes the frequency of the neutral oscillation of a small "parcel" of vertically unstable fluid feeling the effects of buoyancy and gravity inside a uniform background of stably-stratified fluid; see [38, §1.3] for more discussion on this. Substituting these quantities into (1.1) we obtain For a full derivation of this system see [21]. For general discussion on these equations, the many interesting physical phenomenon they model, and its various approximating limits, see [41,43,21]. For a more mathematically focused introduction and review of the subject see [38,35,14,33,39].

Overview of mathematical literature
Due to the wealth of physical applications and the historical inaccessibility of various asymptotic parameter limits in numerical simulations, a sizeable body of mathematical research has been performed on this system and its many variants, with many interesting analytical tools brought to bear on the problem. In addition, the fact that dynamical systems theory has often proven a valuable tool in analyzing the behavior of these equations means that they have featured in Nonlinearity's pages since its earliest days. For instance, papers [34] and [35] both treat models very similar to that treated here, and discuss their applicability to questions of ocean and atmospheric dynamics. In order to discuss some of the key issues in the analysis of this system we set ν = ν ′ = 1; see Remark 1.3 for discussion on the ν = ν ′ and ν = ν ′ = 1 cases. Following [24,3,28], to obtain a skew-symmetric operator we scale θ → √ g N θ, and denote Γ = N √ g. We also combine all the dependent variables into a single-four vector v = (u 1 , u 2 .u 3 , θ). Finally, applying the Helmholtz projection P defined in (2.6) below, we obtain from (1.4) the formulation ∂ t v + P(u ·∇)v − diag(∆)v + PJ Ω,Γ Pv = 0, divu = 0 (1.7) with∇ = (∂ 1 , ∂ 2 , ∂ 3 , 0) T , and PJ Ω,Γ P is the operator formed by the Coriolis, buoyancy, and gravity effects with, J Ω,Γ = ΩJ 0 0 ΓJ , J = 0 −1 1 0 .
As discussed in [28] (and reviewed below), there is one physically relevant eigendirection, which we call the "quasi-geostrophic mode", which undergoes no dispersive smoothing. However, the other eigendirections correspond to rapidly oscillating waves known as "inertial" or "Poincaré waves". While these do not decay in L 2 they do decay in L p spaces with p > 2. This can be quantified with the aid of Strichartz estimates. One then attempts to control such oscillatory behavior and determine the long-time, quasi-geostrophic dynamics of the system. Such dynamics are briefly discussed below. See the book by Chemin et. al [14] for more discussion on the role of inertial wave phenomenon in the context of a rotating incompressible Navier-Stokes system without density effects. There they also illustrate interesting connections of the Ω → ∞ limit with the incompressible limit (low mach number) of compressible Navier-Stokes. Also see the text of Majda which discusses these topics in the context of invisid shallow-water wave equations [38].
Global well-posedness and regularity has been one of the main focuses of mathematical research on this system. Here results date back to the works of Lions, Temam, and Wang [34,35] which proved global existence for weak solutions and estimates on the dimensions of global attractors, for a related model, known as the "primitive" equations, which focused on ocean dynamics. This model includes an equation for water salinity and considered the system in a thin layer on a sphere, S 2 × (−H, 0) where the depth, H was dependent on the spherical variables. The primitive equations are obtained from (1.4) under the hydostatic approximation, where one uses the disparity between the vertical and horizontal length scales in the Earth's ocean, to replace the equation for u 3 with the leading order balance ∂p ∂x3 = −θg. Subsequent works (see [6,50,27,23] and references therein) generalized these results, extending them to strong solutions (both small-and global-time existence) and to more general domains.
The works of Babin, Mahalov, and Nicolenko [2,3,4] considered systems of the form (1.4) with an additional forcing term, posed on fully periodic domains x ∈ [0, 2πa 1 ]×[0, 2πa 2 ]×[0, 2πa 3 ] or horizontally periodic domains with "stress-free" boundary conditions in the vertical direction: (1.8) There, inertial waves cannot "escape to infinity" due to the bounded domain, possibly inducing resonant three dimensional wave interactions. Small divisor techniques are required to prove global existence and regularity of strong solutions, along with the existence of attractors, for certain non-resonant domains in this setting. In these works, they also studied how, in the joint limit Ω = 1/ǫ, N 2 =Ñ /ǫ, ǫ → 0 solutions of (1.4) can be uniformly approximated in time by solutions of the quasi-geostrophic equations, a classic and often-used model first derived by Charney [8] for the slow dynamics of the atmosphere. There also has been an intense study of the full Boussinesq system (1.4) in the whole space x ∈ R 3 . In a series of works [12,13,11,10], Charve has shown convergence to the quasi-geostrophic system in the whole space x ∈ R 3 for both weak and strong solutions of various regularity and initial condition size, in the same limit ǫ → 0 mentioned above. The idea underlying all these approaches is to decompose the solution as a part governed by the quasi-geostrophic system and a "remainder" and to show that the oscillatory effects arising from the fast rotation allow one to prove existence and uniqueness of solutions in three-dimensions, even for large initial data.
Subsequent results such as [31] extend the use of dispersive estimates in R 3 to other L p spaces. The most recent result [28] investigates more deeply the connections between harmonic analysis and the principal curvatures of the linear dispersion relation p η (k); see [22] for example. In particular, degeneracies of the Hessian, (∂ 2 kikj p η ) ij , of the dispersion relation are used to study more general initial conditions, only requiring smallness in the quasi-geostrophic component. We remark that because of the smallness of the quasi-geostrophic initial condition these last two results omit study of the quasigeostrophic limit as in Charve's work. See [28] for a more detailed discussion of these topics and review of the related literature.

The role of dynamical systems
Dynamical systems has also been integral to the study geophysical fluid dynamics. In fact the study of such fluid models has been one of the main motivations for the development of many mathematical tools such as global attractors, inertial manifolds, and slow manifolds. Again, Nonlinearity was a leader in this work and already, in the second volume, dynamical systems methods were exploited to study the stability of stratified fluids [25].
Slow manifolds During the 1970's and 1980's many researchers, including Lorenz [36], posited the existence of a slow invariant manifold in the Boussinesq system which controlled and organized the dynamics. In an effort to understand the interplay between fast, oscillatory gravity-waves and the slow geostrophic dynamics for Ω, Γ ∼ O(1/ǫ) with small ǫ, researchers worked to discover if there was a manifold in the (infinite-dimensional) phase space, which was invariant under the time evolution, devoid of fast waves, and contained the slowly evolving geostrophic dynamics which governs much of the long term behavior of the system. Such manifolds would then be of use as they could give information on the global attractor of the system. While the current consensus is that no such exact manifold exists [37], there has been rigorous work characterizing approximately invariant, or "fuzzy", manifolds when the Boussinesq system undergoes certain types of forcing [49,48,42]. These works have shown that in finite-time the fast dynamics decay and stay O(ǫ)-small, so that the dynamics stay near geostrophic balance. For a nice review of this topic, the important work in the area, and some of the most pressing open questions, see [49].
Scaling variables, invariant manifolds, and omega-limit sets In work that is the closest to our own in both technique and aims, Gallay and Roussier-Michon have used dynamical systems techniques to study long-time asymptotics of the rotating incompressible Navier-Stokes system In [17], they considered the above system posed on R 2 × [0, 1] with periodic boundary conditions in the vertical direction. They showed, for any size initial condition and a correspondingly large Ω, the existence of global, infinite energy solutions which converge to the x 3 -independent diffusively decaying Oseen vortex solution. In other words, they showed for high rotation rates, solutions asymptotically behave like the non-rotating two-dimensional Navier-Stokes equations. Using the barotropic-baroclinic decomposition described in [14], which splits the vector field into an x 3 -independent component,ū, and an x 3 -mean-zero component,ũ, one obtains a system of three equations, one each for the horizontal vorticityω 3 = ∂ 1ū2 − ∂ 2ū1 = 1 0 (∇ × u) 3 dx 3 , the vertical barotropic velocityū 3 , and the baroclinic velocityũ. They obtain exponential decay of the baroclinic component using Poincare's inequality (due to vertical mean-zero and boundedness of vertical domain) and Strichartz estimates. After using linear convection-diffusion estimates of Carlen and Loss [7] to show the algebraic decay of ∇ū 3 , they then study the barotropic vorticity which satisfies an equation very similar to the vorticity formulation of two-dimensional Navier Stokes. Taking the approach of Gallay and Wayne [19], they can employ scaling variables and compactness arguments to determine the omega-limit set of the system, showing thatω, and hence u itself, asymptotically converge to the Oseen vortex solution as t → ∞. However, they do not consider the effects of the coupling between the velocity field and the temperature, nor do they derive detailed asymptotics for a solution near a vortex, like those in Theorem 1.

Un-stably stratified fluids
The un-stably stratified system has also been the focus of much recent study. Here the Prandtl number ν/ν ′ is assumed to be one for simplicity, and R is the Rayleigh number, measuring the effects of gravity, thermal expansion, and the density difference in the background linear density profile. Physical boundary conditions once again take the form (1.8). Such models typically arise when studying more extreme systems, such as rapidly-rotating convective atmospheric layers and are contained in the more complicated magneto-hydrodynamic models used to study stellar plasmas. See [32,5] for a review of these applications. The convective nature of the unstable stratification leads to a rich family of structures and dynamics for various parameter ranges. Experimental and numerical studies of (1.10), with various geometries and boundary conditions, have exhibited a wealth of interesting behaviors, such as convection cells, localized plumes, and large-scale turbulent vortices as the Rayleigh number is varied. See the introduction of [16] for a nice review of the numerical and experimental literature and more discussion on this topic.
In the case of rapid rotation, Ω ≫ 1, the Coriolis term suppresses behaviors characteristic of threedimensional incompressible fluids, such as direct turbulence cascades. This causes the system to behave similarly to a forced two-dimensional fluid, where an inverse cascade dominates turbulent dynamics; see once again [32,5]. Assuming large enough Rayleigh number, R, the thermal forcing creates small-scale eddies which become vertically "aligned" and coalesce into large domain-scale turbulent vortices.
Mathematically, much less is known about these equations compared with the stably-stratified case. While there are some results on existence and bifurcation of finite dimensional attractors in fully periodic domains T 3 [33,26], little work has been done to rigorously characterize coherent structures which arise in these systems. It would be interesting to determine if an invariant or approximately invariant slow-manifold existed in the system. Given the previous literature and discussion above, a reasonable candidate for a slow variable is the the barotropic component of the full three-dimensional velocity, The boundary conditions on θ and u 3 immediately give thatū 3 =θ = 0. Then denoting the projection associated with this decomposition as Q, the equation for the barotropic componentv = (ū h , 0, 0) T can readily be found to be whereṽ = (1 − Q)v is the baroclinic part of the system. Thusū h satisfies a forced 2-dimensional Navier-Stokes equation with forcing term Q [(ũ h · ∇ h )ũ h ], coming from the baroclinic system (which we have not written down here). One would hope to show that in the high rotation limit, Ω = ǫ −1 , 0 < ǫ ≪ 1 (possibly also scaling the Rayleigh number R =R/ǫ as well), the Coriolis force overpowers the unstable stratification term causing the baroclinic variablesṽ to decay to O(ǫ) sizes after a finite-time and remain small for all subsequent times, as in [49]. Hence, the barotropic subspace {ṽ = 0} would form an approximately invariant manifold for the full dynamics. One would then hope to characterize the dynamics ofū h , either explicitly characterizing large-scale vortices or at the very least showing the existence of, and characterizing the attractor. As forced two-dimensional turbulence is in general only understood statistically (see [5]) this approach seems to be an interesting line of research, which we shall pursue in the future.
The system (1.10) is in general difficult to simulate, with direct numerical simulations for large Rayleigh number only being performed in the last few years; see for example [16]. Julien and Knobloch with various collaborators have derived a formal asymptotic quasi-geostrophic model which eliminates fast inertial waves and certain boundary effects, known as Ekman layers, and is thus much more tractable numerically while still retaining much of the three-dimensional dynamics (as opposed to a 2-D Navier-Stokes system with arbitrary forcing). In a series of works [47,29,45], they have shown this model exhibits various coherent structures and similar statistical behaviors observed in experiments and recent simulations of the full system (1.10). It would be of great interest to study such coherent structures of in this asymptotic model rigorously and consequently to rigorously study the asymptotic convergence of (1.10) to this formal model in a way similar to that of the quasi-geostrophic approximation of the stably-stratified system discussed above.

Our results
This paper focuses on determining leading order asymptotics of (1.7). We consider the stably-stratified Boussinesq system (1.7) above, posed on D = R 2 × [0, 1] with either stress-free boundary conditions 13) or periodic boundary conditions (1.14) While the latter boundary conditions are for the most part non-physical, they have been often studied as an idealized version of the system.

Local dynamics near a barotropic Oseen vortex
Our main result of this paper (proved in Section 3), is to determine the leading order asymptotics in the case of periodic vertical boundary conditions. As in many previous works, we begin by splitting the evolution into its barotropic (i.e. vertically averaged) and baroclinic parts,v = Qv : The purely barotropic evolution is obviously an invariant manifold within the phase space of the full system. We then identify a family of explicit vortex solutions in the barotropic system. These vortices correspond to an Oseen vortex for the two-dimensional Navier-Stokes equation in the two horizontal components of the velocity, and a coupled pair of vortices in the vertical velocity and temperature fields which oscillate with the Brunt-Väisälä frequency. We then show that regardless of the vortex amplitude, these solutions are stable with respect to the full Boussinesq evolution. From a dynamical systems perspective, this shows that these solutions are in the barotropic manifold and are at least locally attractive. We also note, that for this local analysis, we need not assume the rotation rate is large, even to obtain the stability of large vortex solutions.
Corollary 1.1. One can use the Biot-Savart relationships in Proposition 3.4 to readily conclude the following decay estimates forv, 1/q = 1/p − 1/2 with p ∈ [1,2], and t > 0, Global existence and asymptotics To study more general initial data, one can then use a high rotation number limit, |Ω| ≫ 1, in dispersive estimates and energy-type estimates to obtain global existence and asymptotics, albeit with less information on the rate of convergence. In particular, requiring smallness of only the quasi-geostrophic part of the initial condition, that is initial data v 0 satisfying R 3 | ṽ(k, 0), a g (k) C 4 | 2 dk ≪ 1, with a g (k) defined in (2.12), we prove global existence of mild solutions. We also once again show that at leading order, thermal fluctuations decay to zero and the velocity field converges to the two dimensional vortex solution with u 0 h defined above. Depending on the boundary conditions, we consider solutions in one of the two the Banach spaces with the norms 20) where once againv = Qv := 1 0 v(x h , x 3 )dx 3 , andṽ = (1 − Q)ṽ. Also, we define a projection operator S on L 2 as where F is the Fourier transform, defined on L 2 (D). Our results for these two boundary conditions are as follows: Theorem 2. For all Γ ∈ R {0}, and initial conditions v 0 ∈ X i with either i = p or i = sf, and Sv 0 L 2 sufficiently small, there exists Ω 0 > 0 such that, for all |Ω| ≥ Ω 0 , the system (1.7) with either periodic boundary conditions (1.14) or stress-free boundary conditions (1.13) respectively, has a mild solution v ∈ C 0 ([0, ∞), X i ) satisfying v(·, 0) = v 0 . Furthermore, there exists a C > 0 such that ||v(t)|| Xi ≤ C for all t > 0 and given A = D (curl u 0 ) 3 dx, this solution satisfies , this term has non-zero X i norm as t → ∞. However, the remaining terms in the asymptotics in Theorem 1 all vanish in the X i -norm -they are "invisible" in this theorem. Thus, while Theorem 2 has the advantage of treating global initial data, it gives far less detailed information about the long-time behavior of solutions than Theorem 1.
The proof of Theorem 2 for each type of boundary condition follows the approach of [17] with modifications to account for the inclusion of temperature effects. Namely, we use a barotropic/baroclinic decomposition to split v into a x 3 independent partv and a x 3 -mean zero partṽ. We shall further decompose the barotropic vectorv studying a system of equations forω 3 = (∇ ×ū) 3 ,ū 3 , andθ. Note stress-free boundary conditions forceū 3 ≡θ ≡ 0 identically.
We use dispersive estimates to showṽ decays exponentially fast and then use energy methods and Gronwall's inequality to prove global existence. Diffusive estimates and dynamical systems techniques can then be used to determine the leading order asymptotics ofū h , the barotropic part, hence determining the leading order dynamics of the system. Since the proofs of the two cases in Theorem 2 are quite similar, we only give the proof in the case of stress-free boundary conditions. Note we do not obtain as detailed information on the decay rates of different solution components compared to that of Theorem 1.

Remark 1.3.
Qualitatively similar results should hold for the more general situation of differing viscosities ν ′ = ν. In this situation, while the expressions for the eigenvalues and eigenvectors of the linear system in Fourier space are more complicated and the collection of eigenvectors are not orthogonal, three of these vectors are still orthogonal to the Fourier vector (k 1 , k 2 , k 3 , 0) T . Thus, one should be able to push through the required dispersive estimates, as done in [9, §4] for example, and obtain the same asymptotics forṽ up to an O(Ω −1 )-sized correction. The changes we expect if ν = ν ′ = 1 are also slightly more complicated. On the whole space, one can change the length scale so that the viscosity is always equal to one. However, if we make such a rescaling in our situation, it would change the thickness of the fluid layer. A different value of ν would probably result in a different decay rate for the baroclinic components of the solution, but we still expect on heuristic grounds that they would decay more rapidly than any inverse power of t. For the barotropic components of the solution, we expect the decay rates would be the same as those in Theorem 1, since these parts of the solution depend only on the unbounded variables (x 1 , x 2 ).
While such results seem to have been implicitly known in the literature, (namely that for sufficiently small quasi-geostrophic part the fast-oscillatory component decays as t → +∞, or in other words inertial waves escape to infinity, and the system asymptotes to solutions of 2-D Navier-Stokes equations), to our knowledge it has never been rigorously stated or proven, nor has the detailed asymptotic behavior obtained in Theorem 1 been derived, especially not for perturbations of large, barotropic vortices. Furthermore, we believe the view point taken here could be useful in future studies, especially in the study of unstably-stratified convective rotating systems where a wealth of coherent structures arise; see Section 1.4.
A natural next line of study extending from this work would be to investigate how much one can say about asymptotics for fully arbitrary initial conditions v 0 ∈ X, where there is no restriction on the quasi-geostrophic component. As mentioned above, Charve has shown in R 3 that general initial data evolves towards solutions of the quasi-geostrophic equation. Thus, a useful question to frame such an inquiry could be: "Can dynamical systems techniques allow for a more refined characterization and classification of the asymptotic dynamics of the stably-stratified system with realistic boundary conditions?" Finally, we believe it would of great interest to investigate whether dynamical systems techniques could be used to characterize the physical phenomenon observed in unstably-stratified rotating convection, described in Section 1.4 above. This seems natural as such approaches have been very successful in describing coherent structures in other contexts (see for example [30]), so one expects that they could be successful for the cellular and plume like dynamics of the rotating system. Furthermore, it would be of interest whether the scaling variables/invariant manifold approach of [18,19] could be used to describe the large-scale turbulent vortices which arise and grow out of small-scale turbulent eddies.
Outline of rest of work The rest of the work is organized as follows. In Section 2, we collect some facts about (1.7) with both periodic and stress-free boundary conditions and set-up our framework. In Section 3, we prove our main result Theorem 1, leaving proofs of several technical propositions to Appendix A. In Section 4.1 we collect estimates on the linear dispersive equation associated with our system. We then use this information in Section 4.2 to prove global existence in the stress-free case, while in Section 4.3 we prove the leading order asymptotics under our assumptions, completing the proof of Theorem 2.

Preliminaries
To recall, we will study the system We denote v = (u 1 , u 2 , u 3 , θ) T to be the combined vector of velocity and thermal fluctuations so that (2.1)-(2.3) takes the form As it will ease computations, we also set J η = Γ −1 J Ω,Γ with η = Ω/Γ. Applying the Helmholtz projection, denoted as P, onto divergence free vector fields to the velocity component of (2.4) we obtain For periodic boundary conditions, P is defined in Fourier space as with δ ij the Kroenecker delta function and k h = (k 1 , k 2 ) T ∈ R 2 , k 3 = 2nπ, n ∈ Z. For stress-free boundary conditions P takes a similar form where one must consider sine and cosine series in the vertical direction.
With this decomposition, the x 3 -dependent terms have the following Fourier decomposition, where the ũ i and ω i terms can be obtained from the inverse Fourier transform. The Biot-Savart law for these boundary conditions can be found to be while the skew-Hermitian term PJ η P takes the form where |k| 2 := |k h | 2 + (πn) 2 and |k η | 2 := |k h | 2 + (ηπn) 2 . The matrix defined above once again has a double eigenvalue at 0 and a pair of eigenvalues ±ip η (k) = ±i|k η |/|k| with corresponding eigenvectors We can then decompose (2.5) into the following system Taking the two-dimensional curl of the equation forū h we obtain andũ andω can be related via the Biot-Savart law (2.19). We shall use this formulation to prove Theorem 2 in Sections 4.1 -4.3 below.

Existence and asymptotics in algebraically weighted spaces
Our main result in this section concerns the existence and asymptotics of solutions of the Boussinesq equations (2.4) with periodic boundaries conditions (1.14) for initial data that lie in algebraically weighted spaces. The existence result follows in a fairly standard fashion by rewriting the equation as an integral equation, coupled with estimates on the linear evolution. We note that the coupling of the vorticity formulation with temperature fluctuations in our system requires a more subtle analysis of the linear evolution; see Appendix A. The analysis of the asymptotics makes use of scaling variables which have been very useful in previous studies of the asymptotics of the Navier-Stokes equations. We work with the baroclinic-barotropic decomposition (2.13) described in Section 2.

Dynamics on barotropic invariant subspace
To motivate our results, we first consider the system on the subspaceṽ ≡ 0, which is invariant under the evolution of the system of equations (2.13). In this subspace, the equations reduce to: Finally note that the equation forū h is independent of the evolution of (ū 3 ,θ) (and also of Ω -the horizontal components of the barotropic velocity are unaffected by the rotation of the system). Furthermore, this is just the two-dimensional Navier-Stokes equation for which the long-time asymptotics are well understood [18], [19]. In particular, all solutions with integrable initial vorticity converge to an Oseen vortex, an explicit solution of the equation with Gaussian vorticity profile. We will study dynamics using a vorticity formulation. Since velocity and vorticity have different spatial decay and regularity properties, we consider the vorticityω 3 = ∂ 2ū1 − ∂ 1ū2 of the horizontal velocity components as well as the skew- T , of the vertical velocity and temperature respectively. With these quantities, one readily obtains from (3.1), the following system Here,ū h ,ū 3 , andθ are obtained fromω 3 ,ω h , andΘ respectively via the Biot Savart laws Note also that the coupling terms in the equations for w h andΘ create an oscillation between these two components of the solution. The next natural step would be to remove this oscillation by going to a rotating coordinate system, but while this simplifies the linear dynamics, it complicates the baroclinic nonlinearity in the resulting system. Hence, we will use a rotating frame to study the linear system, and then use the stationary frame in the nonlinear system. With this in mind we introduce: Rotating coordinates and vorticity formulation Note that if we ignore all terms with derivatives on the right hand side of the equations forw h andΘ in (3.2), we obtain a linear oscillator with frequency Γ. With this in mind, we introduce a rotating system of coordinates: with I 2 the two-dimensional identity matrix, so that (3.2) now takes the form To analyze the asymptotic behavior of the solutions of these equations it is convenient to introduce "scaling variables," i.e. to rescale both the dependent and independent variables in the equations as: In terms of these variables (3.2) takes the form and u h is the velocity field associated to the vorticityw 3 via the two-dimensional Biot-Savart law (3.3) (which somewhat remarkably is unchanged by the introduction of scaling variables). As noted in Section 1.5, we study these equations in the weighted Hilbert spaces L 2 (m) and L 2 2D (m). We also need the associated weighted Sobolev spaces We also define the weighted L q spaces in a similar way The operator L then has the following spectral properties when posed on L 2 2D (m). For a more detailed account of these facts see [18, §A]

Spectral information
(i) For any m ≥ 0, L has eigenvalues λ k = −k/2 for all k ∈ N with eigenfunctions (ii) Each eigenvalue has a spectral projection P n : Hence P 0 is the projection onto the 0-eigenspace and takes the form We also denote the complementary projections as well as for completeness Q −1 = Q. We also denote Next, take the weight m > 3 in the Hilbert spaces L 2 (m) (and L 2 2D (m)) for which the spectrum of L consists of a simple eigenvalue 0, with eigenfunction the Gaussian ϕ 0 and all of the rest of the spectrum in the complex half-plane with real part less than or equal to −1/2. Thus, at least for small initial data, solution components,ω 3 (τ ), of (3.9) should converge toward some multiple of the Gaussian ϕ 0 with a rate ∼ e −τ /2 . We next turn to the equations of w h and Φ to see how u h affects their evolution. Looking at the first equation in (3.9), and replacing u h with its asymptotic limit Au 0 T is the velocity profile obtained from the gaussian vorticity ϕ 0 , we obtain the equation While we are not able to explicitly compute the entire spectrum of the operator on the right hand side of this equation, a very similar operator arises in the study of the linearization of the two-dimensional Navier-Stokes equation about an Oseen vortex (see [[19], Section 4]) and using the insights gained there, one can analyze the leading eigenvalue of this operator. One has Proof. The fact that ∇ ⊥ h ϕ 0 is an eigenfunction with eigenvalue −1/2 follows by direct computation, while the fact that the remainder of the spectrum lies in the half plane with real part less than or equal to −1 is proven in Lemma A.7 Thus, the leading order asymptotics of the w h and Φ components in (3.9) are also given by derivatives of a Gaussian If we now revert to our original, unscaled variables, and reexpress the barotropic motion in terms of the velocity, rather than the vorticity, we find that solutions of (3.1), for small initial data, should behave asū We note a few somewhat surprising facts about these asymptotics: The vertical component of the barotropic velocity,ū 3 , and the barotropic termperature variationsθ are much more strongly localized than the horizontal components of the barotropic velocitȳ u h . From the Biot-Savart law, one knows that |ū 0 h (ξ)| ∼ |ξ| −1 , as |ξ| → ∞ whereasū 3 andθ decay as Gaussians.
Remark 3.3. Note that from a dynamical systems point of view, we will show that, at least locally, any of the vortex solutions defined in the barotropic subspace {ṽ = 0} are at least locally attractive for the full equations, much as if we had a center-manifold in a finite dimensional dynamical systems.

Full dynamics
We now show that these same asymptotics hold for small solutions of the full equations (2.13). We begin by taking the curl ofū andũ to eliminate the pressure terms. As mentioned previously, becauseθ has different regularity and spatial decay properties thanω, we consider the evolution of its horizontal skew-gradient,Θ = ∇ ⊥ hθ = (∂ 2θ , −∂ 1θ ) T . Thus we study a system in terms of the variablesω,Θ,ω, and θ, obtaining , and the nonlinearities are defined as Note here thatω = (∂ 2ū3 , −∂ 1ū3 , ∂ 2ū1 − ∂ 1ū2 ) T and the third component of (ω · ∇)ū is zero due to the imcompressibility condition (3.16).
Our approach for this nonlinear system is the following. We first introduce self-similar variables, and then use the properties of the approximate solution to derive a solution for the residual. We then derive an equivalent mild/integral formulation for this residual and use estimates on the linear evolution to show nonlinear existence and asymptotics via a fixed-point argument.
Scaling variables We now convert the system (3.12)-(3.15) into scaling variables. In addition to the definition of the variablesw, u,Φ andφ introduced in (3.8), we also make the change of coordinates where as before ξ = x h √ 1+t , τ = log(1 + τ ). Note also that ∇ ⊥ hφ = Φ. We then obtain from (3.12)-(3.15) the following system of equations with the nonlinear terms Note also that the incompressiblity conditions now take the form withw and u, and Φ and φ are still related via a Biot-Savart laws, (3.3) -(3.5). Also,w andũ are related by a scaled version of the Biot Savart law (2.9) taking into account the e τ /2 factor accompanying x 3 -derivatives. Furthermore, we shall need the following estimates which are derived in [44, §4.2] (see also [18]): (3.21) Ifw ∈ L 2 (D) thenũ ∈ L q (D) for all q ∈ [2, 6] and there exists a C > 0 such that Furthermore, ifw ∈ L 2 (m) 3 thenũ ∈ H 1 (m) and there exists a C > 0 such that Note, that if one desires estimates on u h or u 3 then one simply sets u = (u h , 0) T or u = (0, 0, u 3 ) T respectively. We consider the existence, and asymptotic behavior of solutions of this equation in the space 24) with the natural norm on this space. We expect, from the form of the system (2.13), that the temperature should have the same regularity as the velocity. If the vorticity is in L 2 , then we expect (as in (3.23)) that the velocity would be in H 1 . For this reason we require increased regularity in the temperature componentφ.
(Note that the due to the way in which the scaling variables are defined, the initial values forw are the same as that forω and similarly for the remaining variables.) Motivated by our heuristic discussion in Section 3.1 above, we then define where c(τ ) := cos(Γ(e τ − 1)), s(τ ) := sin(Γ(e τ − 1)), (u h ) app is the velocity profile associated with ϕ 0 from the Biot-Savart law (3.3), and the constants are determined by the initial data Direct computation then readily gives that the solutionW app = (w app , Φ app ) T solves the barotropic subsystem (3.17) -(3.18) with (w,φ) ≡ 0. One also readily finds that the subspace L 2 2D,1 (m) 2 × L 2 2D,0 (m) 0 × L 2 2D,1 (m) 2 is invariant under the evolution of this subsystem. We can thus make the following decomposition, which splits off the leading order approximate dynamics where the residuals satisfȳ Note u R and φ R can also be obtained fromw R andΦ R respectively using the two-dimensional Biot-Savart laws (3.3) and (3.4) above. We thus consider the residual on the space and the norms (w,Φ) m := w m + Φ m , (w,φ) * ,m := w m + φ H 1 (m) .
In order to simplify the notation, we also defineW R = (w R , Φ R ) T andW = (w,θ) T . Inserting this decomposition into the full system (3.17) - (3.20), and using the aformentioned facts of the approximate solution, we find that the residualsW R andW satisfy the following system

27)
W τ = (L(τ ) +Λ)W +Ñ(W R ,W ), (3.28) with the linear operators and with the nonlinear components Here subscripts in the barotropic nonlinearities denote the first two, third, and last two components of vectors in Y , and the projectionQ : Before considering the full nonlinear equation (3.27)-(3.28), we list some results on the barotropic and baroclinic components of the linear evolution. Since the linear operators are time-dependent, we must characterize the linear evolution in terms of evolutionary families of operators in the sense of Pazy [40,Ch. 5]. We leave the details of their proof to the Appendix. We define a(τ ) = 1 − e −τ , the standard basis vectors e i ∈ R 5 , and the projections e ⊥ i = 1 − e i e T i Proposition 3.5. (Barotropic evolution) Let m > 3, 0 < µ < 1/2, α ∈ N 2 with |α| ≤ 1, andW 0 ∈ Y. Then the operatorL(τ ) +Λ generates an evolutionary family of operators, S(τ, σ), with τ ≥ σ ≥ 0, on Y . If also b mW 0 ∈ L q (R 2 ) 5 , then it satisfies the following decay estimates for any q ∈ [1, 2], Proof. See Appendix A.1.
Proposition 3.6. (Baroclinic Evolution) The linear operatorL(τ ) +Λ generates an evolutionary family of operatorsS(τ, σ), for τ ≥ σ, on L 2 (m) 3   We note here that neither S norS act diagonally on vectors in R 5 and R 4 respectively, so the mild formulation cannot be broken down into components as in [44]. Using this formulation we can then prove the following result about the asymptotic behavior of small solutions of (3.20).
The solution of the residual equation for sufficiently small initial residual data then readily implies existence and asymptotics in the full solution W = W app + W R . Corollary 3.7. There exists K 0 > 0 such that for 0 < µ < 1/2, initial data W 0 ∈ Z with W 0 (·) − W app (·, 0) m + W 0 * ,m < K 0 , there exists a unique solution of (3.17) -(3.20) in C 0 ([0, ∞], Z) which satisfies the following asymptotics for any γ > 3/4. Proof. To begin, we define a Banach space where we recall that e ⊥ 3 = 1 − e 3 e T 3 with e 3 = (0, 0, 1, 0, 0) T . We let S(τ, σ) be the direct sum of the evolution operators S(τ, σ) andS(τ, σ), and F = (F ,F ) T whereF andF are the integral terms in (3.35) and (3.36) respectively. The existence of a global solution then follows by finding a fixed point of the mapping, defined on X µ,γ . Paired with the linear estimates on S andS, this follows in a standard fashion from the estimates The estimate for (3.42) follows in a similar way to (3.41). For (3.41), since the evolutionary operatorsS andS do not act diagonally we cannot consider the nonlinearities in F component-wise as in [44]. We thus use Proposition 3.5 to estimate F and Proposition 3.6 to estimateF, both with q = 3/2.
Estimates on F: We first find where C > 0 is a constant which may change from line to line. To obtain the last inequality, we go termby-term. The two-dimensional Biot-Savart law in Proposition 3.4, the embedding L 2 2D (m) ֒→ L q (R 2 ) for all q ∈ [2/(m + 1), 2], and Hölder's inequality give The estimate ∇u h L p ≤ w 3 L p , which can be derived from [18, Lem. 2.1], along with the Gagliardo-Nirenberg inequality gives For the baroclinic contributions, the three-dimensional Biot-Savart law in Proposition 3.4, the embedding L 2 (m) ֒→ L q (D) for all q ∈ [1, 2], and finally the estimate ∇ũ From this we can then conclude and in a similar manner Here we make a slight abuse of notation, letting W R Xµ,γ and W Xµ,γ denote the X µ,γ norm of (W R , 0) and (0,W ) respectively. Continuing in this way one readily obtains estimates on the gradient terms Xµ,γ . Estimates onF: ForF we use similar estimates and (3.34) of Proposition 3.6 with q = 3/2 to find Here we have used similar estimates as for the bound onF described above as well as the following and other estimates which follow in a similar way. Then, along with similar estimates for a(τ ) 1/2 ∇F m , we obtain Xµ,γ . We thus obtain the quadratic estimate on F in (3.41) and conclude the existence of a fixed point for sufficiently small W 0 ∈ Z. This implies global existence of solutions and the structure of the X µ,γ norm implies the temporal decay prescribed in (3.38).

Global existence and asymptotics 4.1 Dispersive Estimates
The proofs of the different choices of boundary conditions in Theorems 2 are very similar, so we only provide the details of the proof for the stress-free case. Hence for the remaining sections we will consider the boundary conditions (1.13). Indeed, the proofs of these results follow very closely the strategy used by Roussier-Michon and Gallay in [17] so we will mainly just highlight the differences in the proof necessitated by the additional temperature dependence, vis-a-vis the purely rotational problem considered in that reference. In contrast to the preceding section, the dispersive nature of the linearized problem is paired with large-rotation rate, |Ω| ≫ 1, to obtain global existence for initial data in which the smallness assumptions are imposed only on the quasi-geostrophic part of the initial data.
To characterize such dispersive effects, it will be important to consider the linear Rossby-type equation ∂ tṽ + ΓPJ η Pṽ = ∆ṽ, divũ = 0, (4.1) where PJ η P was defined in (2.20) above. We note that since the spatial domain D is bounded in the vertical direction, and since this equation only acts on the baroclinic part of the solution (i.e. the part with non-trivial x 3 -dependence), all eigenvalues of the linear operator will have a negative real part ∼ −4πn 2 , where n is the Fourier index in the x 3 direction. This immediately leads to the fact that for all s ≥ 0 andṽ 0 ∈ (1 − Q)H s (D) 4 with divṽ 0 = 0, a solutionṽ(t) to (4.1) satisfies In addition, recalling that PJ η P is anti-symmetric, we see that From this, we immediately conclude Remark 4.2. Due to the form of the eigenvalues of the linearized equation, solutionsv(k, t) =v(k h , n, t) with |k| = k 2 1 + k 2 2 + 4π 2 n 2 ∼ |k 1 | + |k 2 | + |n| ≥ R, decay like ∼ e −R 2 t , so in order to understand the dispersive properties of the solutions, it suffices to study the part of the solution localized in a neighborhood of zero in Fourier space.
With this in mind, let B R = {k = (k h , n) ∈ R 2 × Z | |k| ≤ R} and S be the projection onto the geostrophic eigenspace, defined as above Because the quasi-geostrophic eigenvalue −|k| 2 corresponding to a g is independent of Ω and Γ, we cannot expect any dispersive smoothing in this mode. However for the other modes, we can prove the following estimates on solutions of the linear Rossby equation (4.1) for initial data compactly supported in Fourier space perpendicular to the geostrophic mode. (4.5) Proof. In Fourier space, equation (4.1) takes the form We decomposeṽ =ṽ + +ṽ − +ṽ g , where in Fourier space where we recall that the explicit expressions for the eigenvectors a g and a ± are given in (2.12). Note here that the incompressibility condition gives ( ṽ 0 )(k), a 0 (k) C 4 = 0. Also note that the forms of the modes ṽ ± are almost identical to those of the dispersive modes in [17,Eqn. (87)]. (see also [9, §4.2] for a discussion of similar dispersive estimates, albeit in an unbounded domain.) Following the arguments of ( [17]; App.B) more-or-less line for line, one obtains Thus, we need only estimateṽ g . However this term can be estimated directly without recourse to the duality methods used in [17]. Note that Thus, using the fact that all terms have n = 0, and the support condition onṽ 0 we find Integrating this estimate with respect to t and combining it with the estimates onṽ ± immediately yields Proposition 4.3.
Note that the restriction on the support ofv in Proposition 4.3 means that the solution of (4.1) lies in any Sobolev space H s with s ≥ 0. From this, we immediately obtain the following: Corollary 4.4. With the assumptions of Proposition 4.3, for any 1 ≤ p ≤ ∞ and 2 ≤ q ≤ ∞, with 1 p + 2 q ≤ 1, the solution of (4.1) obeys:

Global Existence
To make use of these dispersive properties in the nonlinear problem we make the following decompositioñ where λ satisfies ∂ t λ + ΓPJ η Pλ = ∆λ, divλ = 0, (4.11) with initial data λ 0 = P Rṽ0 where P R is the multiplier defined by (4.12) and χ is a smooth function with χ(k) = 1 for all |k| < 1 and χ(k) = 0 for |k| > 2. The remainder r must then solve ∂ t r + ΓPJ η Pr + N 3 = ∆r, divr = 0, (4.13) with The estimates of the previous section control the evolution of λ, while we expect r will decay exponentially fast at the linear level (see Remark 4.2). Such linear estimates are then sufficient to show that solution r(t) with r 0 small will remain so. Also note, we obtain the following local existence result using a standard fixed point argument.
Proposition 4.5. (Local Existence) For any R > 0, there exists T R > 0 such that for all Ω, Γ ∈ R and initial data v 0 ∈ X sf with ||v 0 || X sf < R the equation (2.5) has a unique local solution v ∈ C 0 ([0, Note that this local existence result imposes no restriction on the size of the initial data. We now show that if the initial size of the geostrophic component of the solution is sufficiently small, the solution can be extended for all time, by deriving a bound on the solution in the X sf -norm which is uniform in time. We will prove Theorem 4. For all v 0 ∈ X sf with S(1 − Q)v 0 sufficiently small, that is the projection onto the baroclinic-geostrophic portion of u 0 is small, there exists Ω 0 such that for all Ω ∈ R with |Ω| ≥ Ω 0 the stably-stratified system (2.5) with stress-free boundary conditions has a unique global solution v ∈ C 0 ([0, ∞), X sf ) with v(0) = v 0 . Furthermore, there exists a constant C > 0 such that ||v(t)|| X sf ≤ C for all t > 0.
As noted above, the proof of this theorem follows by showing that there exists a uniform bound of the functionalΨ (4.14) Such a bound implies that the X sf -norm of the solution is uniformly bounded and hence we can extend the local existence theorem indefinitely. Using the decomposition ofṽ defined above, a bound onΨ is equivalent to a bound on From the estimates in Section 4.1, we conclude that sup t≥0 λ(t) 2 H 1 (D) ≤ λ 0 2 H 1 (D) . In addition, Poincaré's inequality implies that r(t) 2 H 1 (D) ≤ C ∇r(t) 2 L 2 (D) . Thus, our bound on the X sf -norm of the solution will follow from a bound on We control the evolution of the various terms in Ψ(t) using the following energy estimates modeled on [17]: Proposition 4.6. There exists a constant C 1 such that if v ∈ C 0 ([0, T ], X sf ) is a solution of (2.5) for some Γ, Ω ∈ R and if v 0 as above for some R > 0, then the solutions of (4.11), (4.13), and (2.24) satisfy for any t ∈ (0, T ] : Proof. Note that for stress-free boundary conditions, the equations forū h andω 3 are exactly the same as they are in the case of rotating fluids, so the proofs of (4.17) and (4.18) follow exactly as in Proposition 2.5 of [17]. For (4.19) we must make a few small changes to account for the presence of the temperature term in our equation. If we compute 1 2 d dt ∇r 2 L 2 , the dissipative term in (4.13) gives rise to − ∆r 2 L 2 , while the term ΓPJ η Pr makes no contribution due to anti-symmetry. Thus, we need only estimate the contributions of the nonlinear term: (4.20) and the Helmholtz projector has vanished due to the fact that ∇ · r = 0. These terms are broken up and estimated in turn. For example: | (∆r) · (ū · ∇)ṽdx| = | (∆r) · (ū · ∇)λdx + (∆r) · (ū · ∇)rdx| , (4.21) and each of these terms is estimated in a fashion analogous to that used in [17], leading to the bound The remaining terms are estimated in a similar way leading to (4.19).
Note that (4.18) gives us control of the ω 3 L 1 term in Ψ(t) provided we can control the evolution of u and ∆ũ and these are in turn controlled by the evolution of r and λ. The evolution of λ is controlled by the estimates of the previous section, and thus, we turn our attention to the "reduced" functional Remark 4.7. Note that we expect the two terms in Φ to have different properties -ω 3 may be large, but is not expected to grow much, while we can make ∇r arbitrarily small (at least initially) by choosing R (in (4.12)) sufficiently large.
Differentiating with respect to t and using the estimates of Proposition 4.6, we obtain We bound the term while we boundū with the aid of the Biot-Savart law (see [18,App. B]) Recalling Remark 4.7 about the expected relative sizes ofω 3 and r, we see that (4.25) implies Lemma 4.8. There exist constants C 2 , C 3 , C 4 such that the following holds. Let v ∈ C 0 ([0, T ], X sf ) be a solution of (2.5) which is decomposed v =v + λ + r as above for some R > 0. Assume as well that there exist K ≥ 1 and ǫ ∈ (0, 1] such that the corresponding components satisfy for all t ∈ (0, T ]. Here ). (4.29) We now prove that with the aid of Gronwall's inequality, the X sf -norm of the solution of (2.5) remains uniformly bounded for all time, and hence, the local existence theorem can be extended without limit. More precisely, we show Proposition 4.9. For any initial conditions, v 0 ∈ X sf of (2.5), there exists K, Ω 0 > 0 and ǫ ∈ (0, 1), such that if |Ω| > Ω 0 , and if the projection onto the geostrophic mode, Sṽ 0 , is sufficiently small, then for any T > 0, sup 0≤t≤T ∇r(t) L 2 ≤ ǫ . Proof. Choose K and R large enough that and ∇r| t=0 L 2 ≤ 1 16 ǫ . (4.33) By the local existence theorem, there exists T * * > 0 such that (4.30) and (4.31) hold for 0 ≤ t ≤ T * * . Let T * be the supremum over the set of values of T for which these estimates hold. We claim that T * = ∞. Suppose instead that T * < ∞. Note that from the estimates of Section 4.1, if we choose Ω 0 sufficiently large, and geostrophic projection sufficiently small, we can insure that for any T > 0. Furthermore, using the estimates on λ from the previous section, plus (4.31) and (4.18), we can insure that both for all 0 ≤ t ≤ T * . If one then applies Gronwall's inequality to (4.28), one finds that for all 0 ≤ t ≤ T * . If we then apply a similar argument with Gronwall's inequality to ∇r(t) L 2 , we also find that for Ω 0 sufficiently large, and the geostrophic projection sufficiently small, we have for all 0 ≤ t ≤ T * . However, these two estimates imply that we could extend the time for which (4.30) and (4.31) hold beyond T * , contradicting its definition. Hence, T * = ∞ as desired. Proof. This once again follows in a similar way as in [17]. We already have by Poincare's inequality that
It follows from this proposition that w 3 (τ ) − AG L 1 (R 2 ) → 0 as τ → ∞. Translating this result into unscaled variables one can obtain leading order asymptotics for the solutionω(x h , t) (and thusū h (x h , t)) to complete the proof of Theorem 2.
Since we must study different components of the evolution we define the following projection operators where e i are the standard basis vectors in R 5 . Also recall that b = (1 + |ξ| 2 ) 1/2 . Next, recall that with the standard L 2 2D (m) norm. We first state asymptotics for the strongly continuous semi-group e τ L generated by the operator L on L 2 2D (m).
One can then readily conclude the following temporal estimates on S 0 .
Transforming the nonlinear equation (3.7) on the barotropic subspace {ṽ = 0} into scaling variables, we obtain an equation in the rotating frame for the (w, Φ).
We use the ansatz with the approximate solution 2D,0 (m). Note that w app and Φ app just correspond to the explicit vortex solution discussed in Section 3.2. Alternatively, it can be obtained from the stationary frame approximate solution by re-writing (3.25) in the rotating frame. Using the fact that the approximate solution solves (A.8) we find that the perturbation of the linear evolution satisfies We can then split Λ into three parts, where Λ 0 is exactly the same as in the 2-D Navier Stokes case discussed in [19], and Λ 2 consists entirely of terms decaying like e −τ /2 . For the rest of the section we shall omit the R subscript.
A.1.1 Evolution generated by L + Λ 0 + Λ 1 We then use this decomposition to obtain the following proposition which gives temporal decay estimates on the linear evolution generated by L + Λ 0 + Λ 1 defined on Y . Recall that a(τ ) = 1 − e −τ .
Proposition A.3. The linear operator L + Λ 0 + Λ 1 defines a strongly continuous semi-group T 1 (τ ) on Y which, for f ∈ Y , |α| ≤ 1, q ∈ [1, 2], m > 3, satisfies the following estimates, We shall split the proof of this result into a series of results. We begin by using Duhamel's formula to write solutions of the evolution equation Estimating the integral terms as we did the corresponding expressions in Section 3 one readily shows that T 1 (τ ) defines a strongly continuous semigroup such that for any T > 0, there exists C T > 0 such that for any 0 ≤ τ ≤ T , and multi-index α with |α| = 1, The more refined estimates of the evolution given in Proposition A.3 will follow from a more detailed analysis of the semigroup. We begin the estimates of T 1 (τ ) by noting that because of the form of Λ 0 and Λ 1 , it splits into three independent parts which govern the evolution of W R,h = (W R,1 , W R,2 ) T , W R,3 and W R,p = (W R,4 , W R,5 ) T . We refer to these three parts of T 1 as T 1,h , T 1,3 and T 1,p and will estimate each of them separately. For later use, we will refer to the pieces of L + (Λ 0 + Λ 1 ) which generate each of these semigroups as (L + (Λ 0 + Λ 1 )) h , and so forth.
As a first step, note that T 1,3 (τ ) is exactly the semigroup studied in Section 4 of [19], as it is the linearization of the two-dimensional Navier-Stokes equation. Thus, estimates (A.11) follow from the results of Proposition 4.13 of that work.
To prove the estimates in (A.10) and (A.12) we use methods similar to those in Section 4 of [19]. We will give the details for (A.10), as the case of (A.12) follows in the exact same way. The proof of (A.10) begins by noting that these estimates are exactly those we would obtain for the horizontal components of the diagonal semigroup S 0 (τ ) (which we will denote S 0,h .) Thus, we proceed in two steps: (i) We first show that the difference between S 0,h (τ ) and T 1,h (τ ) is a compact operator for any τ > 0.
As a consequence, the essential spectral radius of T 1,h (τ ) is the same as that of S 0,h (τ ).
(ii) The first point implies that the only way that the decay rates of the two semigroups can differ is if T 1,h (τ ) has an isolated eigenvalue lying outside the disc which gives the essential spectral radius. Thus, the second step in our analysis is to determine the location of isolated eigenvalues of (L + Λ 0 + Λ 1 ) h .
Thus we shall prove the following two results Proposition A.4. Assume the hypotheses of Proposition A.3. Then for any τ > 0, the linear operator Proposition A.5. Assume that m > 3. Then any eigenvalue, λ, of L + Λ 0 + Λ 1 h on Y h,1 satisfies Combining these two propositions, the estimates (A.10) follow exactly as the proof of Proposition 4.13 in [19]. The proof of Proposition A.4 follows immediately from Rellich's criterion and the following lemma.
Lemma A.6. Assume the hypotheses of Proposition A.3. Then for any T > 0, there exists a constant C T > 0 such that for 0 ≤ τ ≤ T , for any τ > 0, the linear operator Proof. From Duhamel's formula we can write The estimates of the lemma now follow by estimating each of the integral terms. We first have Here, the first inequality used our estimates on the semigroup S 0 (τ ), the second inequality used the fact that u 0 h f m+1 ≤ C f m due to the decay of u 0 h (ξ) as |ξ| → ∞, the third inequality used (A.15), and the last one the fact that the singularity in a(τ ) −1/2 is integrable.
Estimating the · m+1 -norm of the first integral term in (A.18) is similar, but even easier since derivatives of u 0 h (ξ) decay more rapidly as |ξ| → ∞ than u 0 h (ξ) itself, and we leave this estimate to the reader. Now we turn to the second estimate in the lemma, namely the estimate of the derivatives of K h . Differentiating the expression in (A.18) we must estimate (A.22) Once again, because of the explicit formulas for the derivatives of u 0 h (ξ), the estimates of the first term are easier than the second, so we focus on estimating the second. Using the estimates on the semi-group S 0,h , plus the estimates in (A.15), we have Now, using the fact that the estimate, and the lemma follows.
We now turn to the proof of Proposition A.5. We begin with: Lemma A.7. λ is an eigenvalue of the operator posed on Y h,1 , if and only if λ is an eigenvalue of the operator 8. Note that our estimates of the Biot-Savart kernel imply that the velocity field u 3 ∈ L 2 2D,0 (m ′ ).
Proof. The proof of the lemma follows by noting that if W and λ satisfy then we set U = u 3,BS (W ) where u 3,BS is the Biot-Savart mapping (3.4) defined above so that (∂ 2 U, −∂ 1 U ) = W , and U ∈ L 2 2D,0 (m ′ ). Indeed our choice of m, m ′ , Proposition B.1 of [18], and the generalized Holder inequality gives where 1/2 = 1/p + 1/q ′ with q ′ ∈ (2, ∞), chosen so that m ′ − m ′′ > 2/p (and hence b m ′ −m ′′ ∈ L p ) and m ′′ = m − 2/q ′ , (so [18, Prop. B.1] holds). Now, we compute the derivatives of h := L 3 U to find Applying the Biot-Savart law we then find by linearity that h = λU and thus Proposition A.5 now follows from Lemma A.9. Let m ′ be as in Lemma A.7 . Then any eigenvalue, λ, of the operator L 3 whose eigenfunction lies in the space L 2 2D,0 (m ′ ), satisfies Proof. To prove this lemma, we first note that any eigenfunction of L 3 decays as a Gaussian as |ξ| → ∞. This is proven in a fashion almost identical to the analogous result in [19]. More precisely, if we write the eigenfunction ψ in polar coordinates we find that ψ = Ψ(r)e iℓθ where Ψ satisfies the ODE where η(r) = 1 2πr 2 (1 − e −r 2 /4 ). This is almost identical to equation (70) of [19] (in fact, significantly simpler than that case, because the nonlocal term denoted Ω(r) in that reference is absent) and using the same sort of ODE estimates applied there one finds that there exists γ > 0 such that |Ψ(r)| ≤ C(1 + r 2 ) γ e −r 2 /4 . (A.28) To localize the eigenvalues of L 3 we conjugate the two pieces of the operator with the Gaussian ϕ −1/2 0 . Alternatively, we study the operator on the weighted Hilbert space From Lemma 4.7 of [19], we know that L − 1/2 is a selfadjoint operator on this space with pure point spectrum −1, −3/2, −2, . . . . A direct calculation shows that the operator f → u 0 h · ∇f is anti-symmetric in this space. Explicitly, we have where we have used the fact that div (ϕ −1 0 u 0 h ) = 0. But now we use the fact that if an anti-symmetric operator is added to a self-adjoint operator with pure point spectrum, the real part of the spectrum of the resulting operator is less than the largest eigenvalue of the self-adjoint part -in this case −1. (See the proof of Proposition 4.1). in [19].

A.1.2 Evolution generated by L + Λ
We now can find and estimate solutions of the full linear equation (A.9), which being non-autonomous, generates an evolutionary family of operators S(τ, σ), τ ≥ σ ≥ 0 (in the sense of Pazy [40,Ch. 5]) instead of a semi-group.
Proposition A.10. The operator L+Λ defines an evolutionary family of operators S(τ, σ) which satisfy the following temporal estimates for q ∈ [1, 2], |α| ≤ 1, m > 3, µ Proof. First, the existence and uniqueness of S(τ, σ) can be obtained using a fixed point argument with the mapping where we have set W (σ) = f . To obtain the temporal estimates we use a Gronwall type argument which uses the fact that the coefficients of Λ 2 decay with rate e −τ /2 . Due to the uniqueness of the above fixed point, S(τ, σ) can also be represented by the following fixed point formula We then wish to show that is bounded for all τ ≥ σ. Here we weaken the temporal decay rate on the W h and W p components in order to be able to close the Gronwall argument for Ψ. We use Proposition A.3 to obtain Note in the second line we used the estimate ∇u R,h L q ≤ C W 3 L q which holds for all q ∈ [1, ∞] (see [44]) as well as the Biot-Savart estimate (3.21) to obtain u h L 6 ≤ C W 3 L 3/2 , while in the third line we used the fact that L 2 2D (m) ֒→ L q (R 2 ) for all q ∈ [1, 2] to conclude W 3 L 3/2 ≤ C W 3 m . Also, C B > 0 is a constant dependent on the amplitudes B i for i = 1, 2. We can then conclude Using the fact that the third componentW 3 is unaffected by Λ 2 , we similarly obtain as well as the gradient estimates Combining these together we obtain the inequality Then, applying Gronwall's inequality, we obtain that which since µ < 1/2 implies that Ψ(τ, σ) is bounded for all τ ≥ σ. We thus obtain the temporal estimates for q = 2 on each component of W (τ ) by unraveling the definition of Ψ and taking note of the invariance properties of L + Λ. A similar argument then gives the desired L q estimates.

A.1.3 Evolution generated byL +Λ
We now translate the estimates of Proposition A.10 back into the stationary frame in order to conclude estimates on the linear flow forL(τ ) +Λ.
Proposition A.11. Let m > 3, α ∈ N 2 with |α| ≤ 1, andW 0 ∈ Y, then the operatorL(τ ) +Λ generates a evolutionary family of operators,S(τ, σ) for τ ≥ σ ≥ 0 which satisfies the following decay estimates Proof. We setW (ξ, σ) = W 0 (ξ) as the initial condition of the linear evolutionW (ξ, τ ) = S(τ, σ)W (ξ, σ) which solvesW τ = (L(τ ) +Λ)W . In order to use Proposition A.10, we translate to the corresponding evolution in the rotating frame. We first set w 3 (ξ, σ) = W 0,3 (ξ), Since it is untouched by the rotating frame transformation. For the other components, by translating back to unscaled variables, we set with J 2 as in (A.1). We then let W (τ ) := S(τ, σ)W (σ), be the solution of (A.9) with initial data W (σ) as defined above. Note w 3 (τ ) := w 3 (τ ) as the third vorticity component is not affected by the rotating frame. Next, defining R(τ ) = e (e τ −1)ΓJ2 , we have where the subscripts at the end of the first and third lines denote the component of the vector and we recall that c(τ ) := cos(Γ(e τ − 1)), s(τ ) := sin(Γ(e τ − 1)). Next we use the estimates in the rotating frame to estimate different components of W . Using the estimates of Proposition A.10 we find since the rotational transformation does not affect the norm of W 0 . The estimates for the gradients, as well as L q data, follow in a similar way.

A.2 Baroclinic evolution
A.2.1 Unperturbed evolution ofL(τ ) First we estimate the unscaled velocityṽ = (ũ,θ) T under the linear evolution,ṽ t = Mṽ := ∆ṽ+PJ Ω,Γ Pṽ, in algebraically weighted spaces. One readily finds that this equation generates a C 0 -semigroup which we denote by e Mt . The approach used in [44,Prop. 4.4.2] can be used to readily find Proposition A.12. Let m ≥ 0, α ∈ N 3 , q ∈ [1,2], and f ∈ (1 − Q)L 2 (m) 4 . If also b m f ∈ L q (D) then ∂ α e M(t−s) f m ≤ Ce −4π 2 (t−s) (t − s) 1/q−1/2+|α h |/2 a(t − s) (1/q−1/2)/2+α3|/2 b m f L q (D) . (A.39) To transfer these results to the vorticity formulation, we need to show that the unscaled Biot-Savart law, which mapsω to ∂ jũ for j = 1, 2, 3 is a bounded operator on the b m -weighted L q spaces. We view this operator as a matrix Fourier multiplier and aim to apply the theory of weighted L q multipliers developed over the past several decades; see [20] for a review. In particular, a weight function ρ ≥ 0 with ρ ∈ L 1 loc (R n ) is said to be Muckenhaupt class A q (R n ) for some q ∈ (1, ∞) if for all bounded cubes Q with faces parallel to coordinate axes, where |Q| denotes the Lebesque measure of Q. For our weight b m considered in R 2 we first have Lemma A. 13. The weight b m := (1 + |x h | 2 ) m/2 is of class A q (R 2 ) for all q > 1 and 0 ≤ m ≤ 2(q − 1).
Following [46] we can then define a Muckenhaupt class of weights, A q (D), for our spatial domain D = R 2 × T in the same way as above; see [46, pg. 338]. From this we then can conclude the following lemma: Lemma A.14. For all q > 1 and 0 ≤ m < 2(q − 1), the weight b m := (1 + |x h | 2 ) m/2 is of class A q (D).
Proof. Let Q = I 1 × I 2 × I 3 , for some bounded intervals I 1 , I 2 ⊂ R, I 3 ⊂ T. Then we use the fact that b m does not depend on x 3 and Lemma A.13 to obtain for some positive constant C > 0.
For this class of weights, the results of [46] can be used to show that ∂ jũ , for j = 1, 2, 3, are bounded operators in the b m weighted L q space. Namely we obtain the following: Lemma A.15. Letω ∈ (1 − Q)L 2 (m) satisfy b mw ∈ L q (D) for some m ≥ 0 and q ∈ [3/2, ∞). Also let u be determined by the Biot-Savart law given in (2.9) above. Then there exists a C > 0 such that for for j = 1, 2, 3, Proof. First we find observe that ∂ jũ , when defined on R 3 , is obtained fromω via a Fourier matrix multiplier B j (k) with components of the form kikj |k| 2 , k ∈ R 3 . This implies that B j (k) is bounded for k ∈ R 3 , is smooth away from the origin, and satisfies the Mikhlin condition |k| β |∂ β k B j (k)| ≤ C for some C > 0. Thus B j (k) defines a bounded operator on L q (R 3 ).
Next, letting β ∈ N 2 denote a multi-index and m 1 ∈ N, it can readily be found that there exists constants C, C ′ > 0 such that Denoting F as the Fourier transform on D, m 1 = ⌊m⌋ the greatest integer below m we then have Here we have used equivalencies (A.43) in the second and last lines and the L q (m − m 1 ) boundedness of the multipliers ∂ β k h B j , for all m − m 1 < 1 < 2(q − 1), in the fourth line.