Theoretical and Computational Analysis of the Thermal Quasi-Geostrophic Model

This work involves theoretical and numerical analysis of the thermal quasi-geostrophic (TQG) model of submesoscale geophysical fluid dynamics (GFD). Physically, the TQG model involves thermal geostrophic balance, in which the Rossby number, the Froude number and the stratification parameter are all of the same asymptotic order. The main analytical contribution of this paper is to construct local-in-time unique strong solutions for the TQG model. For this, we show that solutions of its regularised version \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α-TQG converge to solutions of TQG as its smoothing parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \rightarrow 0$$\end{document}α→0 and we obtain blow-up criteria for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α-TQG model. The main contribution of the computational analysis is to verify the rate of convergence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α-TQG solutions to TQG solutions as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \rightarrow 0$$\end{document}α→0, for example, simulations in appropriate GFD regimes.

1. Introduction 1.1.Purpose.The thermal quasi-geostrophic (TQG) equations comprise a mesoscale model of ocean dynamics in a solution regime near thermal geostrophic balance.In thermal geostrophic balance, three forces -the Coriolis, hydrostatic pressure gradient and buoyancy-gradient forces -sum to zero.Numerical simulations of the TQG equations show the onset of instability at high wavenumbers which creates small coherent structures which resemble submesoscale (1-20 km) features observed in satellite ocean color images as seen in figure 1.
Figure 1.Comparison of a computational simulation of TQG solutions with a satellite observation of ocean colour on a section of the Lofoten Vortex, courtesy of https: //ovl.oceandatalab.com/ which illustrates the configurations of submesoscale currents obtained from ESA Sentinel-3 OLCI instrument observations of chlorophyll on the surface of the Norwegian Sea in the Lofoten Basin, near the Faroe Islands.
On the left panel of Figure 1 one sees submesoscale features which are prominently displayed in computational simulations of TQG equations for sea-surface height (SSH).The right panel of Figure 1 shows the surface of the Lofoten Basin off the coast of Norway near the Faroe Islands.In crossing the Lofoten Basin, warm saline Atlantic waters create buoyancy fronts as they meet the cold currents of the Arctic Ocean. Figure 1 displays several features of submesoscale currents surveyed in [25].High resolution (4km) computational simulations of the Lofoten Vortex have recently discovered that its time-mean circulation is primarily barotropic, [32], thereby making the flow in the Lofoten Basin a reasonable candidate for investigation using vertically averaged dynamics such as the TQG dynamical system.Both images show a plethora of multiscale features involving shear interactions of vortices, fronts, plumes, spirals, jets and Kelvin-Helmholtz roll-ups.The submesoscale features persist and interact strongly with each other in Kelvin-Helmholtz roll-up dynamics, instead of simply cascading energy to higher wavenumbers.This observation means that the instabilities which create these submesoscale features quickly regain stability without cascading them to ever smaller scales.The present work aims to understand these features of the TQG solution dynamics, both analytically and numerically.
The contributions of this paper.The main analytical contribution of this paper is the construction of a unique local strong solution of the TQG model.In particular, we show that there is a unique maximal solution of the TQG equation defined on a (possibly infinite) time interval [0, T max ).The solution will be shown to exist in a suitable Sobolev space.More precisely, provided that the initial data resides in a chosen Sobolev space, the solution at time t > 0 will remain in this space as long as t < T max .Should t < T max be finite, then the solution will blow up in the chosen Sobolev norm.
As a second analytical contribution, we show that the TQG model depends continuously on the initial condition.This dependence only holds in a slightly weaker norm than the one corresponding to the space where the solution resides.This is a useful property from a numerical perspective.It implies that initial small errors when simulating the TQG model will stay small at any subsequent time.
The third analytical contribution is to construct a regularized version of the TQG model, termed the α-TQG model.This model is constructed in a similar manner as the α model for the Euler and Navier-Stokes equations, see [11,12,23].The α-TQG equations have a unique maximal solution which is also only continuous with respect to the initial conditions in a larger space with weaker norm.In addition, we show that the α-TQG solution converges to the TQG solution as α → 0 in a norm that depends on two physical parameters, the vorticity and the gradient of buoyancy.We also identify the rate of convergence as a function of the α-parameter on a time interval where both the TQG solution as well as the α-TQG solution are shown to exist for any α > 0.
The blow-up phenomenon is important, particularly if it is observed in numerical simulations.Therefore, blow-up criteria (in other words, criteria required for the solution to blow up) are important.In this paper we state three characterizations of the blow-up time.The fourth analytical contribution of this paper is to show that blow-up occurs in the α-TQG model, if either: (1) the L ∞ (T 2 )-norm of the buoyancy gradient ∇b blows up at T max < ∞; (2) the L ∞ (T 2 )-norm of the velocity gradient ∇u blows up at T max < ∞ or that; (3) the W 1,2 (T 2 )-norm of the buoyancy gradient ∇b blows up at T max < ∞.
The contribution of this paper from a numerical perspective is to describe our spatial and temporal discretisation methods for approximating TQG and α-TQG solutions, and to analyse aspects of numerical conservation properties with respect to the theoretical conserved quantities of the α-TQG system.Example simulation results are included and are used to verify numerically the theoretical convergence results.Additionally, we provide linear stability analysis results for the α-TQG system.Given that we have convergence of α-TQG solutions to TQG, the linear stability results can be viewed as generalisations of those shown in [17] for the TQG system.1.2.Brief history of the TQG model.The history of the TQG model goes back about half a century, perhaps first elucidated by O'Brien and Reid [26] as sketched, in [6].Briefly put, the TQG model generalises the classical QG equations by introducing horizontal gradients of buoyancy which alter the geostrophic balance to include the inhomogeneous thermal effects which influence buoyancy.Indeed, O'Brien and Reid [26] write that their work was inspired by observations that the passage of hurricanes could draw enough heat from the ocean to significantly lower the sea-surface temperature in the Gulf of Mexico.Since O'Brien and Reid [26] introduced their two-layer model, further developments of it have been applied to a variety of ocean processes, particularly to equatorial dynamics.For more details of the theoretical model developments, see Ripa [27,29,30] and for developments of applications in oceanography see [1,3,24,31], as well as other citations in [6].In particular, Ripa refers to the TQ models as inhomogeneous-layer (IL) models and his papers explain rational derivations of theories with increasing vertical structure IL1, IL2, etc.The TQG model analysed here and derived systematically from asymptotic expansions in small dimensionless parameters of the Hamilton's principle for the rotating, stratified Euler equations in [17] is equivalent to Ripa's model IL0QG [28] recently analyzed in [5].1.3.A sketch of the derivation of the TQG model.We have explained that certain thermal effects in the mesoscale ocean have historically been modelled by the thermal quasi-geostrophic (TQG) equations.TQG is characterised by several dimensionless numbers arising from the dimensional parameters of planetary rotation, gravity and buoyancy.These are the familiar Rossby number, Froude number and stratification parameter.The Rossby number is the ratio of a typical horizontal velocity divided by the product of the rotation frequency and a typical horizontal length scale.The Froude number is the ratio of a typical horizontal velocity divided by the velocity of the fastest propagating gravity wave, which in turn is given by the square root of the gravity times the typical vertical length scale.The final dimensionless number is the stratification parameter, which specifies the typical size of the buoyancy stratification.
The regime in which the thermal quasi-geostrophic equations are derived is characterised by the thermal geostrophic balance.This balance arises because the Rossby number, the Froude number and the stratification parameter all have a similar amplitude.Preserving this three-fold balance requires simultaneously adapting the Froude number and the stratification parameter to match any change in Rossby number, for example, so that the dimensionless parameters will still have the same size.We consider the non-dissipative case, because of the large scales of mesoscale ocean dynamics.As mentioned earlier, the mesoscale dynamics has high wavenumber instabilities, which in principle can generate submesoscale effects.At smaller scales, viscous dissipation and thermal diffusivity will come into play as well.However, in what follows, dissipative effects will be neglected.The derivation of the TQG model involves a series of coordinated approximations to the rotating, stratified Euler model, as is illustrated in the diagram below.A discussion of the right-most three columns tree in Figure 2 can be found in [16].The derivation of the thermal quasi-geostrophic equations treated here from the thermal rotating shallow water equations on the middle left of the figure can be found in [17].The derivation of the thermal quasi-geostrophic equations following the solid red arrows in Figure 2 starts with the Lagrangian for the rotating, stratified Euler equations at the top of the figure.After identifying the small dimensionless parameters in the Euler Lagrangian, the Lagrangians for the successive approximate models can be derived by inserting asymptotic expansions into the Euler Lagrangian.In the ocean, the buoyancy stratification is typically small.Hence, it makes sense to apply the Boussinesq approximation in the Lagrangian for the rotating, stratified Euler equations.This approximation yields the Lagrangian for the Euler-Boussinesq equations.At this point, one has a choice of two routes.The first route begins by making the hydrostatic approximation, which leads to the Lagrangian for the primitive equations.Then, upon vertically integrating the Lagrangian for the primitive equations, one finds the Lagrangian for the thermal rotating shallow water equations.The alternative route first vertically integrates the Lagrangian for the Euler-Boussinesq equations to obtain the Lagrangian for the thermal and rotating version of the Green-Naghdi equations.By subsequently making the hydrostatic approximation in the Lagrangian for the thermal rotating Green-Naghdi equations, the alternative route arrives at the Lagrangian for the thermal rotating shallow water equations.
The Lagrangian for the thermal rotating shallow water (TRSW) equations is the starting point in [17] for the derivation of the thermal quasi-geostrophic (TQG) equations.The typical values of the dimensionless numbers for mesoscale ocean problems lead to thermal geostrophic balance.This balance implies an algebraic expression for the balanced velocity field in terms of the horizontal gradients of the free surface elevation and the buoyancy.Upon expanding the TRSW Lagrangian around this balance and truncating, one obtains the Lagrangian for the thermal Lagrangian 1 (L1) model.The Lagrangian for the thermal L1 model is not hyperregular, though.Hence, the Legendre transformation is for thermal L1 is not invertible.Hence, no Hamiltonian description of this model is available via the Legendre transformation.Nonetheless, an application of the Euler-Poincaré theorem [18] to the thermal L1 Lagrangian yields the corresponding equations of motion.By identifying the leading order terms in the thermal L1 equations and truncating the asymptotic expansion, one finally obtains the TQG equations.However, the truncations of the asymptotic expansions in the thermal L1 equations also prevent the resulting TQG equations from possessing a Hamilton's principle derivation.This feature is unlike the other models in figure 2, which all arise from approximated Lagrangians in their corresponding action integrals for Hamilton's principle.However, as it turns out, the thermal quasi-geostrophic equations do possess a Hamiltonian formulation in terms of a non-standard Lie-Poisson bracket.The Lie-Poisson bracket for the quasi-geostrophic equations can be obtained via a linear change of variables from the usual semi-direct product Lie-Poisson bracket for fluids.The resulting Hamiltonian formulation for the thermal quasi-geostrophic equations is important for the future application of stochastic advection by Lie transport (SALT), introduced in [15], which requires either a Lagrangian, or a Hamiltonian interpretation of the equations of motion.
By applying the method of Stochastic Advection by Lie Transport (SALT) to the Hamiltonian formulation of the TQG equations, one obtains a stochastic version of them which preserves the infinite family of integral conserved quantities.More details on the stochastic version can be found in the conclusion section as well as in [17].A discussion of the remaining models in Figure 2 can be found in [16].
The TQG model on the two dimensional flat torus T 2 can be formulated in two equivalent ways.The first formulation is the advective formulation, in which the equations of motion are given by Here b is the vertically averaged buoyancy, u is the thermal geostrophically balanced velocity field, ω is the potential vorticity, ψ is the streamfunction, h 1 is the spatial variation around a constant bathymetry profile and f 1 is the spatial variation around a constant background rotation rate.The velocity u and the streamfunction ψ are related by the equation where ∇ ⊥ = (−∂ y , ∂ x ).The vector field u h1 is defined by Alternatively, one can formulate the thermal quasi-geostrophic (TQG) equations in vorticity-streamfunction form.This formulation is given by a y is the Jacobian of two smooth functions a and b defined on the (x, y) plane.The equivalence of the two formulations (1.1)-(1.2) and (1.6)-(1.7)follows from the Jacobian operator relation J(ψ, a) = u • ∇a, in which ψ is the streamfunction associated to the velocity vector field u.The scalar functions f 1 and h 1 relate to the usual Coriolis parameter and bathymetry profile in the following way where Ro = U (f 0 L) −1 is the Rossby number, expressed in terms of the typical horizontal velocity U , typical rotation frequency f 0 and typical horizontal length scale L. This means that the bathymetry and Coriolis parameter become constant as the Rossby number tends to zero.Equations (1.9) are necessary to derive the thermal quasi-geostrophic equations from the thermal L1 equations, as shown in [17].The expansion (1.9) contains the β-plane approximation provided that the boundary conditions are appropriate.Namely, on the β-plane, one requires that βf −1 0 = O(Ro) and f 1 (x) = y.An additional relation can be helpful when using the thermal quasi-geostrophic equations as a model for mesoscale ocean dynamics.In particular, the streamfunction is related to the free surface elevation ζ and buoyancy b via the definition ψ := ζ + 1 2 b.This definition of the streamfunction is useful to relate to observational data.The free surface elevation is a quantity that can be measured with satellite altimetry.These measurements can then be used for data assimilation and model calibration.However, the definition of the streamfunction in terms of the velocity field is not necessary to formulate the model as a closed set of equations, since the system (1.6)-(1.8) is already a closed set of equations.In the expansion (1.9), we will henceforth drop the subscript 1 on h 1 (x) and f 1 (x) for notational convenience.
1.4.TQG versus Rayleigh-Bénard convection with viscosity and thermal diffusivity.The TQG equations in (1.6)-(1.8)can be compared to the equations for Rayleigh-Bénard convection in a vertical plane, as remarked in [17].Recent results of [9] show that bounds exist for the enstrophy and temperature gradient of solutions lying in the attractor for the planar Rayleigh-Bénard convection problem which are algebraic in the viscosity and thermal diffusivity.This result is a significant improvement over previously established estimates.This result also provides a motivation to study the TQG equations with viscous dissipation and thermal diffusivity.The equations for non-dissipative Rayleigh-Bénard convection in the vertical plane are given by (1.10) where ω denotes the vorticity, ψ is the stream function, T is the temperature, g is gravity and α is the thermal expansion coefficient.The Rayleigh-Bénard equations (1.10) have a long and illustrious history in mathematical analysis.The resemblance of the Rayleigh-Bénard equations to the TQG equations suggests that perhaps the TQG equations will also provide a fruitful challenge to mathematical analysis.
In TQG, the buoyancy plays the role of the temperature in the Rayleigh-Bénard equations.By rearranging the TQG equations (1.8) such that the buoyancy terms all appear on the right hand side, one has (1.11) There are similarities and also several differences between the Rayleigh-Bénard equations (1.10) and the TQG equations (1.11).For example, the equation that relates the vorticity to the stream function is a Helmholtz equation for TQG and a Poisson equation for Rayleigh-Bénard.In addition, in Rayleigh-Bénard convection the forcing term on the right hand side of the vorticity equation depends only on the derivative of the temperature in the vertical direction.In TQG the forcing terms depends on the derivatives of the buoyancy in both directions.Thus, TQG arises as an interesting and challenging extension of the celebrated Rayleigh-Bénard problem.
1.5.Plan for the rest of the paper.We now give the plan for the rest of the paper.We collect preliminary tools in Section 2. This includes notations and analytical properties of function spaces used throughout this paper.We also collect useful estimates that will be used at various stages and give precise definitions of the concept of solutions used in our analysis.We finally end Section 2 with statements of the main results.In particular, we state that both the TQG and α-TQG equations admit unique strong solutions for a finite period of time and these solutions are stable in a larger space with weaker norm.Furthermore, a maximum time for these solutions exists.
Since the proof of local well-posedness is the same for the TQG and α-TQG, we will avoid duplication by devoting Section 3 to the construction of solutions for the less regular TQG.The construction relies heavily on the standard energy method.Since we are constructing strong solutions (rather than weak ones), we differentiate the equations in space and then we test the resulting equations with the required differential of the solution to obtain the required bounds.We then end Section 3 by showing that the unique solution constructed has a maximum time of existence and hence, is a maximal solution.
In Section 4 we show that any family of maximal solutions of the α-TQG models converges strongly with α → 0 to the unique maximal solution of the TQG on a common existence time, provided they share the same data.
Next, since our solutions are local in nature, we establish in Section 5, conditions under which this solution may blow up in the sense of Beale-Kato-Majda [2].In particular, we show that in order to control the solution of α-TQG, the essential supremum in space of both the buoyancy gradient and the velocity gradient should be integrable over the anticipated time interval.Once either of these gradients blow up, the solution ceases to exist.Alternatively, in order to control the solution, it suffices to control the H 1 -Sobolev norm of the buoyancy gradient.
Section 6 is devoted to numerical methods and simulation results.We begin by describing the finite element method we use for the spatial derivatives, and aspects of its numerical conservation properties with respect to theoretical results.We then describe the finite difference discretisation method used for the time derivative.Next, we discuss α-TQG linear thermal Rossby wave stability analysis, which can be seen as generalising the results shown in [17] for the TQG system.Then in the last part of the subsection, we discuss our numerical simulation setup and its results.In particular, we numerically verify the theoretical convergence rate for α-TQG → TQG derived in Section 4.

Preliminaries and main results
In this section, we fix the notation, collect some preliminary material on function spaces and present the main analytical results.
Remark 2.1.Although we will be working on the 2-dimensional torus, with minimal effort, the same analysis will work on the whole plane R 2 subject to a far-field condition.The case of a bounded domain with boundary conditions is however outside the scope of the analytical aspect of this work.See Section 6 for numerical works in this regard.
2.1.Notations.Our independent variables consists of spatial points x := x = (x, y) ∈ T 2 on the 2-torus T 2 and a time variable t ∈ [0, T ] where T > 0. For functions F and G, we write F G if there exists a generic constant c > 0 such that F ≤ c G. We also write F p G if the constant c(p) > 0 depends on a variable p.The symbol | • | may be used in four different context.For a scalar function f ∈ R, |f | denotes the absolute value of f .For a vector f ∈ R 2 , |f | denotes the Euclidean norm of f .For a square matrix we denote by W k,p (T 2 ), the Sobolev space of Lebesgue measurable functions whose weak derivatives up to order k belongs to L p (T 2 ).Its associated norm is where β is a 2-tuple multi-index of nonnegative integers of length |β| ≤ k.The Sobolev space W k,p (T 2 ) is a Banach space.Moreover, W k,2 (T 2 ) is a Hilbert space when endowed with the inner product where • , denotes the standard L 2 -inner product.In general, for s ∈ R, we will define the Sobolev space H s (T 2 ) as consisting of distributions v defined on T 2 for which the norm defined in frequency space is finite.Here, v(ξ) denotes the Fourier coefficients of v. To shorten notation, we will write • s,2 for • W s,2 (T 2 ) and/or • H s (T 2 ) .When k = s = 0, we get the usual L 2 (T 2 ) space whose norm we will denote by • 2 for simplicity.We will also use a similar convention for norms • p of general L p (T 2 ) spaces for any p ∈ [1, ∞] as well as for the inner product •, • k,2 := •, • W k,2 (T 2 ) when k ∈ N. Additionally, we will denote by W k,p div (T 2 ; R 2 ), the space of weakly divergence-free vector-valued functions in W k,p (T 2 ).Finally, we define the following space (2.5) 2.2.Preliminary estimates.We begin this section with the following result which follow from a direct computation using the definition (2.3) of the Sobolev norms.Lemma 2.2.Let k ∈ N ∪ {0} and assume that the triple (u, w) satisfies Let us now recall some Moser-type calculus.See [22,19,21].
), we have and 2.3.Main results.Our current goal is to construct a solution for the system of equations (1.1)-(1.5).
To do this, we first make the following assumption on our set of data.
Assumption 2.4.Let u h ∈ W 3,2 div (T 2 ; R 2 ) and f ∈ W 2,2 (T 2 ) and assume that (b 0 , ω 0 ) ∈ M. Unless otherwise stated, Assumption 2.4 now holds throughout the rest of the paper.We are now in a position to make precise, exactly what we mean by a solution.
Remark 2.6.We remark that the regularity of the solution (b, ω, T ) and its data together with the integral equations immediately imply that b and ω are differentiable in time.Indeed, by using the fact that W 3,2 (T 2 ) and W 2,2 (T 2 ) are Banach algebras, we immediately deduce from the integral equations for b and ω above that It also follow from the integral equations for the buoyancy and potential vorticity above that the initial conditions are b(0, x) = b 0 (x) and ω(0, x) = ω 0 (x).The corresponding differential forms (1.1)-(1.2) are clearly immediate from the integral representations.
Remark 2.7.Since we are working on the torus, and the velocity fields are defined by (1.4)-(1.5),we have in particular, ´T2 u h dx = 0 and ´T2 u dx = 0. From the latter, we get that q and f have zero averages.Consequently, we will assume that all functions under consideration have zero averages.
Definition 2.8 (Maximal solution).We call (b, ω, T max ) a maximal solution to the system (1.1)-(1.5)if: • there exists an increasing sequence of time steps (T n ) n∈N whose limit is We shall call T max > 0 the maximal time.
Remark 2.9.Condition (2.9) means that the solution breaks down at the limit point T max .
We are now in a position to state our first main result.
Once we have constructed a local solution, we can show that this solution is continuously dependent on its initial state in a more general class of function space.This choice of class appears to be the strongest space in which the analysis may be performed.More details will follow in the sequel but first, we give the statement on the continuity property of strong solutions with respect to its data.
Remark 2.12.For the avoidance of doubt, we make clear that the final estimate for the differences in Theorem 2.11 above is stated in terms of a larger space W 2,2 (T 2 ) × W 1,2 (T 2 ) with smaller norms than the space of existence W 3,2 (T 2 ) × W 2,2 (T 2 ).In order to obtain this bound however, in particular, we require the initial conditions of one of the solution to be bounded in the stronger space of existence, i.e, boundedness of b 1 0 3,2 , ω 1 0 2,2 rather than in the weaker space for which we obtain our final stability estimate.Furthermore, we also require the potential vorticity ω 2 (at all time of existence) of the second solution to be also bounded in the stronger space of existence W 2,2 (T 2 ).Explicitly, the only flexibility we have has to do with the second buoyancy b 2 for which it appears that we are able to relax to live in the weaker space W 2,2 (T 2 ).Unfortunately however, because of the highly coupled nature of the system of equations under study, needing ω 2 ∈ W 2,2 (T 2 ) for all times automatically requires that b 2 ∈ W 3,2 (T 2 ) for all times.Subsequently, this means that we require boundedness of the initial condition of the second solution in the stronger space of existence, i.e, boundedness of b 2 0 3,2 , ω 2 0 2,2 in addition to that of the first solution.The requirement of needing both pair of initial conditions to live in a stronger space is in contrast to simpler looking models like the Euler equation where it suffices to require just having one initial condition to have stronger regularity.Unfortunately, we can not do better by our method of proof (which is to derive estimates for equations solved by the differences, i.e., the energy method) but we do not claim that other methods for deriving analogous estimates may not yield better result either.
As a consequence of Theorem 2.11, the following statement about uniqueness of the strong solution is immediate.
Finally, we can show that a maximal solution of (1.1)-(1.5), in the sense of Definition 2.8, also exists.
The α-TQG model : In the following, we consider a 'regularized' version of (1.1)-(1.5).Since our notion of a solution to (1.1)-(1.5)involves the pair (b, ω) (and not explicitly in terms of u), henceforth, we will use the triple (b α , ω α , T ) to denote the corresponding solution to the following α-TQG model for the avoidance of confusion.To be precise, for fixed α > 0, we will be exploring the pair (b α , w α ) that solves but where now, Note that at least formally, this implies that We also note that even though the velocity fields u in (1.5) and u α in (2.15) are different, the coupled equations (1.1)-(1.2) and (2.12)-(2.13)are exactly of the same form.In the following, we claim that the exact same result, Theorem 2.14 holds for the α-TQG model (2.12)-(2.14)introduced above.
Remark 2.16.The fact that Theorem 2.15 holds true is hardly surprising considering that we are basically treating the exact same model as (1.1)-(1.2) albeit different velocity fields.Heuristically, one observes that for very small α, ψ α ∼ ψ ∼ ψ and thus, both models are the same.Nevertheless, rigorously, the key lies in the estimate (2.7) which still holds true when one replaces u with u α and where now, the triple (u α , ω α , f ) satisfies for a given f having sufficient regularity.Indeed, Theorem 2.15 follow from the following lemma.Lemma 2.17.Under Assumption 2.4, there exists a unique solution (b α , ω α , T α ) of (2.12)-(2.14)for any α > 0.Moreover, (b α , ω α , T α max ) is a unique maximal solution of (2.12)-(2.14).
Proof.The existence of a unique solution (b α , ω α , T α ) of (2.12)-(2.14)for any α > 0 follow from Proposition 3.1, Lemma 3.5 and Corollary 2.13 with T α = T .We will prove all these results later for the TQG and one will observe that nothing changes for the α-TQG.Note that since (1.1)-(1.5)and (2.12)-(2.14)share the same data, Assumption 2.4, the local existence time T is accordingly independent of α > 0.
The construction of the maximal solution (b α , ω α , T α max ) of (2.12)-(2.14)also follow the same gluing argument used in showing the existence of the maximal solution (b, ω, T max ) for (1.1)-(1.5) in Section 3.4.

Construction of strong solutions
In this section, we give a construction of a solution, in the sense of Definition 2.5, to the thermal quasigeostrophic system of equations (1.1)-(1.5).We will achieve this goal by rewriting our set of equations (1.1)-(1.5) in an abstract form and show that the resulting operator acting on the pair (b, ω) satisfies the assumptions of Theorem 8.1 in the appendix.The details are as follows.We consider where the operator Estimates for convective terms.In the following, we let k ∈ N ∪ {0} and recall that W k,p (T 2 ) is a Banach algebra when kp > 2. By using Sobolev embeddings and (2.7), we obtain the following estimates. (1b ) and u ∈ W 3,2 div (T 2 ; R 2 ) solves (1.5) for a given f ∈ W 2,2 (T 2 ) and a given u h ∈ W 3,2 div (T 2 ; R 2 ), then With these estimates in hand, we can now proceed to prove the existence of a local strong solution of (3.1).

3.2.
Proof of local existence.Our proof of a local strong solution to (1.1)-(1.5)will follow from the following proposition.Its proof will follow the argument presented in [20] (See Appendix, Section 8) for the construction of a local solution to the Euler equation.
Proposition 3.1.Take Assumption 2.4.There exists a solution (b, ω, T ) of (3.1) such that: (1) the time T > 0 satisfy the bound (3) the pair (b, ω) satisfies the bound Proof.Let us consider the following triplet {V, H, X} given as follows.Take X := W 1,2 (T 2 ) × L 2 (T 2 ) and let H := W (3,2),2 (T 2 ) × W (2,1),2 (T 2 ) be a Hilbert space endowed with the inner product We note that for (b, ω) ∈ M, we have that and thus, H is equivalent to M. Also, by virtue of the estimates shown in Section 3.1, we can conclude that the operator A is weakly continuous from H into X.
Next, we let V be the domain of an unbounded selfadjoint operator S ≥ 0 in X with domain(S) ⊂ H and Su, v = u, v H for u ∈ domain(S).More precisely, S = |β|≤3 (−∂) β ∂ β subject to periodic boundary conditions.With this definition of V in hand, we can again conclude from the estimates in Section 3.1 that A maps V into H.Our goal now is to show that for any (b, ω) ∈ V , there exists an increasing function ρ(r) ≥ 0 of r ≥ 0 such that To achieve this goal, we further refine the following estimates from Section 3.1.In particular, we recall that if b ∈ W k+1,2 (T 2 ) and ω ∈ W 2,2 (T 2 ) and that u ∈ W 3,2 div (T 2 ; R 2 ) solves (1.5) for a given f ∈ W 2,2 (T 2 ) and a given u h ∈ W 3,2 div (T 2 ; R 2 ), then for k ∈ {0, 1, 2}, In the above estimates, we have used inequalities such as x 1 + x 2 for x ≥ 0 and Young's inequalities.Furthermore, if b ∈ W 3,2 (T 2 ), ω ∈ W 2,2 (T 2 ) and u ∈ W 3,2 div (T 2 ; R 2 ), we also have that We can therefore conclude that for (b, ω) ∈ H (and in particular, for any (b, ω) ∈ V ), (3.12) with a constant depending only on u h 3,2 and f 2,2 .Since ρ(r) = c(1+r) 2 ≥ 0 is an increasing function of r ≥ 0, we can conclude from the result in the Appendix, Section 8, that for a given (b 0 , ω 0 ) ∈ H ≡ M, there exist a solution (b, ω, T ) of (3.1) of class where r is an increasing function on (0, T ).The time T > 0 in (3.13)-(3.14)depends only on u h 3,2 and f 2,2 by way of the function ρ(•) as well as on (b 0 , ω 0 ) M .To be precise, for ρ(r) = c(1 + r) 2 ≥ 0 where c > 0 depends only on u h 3,2 and f 2,2 , we obtain T > 0 by solving the equation which yields We therefore require for a solution of (3.15) to exist.This finishes the proof.
Remark 3.2.Note that r(t) converges to (b 0 , ω 0 ) 2 M as t → 0 + .Therefore, given the bound (3.14) and (3.16), we can conclude that lim sup On the other hand, because the pair (b, ω) is of class (3.13), we can also conclude that and thus, we obtain strong right continuity at time t = 0. Since our system (1.1)-(1.5) is time reversible (see Remark 3.3 below), we also obtain strong left continuity at t = 0 so that we are able to conclude that (3.13) is actually strongly continuous at t = 0. Remark 3.4.We remark that the function ρ and as such, the function r solving (3.15) is not unique.In fact, one can construct countably many (if not infinitely many) of these functions by varying the estimates for the convective terms on the left-hand sides of (3.10)- (3.11).For example, a different choice of Hölder conjugates for estimating the aforementioned terms will lead to a different ρ and r.The importance of this remark lies in the fact that ρ and r determines the longevity of the solution (b, ω, T ) to (1.1)-(1.5).
One can therefore tune T > 0 by modifying ρ and r accordingly.For example, we may also obtain the following for which the corresponding solutions to (3.15) are , where s = s(t) = r(0) respectively.Again, the constants depends only on u h 3,2 and f 2,2 .The second solution r 2 is unsuitable since for one, it is twofold.The first solution r 1 however means that we require a time for a solution of (3.15) to exist.Therefore, if we optimize the constant c > 0 so that they are the same throughout and we compare r 1 with r given in (3.16), we are able to conclude that r yields a longer time of existence of a solution (b, ω, T ) to (1.1)-(1.5)as compared to r 1 since .
In the next section, we show an estimate for the difference of two solutions constructed above after which we are able to strengthen the weak continuity (3.13) to a strong one for all times of existence.
3.3.Difference estimate.In the following, we let (b 1 , ω 1 , T ) and (b 2 , ω 2 , T ) be two solutions of (1.1)-(1.5)sharing the same data u h ∈ W 3,2 div (T 2 ; R 2 ), f ∈ W 2,2 (T 2 ) and (b 0 , ω 0 ).So in particular, sup where We now show the following stability and uniqueness results in the space W 2,2 (T 2 ) × W 1,2 (T 2 ) which is larger than the space W 3,2 (T 2 ) × W 2,2 (T 2 ) of existence of the individual solutions.Indeed, it suffices to show uniqueness in the even larger space W 1,2 (T 2 ) × L 2 (T 2 ) but we avoid doing so since we wish to concurrently show the proof of Theorem 2.11 as well.
Having shown uniqueness, we can conclude that the weakly continuous functions (3.5) are indeed strongly continuous by using time-reversibility of (1.1)-(1.2).
Lemma 3.5.Take the assumptions of Proposition 3.1 to be true.Then Proof.Firstly, we recall that the solution we constructed in Proposition 3.1 is actually strongly continuous at time t = 0, recall Remark 3.2.Now consider a solution (b, ω, T 0 ) of (1.1)-(1.2) where T 0 ∈ [0, T ] is fixed but arbitrary.Then by (3.6), this solution will satisfy the bound ).We can now construct a new solution (b 1 , ω 1 , T 0 + T 1 ) by taking (b, ω)(T 0 ) as an initial condition.A repetition of the argument leading to (3.39) means that this new pair (b 1 , ω 1 ) solving (1.1)-(1.2) will satisfy the inequality where .
Furthermore, (b 1 , ω 1 ) must coincide with (b, ω) on [T 0 , T 0 + T 1 ] ∩ [0, T ] by uniqueness and the fact that they agree at T 0 ∈ [0, T ].Subsequently, we obtain from (3.40), strong right-continuity of (b 1 , ω 1 ) M at time t = T 0 just as was done for t = 0 in Remark 3.2.Again, by uniqueness, this implies that (b, ω) M is also strongly right-continuous at time t = T 0 .Since T 0 was chosen arbitrarily, this means that (b, ω) is strongly right-continuous on [0, T ].Since our system is reversible, we can repeat the above argument for the backward equation from which we obtain strong left-continuity on [0, T ].We have thus shown that weakly continuous solution (3.5) is indeed strongly continuous.
3.4.Maximal solution.We now end the section with a proof of the existence of a maximal solution to the TQG model.
Proof of Theorem 2.14.By Proposition 3.1, we found a time By repeating the argument, we can also find • we obtain an increasing family (T n ))n ∈ N of time steps whose limit is T max ∈ (0, ∞]; • for each n ∈ N, (b, ω, T n ) is a solution of (1.1)-(1.5); since otherwise, we can repeat the procedure above to obtain T > 0 satisfying ), such that (b, ω, T max max ) is a solution of (1.1)-(1.5)with T max max = T max + T .This will contradict the fact that T max is the maximal time.
We will now devote the entirety of this section to the proof of Proposition 4.1.In order to achieve this goal, we first introduce a unique maximal solution (b α , ω α , T α max ) with initial conditions b 0 ∈ W 3,2 (T 2 ) and ω 0 ∈ W 2,2 (T 2 ) of the following "intermediate α-TQG" equation given by where like (2.15), and where ω(and not ω α ) satisfies the original potential vorticity equation (1.2) strongly in the PDE sense.The existence of a unique maximal solution to the above intermediate system, denoted by iα-TQG for short, follow from the proof of existence for the original TQG and it particular, for a uniform-in-α set of data (u h , f, b 0 , ω 0 ), we have the following uniform-in-α estimates for the maximal solution of (4.3)-(4.5).Indeed, exactly as in Theorem 2.15, we also have the following result.
Similar to the uniqueness argument, if we apply ∂ β to the equation for b 12 above where now |β| ≤ 1, we obtain where are such that Collecting the information above yields the estimate For the equation for ω 12 , we aim to derive a square-integrable estimate in space.For this, we first note that u • ∇ω 12 , ω 12 = 0. (4.14) Next, we have that Finally, we also have that (4.16) Therefore, Summing up with the estimate for b 12 then yield (4.18) .Now recall that the velocity fields for the α-TQG, the iα-TQG, and the TQG are given by respectively.As a result, in particular, the difference enjoys two extra order of regularity.Furthermore, since each of these individual velocity fields satisfies the bound (2.7), it immediately follows that (4.23) This means that for the comparison of the TQG and iα-TQG, we can conclude from (4.19) that and by Grönwall's lemma, (note that b 12 0 = 0 and The above gives the decay rate of the difference of the solution to the iα-TQG and TQG equations.Our next goal is to obtain a decay rate for the difference between the iα-TQG and the α-TQG equations.Recall that they share the same initial data.For this, we use the estimate By the triangle inequality, it follows from (4.25) and (4.27) that as α → 0. This ends the proof.
Remark 4.3.We conjecture that Proposition 4.1 may be extended to the stronger space W 2,2 (T 2 ) × W 1,2 (T 2 ) in which we showed the stability and uniqueness of the maximal solution.The cost, however, may be the loss of the corresponding decay rate (4.1).The ideas in the proof of Proposition 6 of Terence Tao's note : local-well-posedness-for-the-euler-equations may suffice.

Conditions for blowup of the α-TQG model
In the following, we give Beale-Kato-Majda [2] type conditions under which we expect the strong solution of the α-TQG to blowup.Theorem 5.1.Fix α > 0. Under Assumption 2.4, let (b α , ω α , T α max ) be a maximal solution of (2.12)- Proof.First of all, fix α > 0 and let T α max > 0 be the maximal time so that lim sup We now suppose that holds yielding a contradiction to (5.1).Before we show the estimate (5.3), we first require preliminary estimates for ∇u α ∞ and ω α 2 .To obtain these estimates, we note that the Fourier multiplier m α (ξ) defined below satisfies the bound for all fixed α > 0 and in particular, the bound may be taken uniformly of all α ≥ 1. Due to (5.4) and the continuous embedding W 3,2 (T 2 ) → W 1,∞ (T 2 ), we can conclude that On the other hand, if we test (2.13) with ω α and use (2.7) for w = (1 − α∆) −1 (ω α − f ) and k = 0, we obtain for a constant depending only on u h 2 and f 2 .It therefore follow from (5.6) and (5.2) that which when combined with the first estimate in (5.5) yields Recall that f ∈ W 2,2 (T 2 ) by assumption.With these preliminary estimates in hand, we can proceed to derive estimates for (b α , ω α ) in the space M of existence of a solution.Since the space of smooth functions is dense in M, it suffices to show our result for a smooth solution pair (b α , ω α ).To achieve our goal, we apply ∂ β to (2.12) for |β| ≤ 3 to obtain where Now since divu α = 0, if we multiply (5.9) by ∂ β b α and sum over the multiindex β so that |β| ≤ 3, we obtain (5.10) M ) where we have used (2.7) for w = (1 − α∆) −1 (ω α − f ) and k = 2. Next, we find a bound for ω α 2 2,2 .For this, we apply ∂ β to (2.13) for |β| ≤ 2 and we obtain where (5.14) holds true.Recall that u h ∈ W 3,2 (T 2 ) by assumption and note that we have used the continuous embedding W 3,2 (T 2 ) → W 2,4 (T 2 ) and the second inequality in (5.5) in order to obtain the term 1+ ω α 2 for the estimates (5.12).Next, by using divu α = 0, we have the following identity Additionally, the following estimates for L 2 inner products holds true (5.17) since u h ∈ W 3,2 (T 2 ).If we now collect the estimates above (keeping in mind that f ∈ W 2,2 (T 2 ) and u h ∈ W 3,2 (T 2 )), we obtain by multiplying (5.11) by ∂ β ω α with |β| ≤ 2 and then summing over |β| ≤ 2, the following Summing up (5.10) and (5.18) and using (5.7)-(5.5)yields so that by Grönwall's lemma and (5.2), we obtain max ,K,f,u h ,ω0,b0,α 1 for all t < T α max contradicting (5.1).Remark 5.2 (Global existence for constant-in-space buoyancy).We note that for the α-TQG model, when the buoyancy is constant in space so that ∇b α = 0, then the equation for the buoyancy decouples from that of the potential vorticity.The potential vorticity equation then reduces to the 2-dimensional Euler equation albeit a different Biot-Savart law.Nevertheless, by inspecting the proof of Theorem 5.1 above, one observes that the same global-in-time result applies.In particular, the residual estimate (5.12) still holds true and we obtain from this, a refine version of (5.7)- (5.8)where exp(cK) = 1.Recall from (5.2) that K = 0 when ∇b α = 0. Remark 5.3 (Global existence for super diffusive buoyancy).Fix α > 0. Note that by adding any super diffusive term Λ β b α , β ≥ 1 to the right-hand side of (2.12) so that holds for all t < T α max , then we obtain a global solution since in this case, we obtain the energy estimate for all t < T α max .Next, we also give an alternative to the blowup condition in Theorem 5.1 above.
Remark 5.6.Note that the original statement in [8] restricted the size of f 1,2 to being at most one.The current form for any size of f 1,2 follows immediately as demonstrated in for example [33,10].With Lemma 5.5 in hand, we can now show the final blowup condition.
Proof.Fix α > 0 and let T α max > 0 be the maximal time so that lim sup holds and thereby yields a contradiction to (5.29).
To show this, we first fix α > 0 and define With this definition in hand, it follows from estimates (5.10) and (5.18) that By using Lemma 5.5 and the monotonic property of logarithms however, we can deduce that Next, we observe that the second estimate in (5.5) yields ).However, if we test (2.13) with ω α , we also obtain for a constant depending only on u h ∞ and f 2 .It therefore follow from (5.30) that ω α (t) holds for t < T α max .In particular, given (5.30), we have shown that T α max ,K,f,u h ,ω0,b0 1 for all t < T α max contradicting (5.29).

Numerical implementation
For numerical implementation, we choose to work in weaker spaces than what the wellposedness theorem dictates.Additionally we allow for boundary conditions in the numerical setup.We recognise these choices create gaps between the theory and the implementation.
6.1.Discretisation methods for α-TQG and TQG.In this subsection, we describe the finite element (FEM) spatial discretisation, and finite difference Runge-Kutta time stepping discretisation methods that are utilised for the α-TQG system.By setting α to zero we obtain our numerical setup for the TQG system.
For D = T 2 , our numerical setup does not change, in which case the boundary flux terms in the discretised equations are set to zero.
6.1.1.The stream function equation.Let H 1 (D) denote the Sobolev W 1,2 (D) space and let .∂D denote the L 2 (∂D) norm.Define the space (6.2) We express (2.14) as two inhomogeneous Helmholtz equations We take ψ, ψ α ∈ W 1 (D).Since the two equations are of the same form, let us first consider (6.3).Using an arbitrary test function φ ∈ W 1 (D), we obtain the following weak form of (6.3), for v, φ ∈ W 1 (D), then (6.5) can be written as And similarly for (6.4), we have (6.9) We discretise equations (6.8) and (6.9) using a continuous Galerkin (CG) discretisation scheme.Let δ be the discretisation parameter, and let D δ denote a space filling triangulation of the domain, that consists of geometry-conforming nonoverlapping elements.Define the approximation space (6.10) is the space of continuous functions on D, and Π k (K) denotes the space of polynomials of degree at most k on element K ∈ D δ .
For (6.8), given f δ ∈ W k δ (D) and ω δ ∈ V k δ (D) (see (6.13) for the definition of V k δ (D)), our numerical approximation is the solution ψδ ∈ W k δ (D) that satisfies (6.11) for all test functions φ δ ∈ W k δ (D).Then, using ψδ , our numerical approximation of ψ α is the solution for all test functions φ δ ∈ W k δ (D).For a detailed exposition of the numerical algorithms that solves the discretised problems (6.11) and (6.12) we point the reader to [13,7].6.1.2.Hyperbolic equations.We choose to discretise the hyperbolic buoyancy (2.12) and potential vorticity (2.13) equations using a discontinuous Galerkin (DG) scheme.For a detailed exposition of DG methods, we refer the interested reader to [14].
Define the DG approximation space, denoted by V k δ (D), to be the element-wise polynomial space, (6.13) We look to approximate b α and ω α in the space V k δ (D).Essentially, this means our approximations of b α and ω α are each an direct sum over the elements in D δ .Additional constraints on the numerical fluxes across shared element boundaries are needed to ensure conservation properties and stability.Further, note that W k δ (D) ⊂ V k δ (D).This inclusion is also needed for ensuring numerical conservation laws, see Section 6.1.3.
For the buoyancy equation (2.12), we obtain the following variational formulation where ν δ ∈ V k δ (D) is any test function, ∂K denotes the boundary of K, and n denotes the unit normal vector to ∂K.Let b α δ be the approximation of b α in V k δ , and let u α δ = ∇ ⊥ ψ α δ for ψ α δ ∈ W k δ .Our discretised buoyancy equation over each element is given by (6.15) be the approximation of ω α , and let h δ ∈ W k δ .We obtain the following discretised variational formulation that corresponds to (2.13), for test function ν δ ∈ V k δ (D).At this point, we only have the discretised problem on single elements.To obtain the global approximation, we sum over all the elements in D δ .In doing so, the ∂K terms in (6.15) and (6.16) must be treated carefully.Let ∂K ext denote the part of cell boundary that is contained in ∂D.Let ∂K int denote the part of the cell boundary that is contained in the interior of the domain D\∂D.On ∂K ext we simply impose the PDE boundary conditions.However, on ∂K int we need to consider the contribution from each of the neighbouring elements.By choice, the approximant ψ α δ ∈ W k δ is continuous on ∂K.And since (6.17) where τ denotes the unit tangential vector to ∂K, u α δ • n is also continuous.This means u α δ • n in (6.15) (also in (6.16)) is single valued.However, due to the lack of global continuity constraint in the definition of V α k (D), b α δ and ω α δ are multi-valued on ∂K int .Thus, as our approximation of ω α and b α over the whole domain is the sum over K ∈ D δ , we have to constrain the flux on the set K∈D δ ∂K \ ∂D.This is done using appropriately chosen numerical flux fields in the boundary terms of (6.15) and (6.16).
Let ν − := lim ↑0 ν(x + n) and ν + := lim ↓0 ν(x + n), for x ∈ ∂K, be the inside and outside (with respect to a fixed element K) values respectively, of a function ν on the boundary.Let f be a numerical flux function that satisfies the following properties: (i) consistency stable in the enstrophy norm with respect to the buoyancy equation, see [4,Section 6].
With such an f , we replace b (6.15).Similarly, in (6.16), we replace (ω For a general nonlinear conservation law, one has to solve what is called the Riemann problem for the numerical flux, see [14] for details.In our setup, we use the following local Lax-Friedrichs flux, which is an approximate Riemann solver, ( where (6.21) Finally, our goal is to find b α δ , ω α δ ∈ V k δ (D) such that for all ν δ ∈ V k δ (D) we have being the numerical approximation to the stream function.Remark 6.2.In (6.22) and (6.23) we do not explicitly distinguish ∂K ext and ∂K int because for the boundary conditions (6.1), the ∂K ext terms vanish.6.1.3.Numerical conservation.Conservation properties of the TQG system was first shown in [17].More specifically, the TQG system conserves energy, and an infinite family of quantities called casimirs.Proposition 6.3 below describes the conservation properties of the α-TQG system.We note that the form of the conserved energy and casimirs are the same as that of the TQG system.Although the result is stated for the system with boundary conditions, it is easy to show that the same result holds for T 2 .Proposition 6.3 (α-TQG conserved quantities).On a bounded domain D with boundary ∂D, consider the α-TQG system (2.12) -(2.14) with boundary conditions We have Proof.For the energy E α (t), we obtain Substitute ∂ t ω α and ∂ t b α using (2.12) and (2.13), the result then follows from direct calculations.
Similarly for the casimirs, we obtain the result by directly evaluating d t C α Ψ,Φ (t).In our FEM discretisation, we impose only Dirichlet boundary conditions.If the integral Neumann boundary conditions are not imposed, then for Proposition 6.3 to hold, we necessarily require ψ α and ∆ψ α on ∂D to be zero.
We now analyse the conservation properties of our spatially discretised α-TQG system.Let •, • H 1 (D) be the H 1 (D) inner product.Define the numerical total energy (6.28)In other words, the semi-discrete discretisation conserves energy in the TQG case.
6.2.α-TQG linear thermal Rossby wave stability analysis.A dispersion relation for the TQG system was derived in [17].There, the authors showed that the linear thermal Rossby waves of the TQG system possess high wavenumber instabilities.More specifically, the Doppler-shifted phase speed of these waves becomes complex at sufficiently high wavenumbers.However, the growth rate of the instability decreases to zero as |k| −1 , in the limit that |k| → ∞.Consequently, the TQG dynamics is linearly wellposed.That is, the TQG solution depends continuously on initial conditions.In this subsection, using the same equilibrium state as in [17], we derive a dispersion relation for the thermal Rossby wave solutions of the linearised α-TQG system.Again the Doppler-shifted phase speed of these waves becomes complex at sufficiently high wavenumbers.However, in this case, the growth rate of the instability decreases to zero as |k| −2 , in the limit that |k| → ∞.In the limit that α tends to zero, one recovers the TQG dispersion relation derived in [17].
Substituting these solutions into the linearised equations we have where k is the first component of k.Upon setting α = 0 in (6.44), we recover the Doppler-shifted phase speed of thermal Rossby waves for the TQG system in [17].When α = 0, B = 0, and U = 0 in (6.44), the remaining dispersion relation differs slightly from the dispersion relation for QG Rossby waves in non-dimensional variables.This is because the Casimirs in the Hamiltonian formulation of the TQG system differ from the those of standard QG.As discussed in [17] for the TQG system, when Y > 0, the Doppler shifted phase velocity in (6.44) for the linearised wave motion becomes complex at high wavenumber |k| 1; namely for (|k| 2 + 1)(α|k| 2 + 1) ≥ X 2 /(4Y ), for both branches of the dispersion relation.However, the growth-rate of the instability found from the imaginary part of C(α) in (6.44) decays to zero as O(|k| −2 ), for |k| 1.Indeed, the linear stability analysis leading to C(α) predicts that the maximum growth rate occurs at a finite wavenumber |k| max which depends on the value of α, beyond which the growth rate of the linearised TQG wave amplitude falls rapidly to zero.One would expect that simulated solutions of α-TQG would be most active at the length-scale corresponding to |k| max , at which the linearised thermal Rossby waves are the most unstable.If this maximum activity length-scale is near the grid truncation size, then numerical truncation errors could cause additional numerical stability issues!We have experience such problems during our testing of the numerical algorithm for the TQG system.Even for equilibrium solution "sanity" check tests, we have found that unless the time step was taken to be incredibly small, high wavenumber truncation errors have caused our numerical solutions to eventually blow up.
The magnitude of α controls the unstable growth rate of linearised TQG waves at asymptotically high wavenumbers.That is, the presence of α regularises the wave activity at high wavenumbers, see Figure 3. Thus, from the perspective of numerics, the α regularisation is available to control numerical problems that may arise from the inherent model instability at high wavenumbers, without the need for additional dissipation terms in the equation.6.3.Numerical example.The numerical setup we consider for this paper is as follows.The spatial domain D is an unit square with doubly periodic boundaries.We chose to discretise D using a grid that consists of 256 × 256 cells, i.e. cardinality(D δ ) = 256 × 256.This was the maximum resolution we could computationally afford, to obtain results over a reasonable amount of time.
We computed each solution for 5000 time steps1 .Figures 5, 6 and 7 show α-TQG solution snapshots of buoyancy, potential vorticity and velocity magnitude respectively, at the 1280'th and 2600'th time steps, for α values 0, 1 128 2 , 1 64 2 and 1 16 2 .The interpretation of the regularisation parameter α is that its square root value corresponds to the fraction of the domain's length scale that get regularised.At the 1280'th time step, the flows are in early spin-up phase.At the 2600'th time step, although the flows are still in spin-up phase, much more flow features have developed.The sub-figures illustrate via comparisons, how α controls the development of small scale features and instabilities.Increasing α leads to more regularisation at larger scales.
In view of Proposition 4.1, we investigated numerically the convergence of α-TQG to TQG using our numerical setup.Consider the relative error between TQG buoyancy and α-TQG buoyancy, and the relative error between TQG potential vorticity and α-TQG potential vorticity, The norms were chosen in view of Proposition 4.1.Figures 8a and 8b show plots of e b (t, α) and e ω (t, α) as functions of time only, for fixed α values 1 16 2 , 1 32 2 , 1 64 2 , 1 128 2 and 1 220 2 , from the initial time up to and including the 2800'th time step.We observe that starting from 0, the relative errors e b and e ω initially increase over time but plateau at around and beyond the 2000'th time step.Up to the 1400'th time step, the plotted relative errors remain less than 1.0, and are arranged in the ascending order of α, i.e. smaller α's give smaller relative errors.
We note that, if we compare and contrast the relative error results with the solution snapshots of buoyancy and potential vorticity, particularly at the 2600'th time step (shown in Figures 5b and 6b), we observe "discrepancies" -for lack of a better term -between the results.For example, in Figure 5b, if we compare the α = 1 128 2 snapshot with the α = 1 16 2 one, the former show small scale features that are much closer to those that exist in the reference TQG α = 0 solution.However, according to the relative error measurements, the two regularised solutions are more or less equivalent in their differences to the reference TQG solution at time step 2600.
If we fix the value of the time parameter in e b (t, α) and e ω (t, α), and vary α, we can estimate the convergence rate of α-TQG to TQG for our numerical setup.Proposition 4.1 predicts a convergence rate of 1 in α, using the H 1 norm on the buoyancy differences and the L 2 norm on the potential vorticity differences.Note that, in (4.1), the left hand side evaluates a supremum over a given compact time interval.We see in Figures 8a and 8b that the relative errors are monotonic up to around the 1500'th time step mark.So for a given time point T between the initial time and the 1500'th time step, we assume we can evaluate e b and e ω at T to estimate the supremum over [0, T ].
Figures 9a and 9b show plots -in log-log scale -of e b and e ω as functions of α only at the 600'th, 800'th, 1000'th and 1200'th time steps, for all the chosen α values.These numbers of time steps correspond to T = 0.3, T = 0.4, T = 0.5 and T = 0.6 respectively.We also plotted in Figures 9a and 9b linear functions of α to provide reference order 1 slopes.Comparing the results to the reference, we see that the theoretical convergence rate of 1 is attained for T = 0.3 (600'th time step), T = 0.4 (800'th time step) and T = 0.5 (1000'th time step).However for T = 0.6 (1200'th time step), the convergence rates are at best order 1/2.

Conclusion and outlook
In conclusion, we have formulated the thermal quasi-geostrophic (TQG) model equations in the regime of approximations relevant to GFD and we have explored their analytical and numerical properties.In respect to the analytical properties, we have shown that the TQG model and its α-regularized version, the α-TQG model, admit unique local strong solutions that are stable in a larger space with weaker norm, that both models have a maximum time of existence, solutions of the latter model converge to solutions of the former model on any time interval in which a solution to the former lives, and we have also determined conditions under which the solution of the α-TQG will blow up.
With respect to numerics, we described our discretisation methods for approximating TQG and α-TQG solutions and showed that the FEM semi-discrete scheme conserves numerical energy in the TQG case.Using example simulation results, we numerically verified that α-TQG solutions converge to that of the TQG system and attains the order 1 in α convergence rate predicted by Proposition 4.1.Additionally, we derived a dispersion relation for the linearised α-TQG system and showed how α-regularisation could be used to control the development of high wavenumber instabilites.Given that we have convergence of α-TQG solutions to TQG, the linear stability results can be viewed as generalisations of those shown in [17] for the TQG system.
We now end this section with some open problems.approach to the Hamiltonian formulation of the TQG equations results in a stochastic Hamiltonian formulation of TQG equations in (1.6)- (1.8).This stochastic version of TQG preserves an infinite family of integral conserved quantities, as is shown in [17].The SALT version of TQG represents an outstanding challenge for uncertainty quantification and data assimilation which will surely spur us on to further investigations of these problems.• As also discussed in Section 1.3, the parallels between TQG and the Rayleigh-Bénard equations should also attract our attention for the consideration of a SALT version of the deterministic Rayleigh-Bénard convection equations.

Appendix
We present a result in this section, which along with its proof, can be found in [20].Our construction of local solutions relies on this result and thus, we present it here for the sake of completeness.We consider the abstract equation ∂ ∂t u + A(t, u) = 0, t ≥ 0, u t=0 = u 0 , (8.1) where A is a nonlinear operator.Theorem 8.1.Let {V, H, X} be real separable Banach spaces such that • the embeddings V → H → X are continuous and dense; • H is a Hilbert space; • There is a continuous, nondegenerate bilinear form ( , ) on X × V such that (u, v) = u, v H for u ∈ H and v ∈ V .Let A be a (sequentially) weakly continuous map on [0, T * ] × H into X such that where ρ(r) ≥ 0 is a monotone increasing function of r ≥ 0. Then for any u 0 ∈ H, there exists T > 0, T ≤ T * , and a solution u of (8.1) in the class T may be any value such that r exists on [0, T ].If the solution to (8.5) is not unique, r should be the maximal soluton.(c) u(t) → u 0 holds strongly in H as t → 0. In other words, the solution (8.3) is strongly continuous at the initial time t = 0.     , and e ω (t, α) as functions of α only, at the fixed time values t = 0.3 (equivalently at the 600'th time step), t = 0.4 (equivalently at the 800'th time step), t = 0.5 (equivalently at the 1000'th time step) and t = 0.6 (equivalently at the 1200'th time step).In view of Proposition 4.1, we assumed that, based on the initial monotonicity of the relative errors (see Figure 8), the supremum over [0, T ] in (4.1) can be estimated by evaluating e b and e ω at T .Proposition 4.1 predicts that up to a certain time, the rate of convergence of e b and e ω to 0 should be no less than order 1 in α.In both sub-figures we have plotted, as the reference for comparison, linear functions of α.Thus, comparing slopes to the reference, we see that numerically we get order 1 convergence for t = 0.3, t = 0.4 and t = 0.5.However, for t = 0.6, the rate of convergence is no more than 1/2.

Figure 2 .
Figure 2. The tree of model derivations can function as a roadmap of geophysical fluid dynamics.The solid red arrows indicate the sequence of approximations that lead to the thermal quasi-geostrophic equations.

Definition 2 . 5 (
Local strong solution).Let T > 0 be a constant.We call the triple (b, ω, T ) a local strong solution or simply, a local solution or a solution to the system (1.1)-(1.5)if the following holds.•The buoyancy b satisfies b ∈ C([0, T ]; W 3,2 (T 2 )) and the equation b

Figure 4 .
Figure 4. Initial conditions -buoyancy (6.46) on the left, and potential vorticity (6.45) on the right.The different colours of the PV plot correspond to positive and negative values of PV, which can be interpreted as clockwise and anticlockwise eddies respectively.

•
Can one construct a global-in-time strong solution (or a global weak solution) of either the TQG model equations (1.1)-(1.5),or the α-TQG model model equations (2.12)-(2.14)?• Can one give a Beale-Kato-Majda condition for the blowup of a strong solution to the TQG in terms of a single unknown?That is, is there a BKM condition in terms of either b or ω (or u) that does not require a combination of both variables?• As discussed in Section 1.3, an application of the Stochastic Advection by Lie Transport (SALT)

Figure 5 .
Figure 5. Comparisons of solution snapshots of the α-TQG buoyancy field that correspond to four different values of α, at two different points in time.Potential vorticity and velocity magnitude snapshots of the same solutions are shown in Figure 6 and Figure 7 respectively.Shown in each sub-figure are the results corresponding to α = 0 (top left), α = 1 128 2 (top right), α = 1 64 2 (bottom right) and α = 1 16 2 (bottom left).When α = 0, the solution is of the TQG system.Sub-figure (A) shows the solutions at the 1280'th time step, or equivalently when t = 0.64.Sub-figure (B) shows the solutions at the 2600'th time step, or equivalently when t = 1.30.The flows in both sub-figures are at different stages of spin-up, with the flow at t = 1.30showing more developed features.The α parameter is interpreted as the fraction of the domain's length-scale squared value at which regularisation is applied.Due to regularisation, the flows do develop differently.Nevertheless, we observe that as α gets smaller, the α-TQG flow features converge to that of the TQG flow.

Figure 6 .
Figure 6.Comparisons of solution snapshots of the α-TQG potential vorticity field that correspond to four different values of α, at two different points in time.See the caption of Figure5for explanations of the arrangements of these plots.Buoyancy and velocity magnitude snapshots of the same solutions are shown in Figure5and Figure7respectively.As in Figure4, the different colours of the PV field correspond to positive and negative values of PV, which can be interpreted as clockwise and anticlockwise eddies respectively.In sub-figure (B), we observe more clearly the regularisation effects of α.As α increases, the flows develop in different ways -features at smaller scales no longer develop, and larger features evolve differently as a result.Nevertheless, we observe that as α gets smaller, the α-TQG flow features converge to that of the TQG flow.

Figure 7 .
Figure 7. Comparisons of solution snapshots of the α-TQG velocity magnitudes that correspond to four different values of α, at two different points in time.See the caption of Figure5for explanations of the arrangements of these plots.Buoyancy and potential vorticity snapshots of the same solutions are shown in Figure5and Figure6respectively.Using the same scale for colouring, we observe in both sub-figures the strength of the colours weaken as α increases.This indicates smoothing of large velocity magnitudes.Additionally, in sub-figure (B), we observe how features at smaller scales get smoothed out.In particular, in the α =1  16 2 plot, we see that essentially only large scale vorticies remain.Hence, considering the velocity field as a part of the system's nonlinear advection operator, these figures help us to better visualise how the flow features developed in Figure5and Figure6.

Figure 8 . 2 .
Figure 8. Sub-figures (A) and (B) show, respectively, plots of the relative error functions e b (t, α) and e ω (t, α) as functions of time only, for the fixed α values 1 16 2 , 1 32 2 , 1 64 2 , 1 128 2 and 1 220 2 .See equations (6.49) and (6.50) for the definitions of e b and e ω respectively.The plots are shown from t = 0 up to and including t = 1.4,which is equivalent to 2800 time steps using ∆t = 0.0005.We observe that, starting from 0, all plots of e b and e ω increase initially, and plateau at around the 2000'th time step.Further, up to the 1400'th time step the plotted relative errors remain less than 1.0 and are arranged in the ascending order of α.Relating these relative error results at the 2600'th time step to the solution snapshots at the same time point (shown in figures 5b and 6b), we see that although the snapshots show convergence of flow features, the relative error values suggest all the α-TQG solutions are more or less equally far away from the TQG flow.

Figure 9 .
Figure 9. Sub-figures (A) and (B) show, respectively, plots in log-log scale of the relative error functions e b (t, α), and e ω (t, α) as functions of α only, at the fixed time values t = 0.3 (equivalently at the 600'th time step), t = 0.4 (equivalently at the 800'th time step), t = 0.5 (equivalently at the 1000'th time step) and t = 0.6 (equivalently at the 1200'th time step).In view of Proposition 4.1, we assumed that, based on the initial monotonicity of the relative errors (see Figure8), the supremum over [0, T ] in (4.1) can be estimated by evaluating e b and e ω at T .Proposition 4.1 predicts that up to a certain time, the rate of convergence of e b and e ω to 0 should be no less than order 1 in α.In both sub-figures we have plotted, as the reference for comparison, linear functions of α.Thus, comparing slopes to the reference, we see that numerically we get order 1 convergence for t = 0.3, t = 0.4 and t = 0.5.However, for t = 0.6, the rate of convergence is no more than 1/2.