# Theoretical and computational analysis of the thermal quasi-geostrophic
model

D Crisan, DD Holm, E Luesink, PR Mensah, W Pan  DEPARTMENT OF MATHEMATICS,
IMPERIAL COLLEGE, LONDON SW7 2AZ, UK<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

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 regularized version $\alpha$-TQG converge
to solutions of TQG as its smoothing parameter $\alpha\rightarrow 0$ and we
obtain blowup criteria for the $\alpha$-TQG model. The main contribution of
the computational analysis is to verify the rate of convergence of
$\alpha$-TQG solutions to TQG solutions as $\alpha\rightarrow 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_{\mathrm{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_{\mathrm{max}}$. Should $t<T_{\mathrm{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 $\alpha$-TQG model. This model is constructed in a
similar manner as the $\alpha$ model for the Euler and Navier–Stokes
equations, see [11, 12, 23]. The $\alpha$-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 $\alpha$-TQG
solution converges to the TQG solution as $\alpha\rightarrow 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
$\alpha$-parameter on a time interval where both the TQG solution as well as
the $\alpha$-TQG solution are shown to exist for any $\alpha>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 $\alpha$-TQG
model, if either:

  1. (1)

the $L^{\infty}(\mathbb{T}^{2})$-norm of the buoyancy gradient $\nabla b$
blows up at $T_{\max}<\infty$;

  2. (2)

the $L^{\infty}(\mathbb{T}^{2})$-norm of the velocity gradient
$\nabla\mathbf{u}$ blows up at $T_{\max}<\infty$ or that;

  3. (3)

the $W^{1,2}(\mathbb{T}^{2})$-norm of the buoyancy gradient $\nabla b$ blows
up at $T_{\max}<\infty$.

The contribution of this paper from a numerical perspective is to describe our
spatial and temporal discretisation methods for approximating TQG and
$\alpha$-TQG solutions, and to analyse aspects of numerical conservation
properties with respect to the theoretical conserved quantities of the
$\alpha$-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 $\alpha$-TQG system. Given
that we have convergence of $\alpha$-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.

${{\begin{matrix}\text{Rotating, stratified}\\\ \text{Euler
equations}\end{matrix}}}$${{\begin{matrix}\text{Primitive}\\\
\text{equations}\end{matrix}}}$${{\begin{matrix}\text{Euler--Boussinesq}\\\
\text{equations}\end{matrix}}}$${{\begin{matrix}\text{Thermal rotating}\\\
\text{shallow water}\end{matrix}}}$${{\begin{matrix}\text{Thermal rotating}\\\
\text{Green--Naghdi equations}\end{matrix}}}$${{\begin{matrix}\text{Thermal
rotating}\\\ \text{great lake equations}\end{matrix}}}$Thermal Lagrangian
1${{\begin{matrix}\text{Thermal rotating}\\\ \text{lake
equations}\end{matrix}}}$${{\begin{matrix}\text{Thermal quasi-geostrophic}\\\
\textbf{equations}\end{matrix}}}$${{\begin{matrix}\text{Saint-Venant}\\\
\text{family}\end{matrix}}}$${{\begin{matrix}\text{Camassa-Holm}\\\
\text{family}\end{matrix}}}$$\scriptstyle{\definecolor{.}{rgb}{1,0.2,0.2}\color[rgb]{1,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.2,0.2}\begin{matrix}\text{small
buoyancy}\\\
\text{stratification}\end{matrix}}$$\scriptstyle{\definecolor{.}{rgb}{1,0.2,0.2}\color[rgb]{1,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.2,0.2}\begin{matrix}\text{vertical
average,}\\\ \text{small wave
amplitude}\end{matrix}}$$\scriptstyle{\begin{matrix}\text{vertical
average,}\\\ \text{very small wave
amplitude}\end{matrix}}$$\scriptstyle{\definecolor{.}{rgb}{1,0.2,0.2}\color[rgb]{1,0.2,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.2,0.2}\begin{matrix}\text{vertical
average,}\\\ \text{small wave amplitude}\end{matrix}}$ hydrostatic
---
approx.
restrict to 1D asymptotic expansion around
---
thermal geostrophic balance
restrict to 1D hydrostatic
---
approx.
rigid lid
---
hydrostatic
---
approx.
strict asymptotics hydrostatic
---
approx.
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.

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 $\mathbb{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

(1.1) $\displaystyle\frac{\partial}{\partial t}b+\mathbf{u}\cdot\nabla b$
$\displaystyle=0,$ (1.2) $\displaystyle\frac{\partial}{\partial
t}\omega+\mathbf{u}\cdot\nabla(\omega-b)$
$\displaystyle=-\mathbf{u}_{h_{1}}\cdot\nabla b,$ (1.3)
$\displaystyle\omega=(\Delta-1)\psi+f_{1}$ $\displaystyle\,,\quad{\rm
div}\mathbf{u}=0\,.$

Here $b$ is the vertically averaged buoyancy, $\mathbf{u}$ is the thermal
geostrophically balanced velocity field, $\omega$ is the potential vorticity,
$\psi$ 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 $\mathbf{u}$ and the
streamfunction $\psi$ are related by the equation

(1.4) $\mathbf{u}=\nabla^{\perp}\psi,$

where $\nabla^{\perp}=(-\partial_{y},\partial_{x})$. The vector field
$\mathbf{u}_{h_{1}}$ is defined by

(1.5) $\mathbf{u}_{h_{1}}:=\frac{1}{2}\nabla^{\perp}h_{1}.$

Alternatively, one can formulate the thermal quasi-geostrophic (TQG) equations
in vorticity-streamfunction form. This formulation is given by

(1.6) $\displaystyle\frac{\partial}{\partial t}b+J(\psi,b)$ $\displaystyle=0,$
(1.7) $\displaystyle\frac{\partial}{\partial t}\omega+J(\psi,\omega-b)$
$\displaystyle=-\frac{1}{2}J(h_{1},b),$ (1.8) $\displaystyle\omega$
$\displaystyle=(\Delta-1)\psi+f_{1}.$

The operator $J(a,b)=\nabla^{\perp}a\cdot\nabla b=a_{x}b_{y}-b_{x}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(\psi,a)=\mathbf{u}\cdot\nabla a$, in
which $\psi$ is the streamfunction associated to the velocity vector field
$\mathbf{u}$. The scalar functions $f_{1}$ and $h_{1}$ relate to the usual
Coriolis parameter and bathymetry profile in the following way

(1.9) $\displaystyle h(\mathbf{x})$ $\displaystyle=1+{\rm
Ro}\,h_{1}(\mathbf{x}),$ $\displaystyle f(\mathbf{x})$ $\displaystyle=1+{\rm
Ro}\,f_{1}(\mathbf{x}),$

where ${\rm 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 $\beta$-plane approximation provided that the boundary conditions
are appropriate. Namely, on the $\beta$-plane, one requires that $\beta
f_{0}^{-1}=\mathcal{O}({\rm Ro})$ and $f_{1}(\mathbf{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 $\zeta$ and buoyancy $b$ via the
definition $\psi:=\zeta+\frac{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}(\mathbf{x})$ and
$f_{1}(\mathbf{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) $\displaystyle\frac{\partial}{\partial t}\omega+J(\psi,\omega)$
$\displaystyle=\alpha gT_{z},$ $\displaystyle\frac{\partial}{\partial
t}T+J(\psi,T)$ $\displaystyle=0,$ $\displaystyle\omega$
$\displaystyle=\Delta\psi,$

where $\omega$ denotes the vorticity, $\psi$ is the stream function, $T$ is
the temperature, $g$ is gravity and $\alpha$ 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) $\displaystyle\frac{\partial}{\partial t}\omega+J(\psi,\omega)$
$\displaystyle=\frac{1}{2}J(\psi,b)+\frac{1}{2}J(\zeta-h,b),$
$\displaystyle\frac{\partial}{\partial t}b+J(\psi,b)$ $\displaystyle=0,$
$\displaystyle\omega$ $\displaystyle=\Delta\psi-\psi+f.$

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 $\alpha$-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
$\alpha$-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 $\alpha$-TQG
models converges strongly with $\alpha\rightarrow 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 $\alpha$-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 $\alpha$-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 $\alpha\hbox{-TQG}\to\hbox{TQG}$
derived in Section 4.

##  2\. 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 $\mathbb{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:=\mathbf{x}=(x,y)\in\mathbb{T}^{2}$ on the $2$-torus $\mathbb{T}^{2}$ and a
time variable $t\in[0,T]$ where $T>0$. For functions $F$ and $G$, we write
$F\lesssim G$ if there exists a generic constant $c>0$ such that $F\leq c\,G$.
We also write $F\lesssim_{p}G$ if the constant $c(p)>0$ depends on a variable
$p$. The symbol $|\cdot|$ may be used in four different context. For a scalar
function $f\in\mathbb{R}$, $|f|$ denotes the absolute value of $f$. For a
vector $\mathbf{f}\in\mathbb{R}^{2}$, $|\mathbf{f}|$ denotes the Euclidean
norm of $\mathbf{f}$. For a square matrix $\mathbb{F}\in\mathbb{R}^{2\times
2}$, $|\mathbb{F}|$ shall denote the Frobenius norm
$\sqrt{\mathrm{trace}(\mathbb{F}^{T}\mathbb{F})}$. Finally, if
$S\subseteq\mathbb{R}^{2}$ is a (sub)set, then $|S|$ is the $2$-dimensional
Lebesgue measure of $S$.
For $k\in\mathbb{N}\cup\\{0\\}$ and $p\in[1,\infty]$, we denote by
$W^{k,p}(\mathbb{T}^{2})$, the Sobolev space of Lebesgue measurable functions
whose weak derivatives up to order $k$ belongs to $L^{p}(\mathbb{T}^{2})$. Its
associated norm is

(2.1) $\displaystyle\|v\|_{W^{k,p}(\mathbb{T}^{2})}=\sum_{|\beta|\leq
k}\|\partial^{\beta}v\|_{L^{p}(\mathbb{T}^{2})},$

where $\beta$ is a $2$-tuple multi-index of nonnegative integers of length
$|\beta|\leq k$. The Sobolev space $W^{k,p}(\mathbb{T}^{2})$ is a Banach
space. Moreover, $W^{k,2}(\mathbb{T}^{2})$ is a Hilbert space when endowed
with the inner product

(2.2) $\displaystyle\langle
u,v\rangle_{W^{k,2}(\mathbb{T}^{2})}=\sum_{|\beta|\leq
k}\langle\partial^{\beta}u\,,\,\partial^{\beta}v\rangle,$

where $\langle\cdot\,,\,\rangle$ denotes the standard $L^{2}$-inner product.
In general, for $s\in\mathbb{R}$, we will define the Sobolev space
$H^{s}(\mathbb{T}^{2})$ as consisting of distributions $v$ defined on
$\mathbb{T}^{2}$ for which the norm

(2.3)
$\displaystyle\|v\|_{H^{s}(\mathbb{T}^{2})}=\bigg{(}\sum_{\xi\in\mathbb{Z}^{2}}\big{(}1+|\xi|^{2}\big{)}^{s}|\widehat{v}(\xi)|^{2}\bigg{)}^{\frac{1}{2}}\equiv\|v\|_{W^{s,2}(\mathbb{T}^{2})}$

defined in frequency space is finite. Here, $\widehat{v}(\xi)$ denotes the
Fourier coefficients of $v$. To shorten notation, we will write
$\|\cdot\|_{s,2}$ for $\|\cdot\|_{W^{s,2}(\mathbb{T}^{2})}$ and/or
$\|\cdot\|_{H^{s}(\mathbb{T}^{2})}$. When $k=s=0$, we get the usual
$L^{2}(\mathbb{T}^{2})$ space whose norm we will denote by $\|\cdot\|_{2}$ for
simplicity. We will also use a similar convention for norms $\|\cdot\|_{p}$ of
general $L^{p}(\mathbb{T}^{2})$ spaces for any $p\in[1,\infty]$ as well as for
the inner product
$\langle\cdot,\cdot\rangle_{k,2}:=\langle\cdot,\cdot\rangle_{W^{k,2}(\mathbb{T}^{2})}$
when $k\in\mathbb{N}$. Additionally, we will denote by
$W^{k,p}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, the space of weakly
divergence-free vector-valued functions in $W^{k,p}(\mathbb{T}^{2})$. Finally,
we define the following space

(2.4) $\displaystyle\mathcal{M}=W^{3,2}(\mathbb{T}^{2})\times
W^{2,2}(\mathbb{T}^{2})$

endowed with the norm

(2.5)
$\displaystyle\|(b,\omega)\|_{\mathcal{M}}:=\|b\|_{3,2}+\|\omega\|_{2,2}.$

###  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\in\mathbb{N}\cup\\{0\\}$ and assume that the triple $(\mathbf{u},w)$
satisfies

(2.6) $\displaystyle\mathbf{u}=\nabla^{\perp}\psi,\qquad w=(\Delta-1)\psi.$

If $w\in W^{k,2}(\mathbb{T}^{2})$, then the following estimate

(2.7) $\displaystyle\|\mathbf{u}\|_{k+1,2}^{2}$
$\displaystyle\lesssim\|w\|_{k,2}^{2}$

holds.

Let us now recall some Moser-type calculus. See [22, 19, 21].

######  Lemma 2.3 (Commutator estimates).

Let $\beta$ be a $2$-tuple multi-index of nonnegative integers such that
$|\beta|\leq k$ holds for $k\in\\{1,2\\}$. Let $p,p_{2},p_{3}\in(1,\infty)$
and $p_{1},p_{4}\in(1,\infty]$ be such that

$\displaystyle\frac{1}{p}=\frac{1}{p_{1}}+\frac{1}{p_{2}}=\frac{1}{p_{3}}+\frac{1}{p_{4}}.$

For $u\in W^{k,p_{3}}(\mathbb{T}^{2})\cap W^{1,p_{1}}(\mathbb{T}^{2})$ and
$v\in W^{k-1,p_{2}}(\mathbb{T}^{2})\cap L^{p_{4}}(\mathbb{T}^{2})$, we have
and

(2.8)
$\left\|\partial^{\beta}(uv)-u\partial^{\beta}v\right\|_{p}\lesssim\left(\|\nabla
u\|_{p_{1}}\|v\|_{k-1,p_{2}}+\|u\|_{k,p_{3}}\|v\|_{p_{4}}\right).$

###  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 $\mathbf{u}_{h}\in W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$
and $f\in W^{2,2}(\mathbb{T}^{2})$ and assume that
$(b_{0},\omega_{0})\in\mathcal{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.

######  Definition 2.5 (Local strong solution).

Let $T>0$ be a constant. We call the triple $(b,\omega,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\in C([0,T];W^{3,2}(\mathbb{T}^{2}))$ and the
equation

$\displaystyle b(t)$
$\displaystyle=b_{0}-\int_{0}^{t}\mathrm{div}(b\mathbf{u})\,\,\mathrm{d}\tau,$

holds for all $t\in[0,T]$;

  * •

the potential vorticity $\omega$ satisfies $\omega\in
C([0,T];W^{2,2}(\mathbb{T}^{2}))$ and the equation

$\displaystyle\omega(t)$
$\displaystyle=\omega_{0}-\int_{0}^{t}\Big{[}\mathrm{div}((\omega-b)\mathbf{u})+\mathrm{div}(b\mathbf{u}_{h})\Big{]}\,\,\mathrm{d}\tau$

holds for all $t\in[0,T]$.

######  Remark 2.6.

We remark that the regularity of the solution $(b,\omega,T)$ and its data
together with the integral equations immediately imply that $b$ and $\omega$
are differentiable in time. Indeed, by using the fact that
$W^{3,2}(\mathbb{T}^{2})$ and $W^{2,2}(\mathbb{T}^{2})$ are Banach algebras,
we immediately deduce from the integral equations for $b$ and $\omega$ above
that

$\displaystyle b\in C^{1}([0,T];W^{2,2}(\mathbb{T}^{2})),\qquad\omega\in
C^{1}([0,T];W^{1,2}(\mathbb{T}^{2})).$

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
$\omega(0,x)=\omega_{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,
$\int_{\mathbb{T}^{2}}\mathbf{u}_{h}\,\mathrm{d}\mathbf{x}=0$ and
$\int_{\mathbb{T}^{2}}\mathbf{u}\,\mathrm{d}\mathbf{x}=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,\omega,T_{\max})$ a maximal solution to the system (1.1)–(1.5) if:

  * •

there exists an increasing sequence of time steps $(T_{n})_{n\in\mathbb{N}}$
whose limit is $T_{\max}\in(0,\infty]$;

  * •

for each $n\in\mathbb{N}$, the triple $(b,\omega,T_{n})$ is a local strong
solution to the system (1.1)–(1.5) with initial condition
$(b_{0},\omega_{0})$;

  * •

if $T_{\max}<\infty$, then

(2.9) $\displaystyle\limsup_{T_{n}\rightarrow
T_{\max}}\|(b,\omega)(T_{n})\|_{\mathcal{M}}^{2}=\infty.$

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.

######  Theorem 2.10 (Existence of local solutions).

There exists a solution $(b,\omega,T)$ of (1.1)–(1.5) under Assumption 2.4.

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.

######  Theorem 2.11 (Continuity at low regularity).

Let $\mathbf{u}_{h}\in W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$
and $f\in W^{2,2}(\mathbb{T}^{2})$. Assume that $(b^{1},\omega^{1},T^{1})$ and
$(b^{2},\omega^{2},T^{2})$ are solutions of (1.1)–(1.5) with initial
conditions $(b_{0}^{1},\omega_{0}^{1})\in\mathcal{M}$ and
$(b_{0}^{2},\omega_{0}^{2})\in\mathcal{M}$ respectively. Then there exists a
constant

$\displaystyle
c=c\big{(}\|b_{0}^{1}\|_{3,2},\|\omega_{0}^{1}\|_{2,2},\|b_{0}^{2}\|_{3,2},\|\omega_{0}^{2}\|_{2,2},\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2}\big{)}$

such that

(2.10)
$\displaystyle\|(b^{1}-b^{2})(t)\|_{2,2}^{2}+\|(\omega^{1}-\omega^{2})(t)\|_{1,2}^{2}\leq\exp(cT)\big{(}\|b_{0}^{1}-b_{0}^{2}\|_{2,2}^{2}+\|\omega_{0}^{1}-\omega_{0}^{2}\|_{1,2}^{2}\big{)}$

holds for all $t\in[0,T]$ where $T=\min\\{T^{1},T^{2}\\}$.

######  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}(\mathbb{T}^{2})\times W^{1,2}(\mathbb{T}^{2})$ with smaller norms
than the space of existence $W^{3,2}(\mathbb{T}^{2})\times
W^{2,2}(\mathbb{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_{0}^{1}\|_{3,2},\|\omega_{0}^{1}\|_{2,2}$ rather than in the weaker space
for which we obtain our final stability estimate. Furthermore, we also require
the potential vorticity $\omega^{2}$ (at all time of existence) of the second
solution to be also bounded in the stronger space of existence
$W^{2,2}(\mathbb{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}(\mathbb{T}^{2})$. Unfortunately
however, because of the highly coupled nature of the system of equations under
study, needing $\omega^{2}\in W^{2,2}(\mathbb{T}^{2})$ for all times
automatically requires that $b^{2}\in W^{3,2}(\mathbb{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_{0}^{2}\|_{3,2},\|\omega_{0}^{2}\|_{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.

######  Corollary 2.13 (Uniqueness).

Let $(b^{1},\omega^{1},T^{1})$ and $(b^{2},\omega^{2},T^{2})$ be two solutions
of (1.1)–(1.5) under Assumption 2.4. Then the difference
$(b^{1}-b^{2},\omega^{1}-\omega^{2})$ satisfies the equation

(2.11)
$\displaystyle\|(b^{1}-b^{2})(t)\|_{2,2}^{2}+\|(\omega^{1}-\omega^{2})(t)\|_{1,2}^{2}=0$

for all $t\in[0,T]$ where $T=\min\\{T^{1},T^{2}\\}$.

Finally, we can show that a maximal solution of (1.1)–(1.5), in the sense of
Definition 2.8, also exists.

######  Theorem 2.14 (Existence of maximal solution).

There exist a unique maximal solution $(b,\omega,T_{\max})$ of (1.1)–(1.5)
under Assumption 2.4.

The $\alpha$-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,\omega)$ (and not explicitly in terms of $\mathbf{u}$), henceforth,
we will use the triple $(b^{\alpha},\omega^{\alpha},T)$ to denote the
corresponding solution to the following $\alpha$-TQG model for the avoidance
of confusion. To be precise, for fixed $\alpha>0$, we will be exploring the
pair $(b^{\alpha},w^{\alpha})$ that solves

(2.12) $\displaystyle\frac{\partial}{\partial
t}b^{\alpha}+{\mathbf{u}}^{\alpha}\cdot\nabla b^{\alpha}=0,$ (2.13)
$\displaystyle\frac{\partial}{\partial
t}\omega^{\alpha}+{\mathbf{u}}^{\alpha}\cdot\nabla(\omega^{\alpha}-b^{\alpha})=-\mathbf{u}_{h}\cdot\nabla
b^{\alpha}$

but where now,

(2.14)
$\displaystyle{\mathbf{u}}^{\alpha}=\nabla^{\perp}\psi^{\alpha},\qquad\mathbf{u}_{h}=\frac{1}{2}\nabla^{\perp}h,\qquad\omega^{\alpha}=(\Delta-1)(1-\alpha\Delta)\psi^{\alpha}+f.$

Note that at least formally, this implies that

(2.15)
$\displaystyle{\mathbf{u}}^{\alpha}=\nabla^{\perp}(1-\alpha\Delta)^{-1}(\Delta-1)^{-1}(\omega^{\alpha}-f).$

We also note that even though the velocity fields $\mathbf{u}$ in (1.5) and
${\mathbf{u}}^{\alpha}$ 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 $\alpha$-TQG
model (2.12)–(2.14) introduced above.

######  Theorem 2.15.

Fix $\alpha>0$. There exist a unique maximal solution
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ of (2.12)–(2.14) under
Assumption 2.4.

######  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 $\alpha$,
$\psi^{\alpha}\sim\tilde{\psi}\sim\psi$ and thus, both models are the same.
Nevertheless, rigorously, the key lies in the estimate (2.7) which still holds
true when one replaces $\mathbf{u}$ with ${\mathbf{u}}^{\alpha}$ and where
now, the triple $({\mathbf{u}}^{\alpha},\omega^{\alpha},f)$ satisfies

(2.16)
$\displaystyle{\mathbf{u}}^{\alpha}=\nabla^{\perp}\psi^{\alpha},\qquad(\omega^{\alpha}-f)=(\Delta-1)(1-\alpha\Delta)\psi^{\alpha}.$

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^{\alpha},\omega^{\alpha},T^{\alpha})$ of (2.12)–(2.14) for any $\alpha>0$.
Moreover, $(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ is a unique maximal
solution of (2.12)–(2.14).

###### Proof.

The existence of a unique solution $(b^{\alpha},\omega^{\alpha},T^{\alpha})$
of (2.12)–(2.14) for any $\alpha>0$ follow from Proposition 3.1, Lemma 3.5 and
Corollary 2.13 with $T^{\alpha}=T$. We will prove all these results later for
the TQG and one will observe that nothing changes for the $\alpha$-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 $\alpha>0$.
The construction of the maximal solution
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ of (2.12)–(2.14) also follow
the same gluing argument used in showing the existence of the maximal solution
$(b,\omega,T_{\max})$ for (1.1)–(1.5) in Section 3.4. ∎

##  3\. Construction of strong solutions

In this section, we give a construction of a solution, in the sense of
Definition 2.5, to the thermal quasi-geostrophic 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,\omega)$ satisfies the assumptions of Theorem 8.1 in the
appendix. The details are as follows. We consider

(3.1) $\displaystyle\frac{\partial}{\partial
t}\binom{b}{\omega}+\mathcal{A}\begin{pmatrix}b\\\
\omega\end{pmatrix}=0,\qquad t\geq 0,\qquad
b\big{|}_{t=0}=b_{0},\qquad\omega\big{|}_{t=0}=\omega_{0}$

where the operator

(3.2)
$\displaystyle\mathcal{A}:=\begin{bmatrix}\mathbf{u}\cdot\nabla&0\\\\[3.00003pt]
(\mathbf{u}_{h}-\mathbf{u})\cdot\nabla&\mathbf{u}\cdot\nabla\end{bmatrix}$

is defined such that

(3.3)
$\displaystyle\mathbf{u}=\nabla^{\perp}\psi,\qquad\mathbf{u}_{h}=\frac{1}{2}\nabla^{\perp}h,\qquad\omega=(\Delta-1)\psi+f.$

###  3.1. Estimates for convective terms

In the following, we let $k\in\mathbb{N}\cup\\{0\\}$ and recall that
$W^{k,p}(\mathbb{T}^{2})$ is a Banach algebra when $kp>2$. By using Sobolev
embeddings and (2.7), we obtain the following estimates.

  1. (1b)

Let $k\in\\{0,1,2\\}$. If $b\in W^{k+1,2}(\mathbb{T}^{2})$, $\omega\in
W^{2,2}(\mathbb{T}^{2})$ and $\mathbf{u}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ solves (1.5) for a
given $f\in W^{2,2}(\mathbb{T}^{2})$ and a given $\mathbf{u}_{h}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, then

$\displaystyle\|\mathbf{u}\cdot\nabla
b\|_{k,2}\lesssim\|\mathbf{u}\|_{2,2}\|b\|_{k+1,2}\lesssim\big{(}1+\|\omega\|_{1,2}\big{)}\|b\|_{k+1,2}$

and

$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,b\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}\|_{3,2}\|b\|_{k,2}^{2}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|b\|_{k,2}^{2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}\|_{2,2}\|b\|_{k+1,2}\|\omega\|_{k,2}\lesssim\big{(}1+\|\omega\|_{1,2}\big{)}\|b\|_{k+1,2}\|\omega\|_{k,2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}_{h}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}_{h}\|_{2,2}\|b\|_{k+1,2}\|\omega\|_{k,2}\lesssim\|b\|_{k+1,2}\|\omega\|_{k,2}.$

  2. (1b)

Let $k\geq 3$. If $b\in W^{k+1,2}(\mathbb{T}^{2})$, $\omega\in
W^{k,2}(\mathbb{T}^{2})$ and $\mathbf{u}\in
W^{k+1,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ solves (1.5) for a
given $f\in W^{k,2}(\mathbb{T}^{2})$ and a given $\mathbf{u}_{h}\in
W^{k+1,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, then

$\displaystyle\|\mathbf{u}\cdot\nabla
b\|_{k,2}\lesssim\|\mathbf{u}\|_{2,2}\|b\|_{k+1,2}+\|\mathbf{u}\|_{k,2}\|b\|_{3,2}\lesssim\big{(}1+\|\omega\|_{k-1,2}\big{)}\|b\|_{k+1,2}$

and

$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,b\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}\|_{3,2}\|b\|_{k,2}^{2}+\|\mathbf{u}\|_{k,2}\|b\|_{k,2}\|b\|_{3,2}\lesssim\big{(}1+\|\omega\|_{k-1,2}\big{)}\|b\|_{k,2}^{2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\big{(}\|\mathbf{u}\|_{2,2}\|b\|_{k+1,2}+\|\mathbf{u}\|_{k,2}\|b\|_{3,2}\big{)}\|\omega\|_{k,2}\lesssim\big{(}1+\|\omega\|_{k-1,2}\big{)}\|b\|_{k+1,2}\|\omega\|_{k,2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}_{h}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\big{(}\|\mathbf{u}_{h}\|_{2,2}\|b\|_{k+1,2}+\|\mathbf{u}_{h}\|_{k,2}\|b\|_{3,2}\big{)}\|\omega\|_{k,2}\lesssim\|b\|_{k+1,2}\|\omega\|_{k,2}.$

  3. (2a)

Let $k\in\\{0,1,2\\}$. If $\omega\in W^{k+1,2}(\mathbb{T}^{2})$ and that
$\mathbf{u}\in W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ solves
(1.5) for a given $f\in W^{2,2}(\mathbb{T}^{2})$, then

$\displaystyle\|\mathbf{u}\cdot\nabla\omega\|_{k,2}\lesssim\|\mathbf{u}\|_{2,2}\|\omega\|_{k+1,2}\lesssim\big{(}1+\|\omega\|_{1,2}\big{)}\|\omega\|_{k+1,2}$

and

$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla\omega\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}\|_{3,2}\|\omega\|_{k,2}^{2}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|\omega\|_{k,2}^{2}.$

  4. (2b)

Let $k\geq 3$. If $\omega\in W^{k+1,2}(\mathbb{T}^{2})$ and $\mathbf{u}\in
W^{k+1,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ solves (1.5) for a
given $f\in W^{k,2}(\mathbb{T}^{2})$, then

$\displaystyle\|\mathbf{u}\cdot\nabla\omega\|_{k,2}\lesssim\|\mathbf{u}\|_{2,2}\|\omega\|_{k+1,2}+\|\mathbf{u}\|_{k,2}\|\omega\|_{3,2}\lesssim\big{(}1+\|\omega\|_{k-1,2}\big{)}\|\omega\|_{k+1,2}$

and

$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla\omega\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\|\mathbf{u}\|_{3,2}\|\omega\|_{k,2}^{2}+\|\mathbf{u}\|_{k,2}\|\omega\|_{k,2}\|\omega\|_{3,2}\lesssim\big{(}1+\|\omega\|_{k-1,2}\big{)}\|\omega\|_{k,2}^{2}.$

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,\omega,T)$ of (3.1) such
that:

  1. (1)

the time $T>0$ satisfy the bound

(3.4) $\displaystyle
T<\frac{1}{2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})},\quad
c=c(\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2});$

  2. (2)

the pair $(b,\omega)$ is of class

(3.5) $\displaystyle b\in C_{w}([0,T];W^{3,2}(\mathbb{T}^{2})),\qquad\omega\in
C_{w}([0,T];W^{2,2}(\mathbb{T}^{2}));$

  3. (3)

the pair $(b,\omega)$ satisfies the bound

(3.6)
$\displaystyle\|(b,\omega)(t)\|_{\mathcal{M}}^{2}\leq\frac{\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}+2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})\,t}{1-2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})\,t},\quad\quad
t\in(0,T).$

###### Proof.

Let us consider the following triplet $\\{V,H,X\\}$ given as follows. Take
$X:=W^{1,2}(\mathbb{T}^{2})\times L^{2}(\mathbb{T}^{2})$ and let
$H:=W^{(3,2),2}(\mathbb{T}^{2})\times W^{(2,1),2}(\mathbb{T}^{2})$ be a
Hilbert space endowed with the inner product

(3.7)
$\displaystyle\big{\langle}(b,\omega),(b^{\prime},\omega^{\prime})\big{\rangle}_{H}=\langle
b\,,\,b^{\prime}\rangle_{3,2}+\langle
b\,,\,b^{\prime}\rangle_{2,2}+\langle\omega\,,\,\omega^{\prime}\rangle_{2,2}+\langle\omega\,,\,\omega^{\prime}\rangle_{1,2}.$

We note that for $(b,\omega)\in\mathcal{M}$, we have that

(3.8) $\displaystyle\|(b,\omega)\|_{\mathcal{M}}\leq\|(b,\omega)\|_{H}\leq
2\|(b,\omega)\|_{\mathcal{M}}$

and thus, $H$ is equivalent to $\mathcal{M}$. Also, by virtue of the estimates
shown in Section 3.1, we can conclude that the operator $\mathcal{A}$ is
weakly continuous from $H$ into $X$.
Next, we let $V$ be the domain of an unbounded selfadjoint operator $S\geq 0$
in $X$ with $\mathrm{domain}(S)\subset H$ and $\langle Su,v\rangle=\langle
u,v\rangle_{H}$ for $u\in\mathrm{domain}(S)$. More precisely,
$S=\sum_{|\beta|\leq 3}(-\partial)^{\beta}\partial^{\beta}$ subject to
periodic boundary conditions. With this definition of $V$ in hand, we can
again conclude from the estimates in Section 3.1 that $\mathcal{A}$ maps $V$
into $H$. Our goal now is to show that for any $(b,\omega)\in V$, there exists
an increasing function $\rho(r)\geq 0$ of $r\geq 0$ such that

(3.9)
$\displaystyle\Big{|}\big{\langle}\mathcal{A}(b,\omega),(b,\omega)\big{\rangle}_{H}\Big{|}\leq\rho\big{(}\|(b,\omega)\|^{2}_{H}\big{)}.$

To achieve this goal, we further refine the following estimates from Section
3.1. In particular, we recall that if $b\in W^{k+1,2}(\mathbb{T}^{2})$ and
$\omega\in W^{2,2}(\mathbb{T}^{2})$ and that $\mathbf{u}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ solves (1.5) for a
given $f\in W^{2,2}(\mathbb{T}^{2})$ and a given $\mathbf{u}_{h}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, then for
$k\in\\{0,1,2\\}$,

(3.10) $\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,b\big{\rangle}_{k,2}\Big{|}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|b\|_{k+1,2}^{2}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}\|b\|_{3,2}^{2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla\omega\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|\omega\|_{2,2}^{2}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}\|\omega\|_{2,2}^{2},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|b\|_{k+1,2}\|\omega\|_{2,2}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}\big{(}\|b\|_{3,2}^{2}+\|\omega\|_{2,2}^{2}\big{)},$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}_{h}\cdot\nabla
b\,,\,\omega\big{\rangle}_{k,2}\Big{|}\lesssim\|b\|_{k+1,2}\|\omega\|_{2,2}\lesssim\big{(}\|b\|_{3,2}^{2}+\|\omega\|_{2,2}^{2}\big{)}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}.$

In the above estimates, we have used inequalities such as $x\lesssim 1+x^{2}$
for $x\geq 0$ and Young’s inequalities. Furthermore, if $b\in
W^{3,2}(\mathbb{T}^{2})$, $\omega\in W^{2,2}(\mathbb{T}^{2})$ and
$\mathbf{u}\in W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, we also
have that

(3.11) $\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla
b\,,\,b\big{\rangle}_{3,2}\Big{|}\lesssim\big{(}1+\|\omega\|_{2,2}\big{)}\|b\|_{3,2}^{2}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}\|b\|_{3,2}^{2}.$

We can therefore conclude that for $(b,\omega)\in H$ (and in particular, for
any $(b,\omega)\in V$),

(3.12)
$\displaystyle\Big{|}\big{\langle}\mathcal{A}(b,\omega),(b,\omega)\big{\rangle}_{H}\Big{|}\lesssim\big{(}1+\|(b,\omega)\|^{2}_{H}\big{)}^{2}$

with a constant depending only on $\|\mathbf{u}_{h}\|_{3,2}$ and
$\|f\|_{2,2}$. Since $\rho(r)=c(1+r)^{2}\geq 0$ is an increasing function of
$r\geq 0$, we can conclude from the result in the Appendix, Section 8, that
for a given $(b_{0},\omega_{0})\in H\equiv\mathcal{M}$, there exist a solution
$(b,\omega,T)$ of (3.1) of class

(3.13) $\displaystyle b\in
C_{w}([0,T];W^{3,2}(\mathbb{T}^{2})),\qquad\omega\in
C_{w}([0,T];W^{2,2}(\mathbb{T}^{2}))$

satisfying the bound (recall (3.8))

(3.14)
$\displaystyle\|(b,\omega)(t)\|_{\mathcal{M}}^{2}\leq\|(b,\omega)(t)\|_{H}^{2}\leq
r(t),\quad\quad t\in(0,T)$

where $r$ is an increasing function on $(0,T)$. The time $T>0$ in
(3.13)–(3.14) depends only on $\|\mathbf{u}_{h}\|_{3,2}$ and $\|f\|_{2,2}$ by
way of the function $\rho(\cdot)$ as well as on
$\|(b_{0},\omega_{0})\|_{\mathcal{M}}$. To be precise, for
$\rho(r)=c(1+r)^{2}\geq 0$ where $c>0$ depends only on
$\|\mathbf{u}_{h}\|_{3,2}$ and $\|f\|_{2,2}$, we obtain $T>0$ by solving the
equation

(3.15) $\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}r=2\rho(r),\qquad\qquad
r(0)=\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}$

which yields

(3.16) $\displaystyle
r(t)=\frac{\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}+2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})\,t}{1-2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})\,t}.$

We therefore require

(3.17) $\displaystyle
T<\frac{1}{2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})}$

for a solution of (3.15) to exist. This finishes the proof. ∎

######  Remark 3.2.

Note that $r(t)$ converges to $\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}$ as
$t\rightarrow 0^{+}$. Therefore, given the bound (3.14) and (3.16), we can
conclude that

$\displaystyle\limsup_{t\rightarrow
0^{+}}\|(b,\omega)(t)\|^{2}_{\mathcal{M}}\leq\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}.$

On the other hand, because the pair $(b,\omega)$ is of class (3.13), we can
also conclude that

$\displaystyle\liminf_{t\rightarrow
0^{+}}\|(b,\omega)(t)\|^{2}_{\mathcal{M}}\geq\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}$

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.3.

By time-reversible, we mean that if $b(t,x)$ and $\omega(t,x)$ (with
$\omega(t,x)$ related to $\mathbf{u}(t,x)$ by (1.5)) solves (1.1)–(1.5) with
data $b_{0}(x)$, $\omega_{0}(x)$, $\mathbf{u}_{h}(x)$ and $f(x)$, then
$-b(-t,x)$ and $-\omega(-t,x)$ (with $-\omega(-t,x)$ related to
$-\mathbf{u}(-t,x)$ by (1.5)) solves (1.1)–(1.5) with data $-b_{0}(x)$,
$-\omega_{0}(x)$, $-\mathbf{u}_{h}(x)$ and $-f(x)$ respectively.

######  Remark 3.4.

We remark that the function $\rho$ 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 (LABEL:convEst1)–(LABEL:convEst2).
For example, a different choice of Hölder conjugates for estimating the
aforementioned terms will lead to a different $\rho$ and $r$. The importance
of this remark lies in the fact that $\rho$ and $r$ determines the longevity
of the solution $(b,\omega,T)$ to (1.1)–(1.5). One can therefore tune $T>0$ by
modifying $\rho$ and $r$ accordingly. For example, we may also obtain the
following

(3.18)
$\displaystyle\rho_{1}(r)=c(r^{2}+r),\qquad\quad\rho_{2}(r)=c(r^{\frac{3}{2}}+r)$

for which the corresponding solutions to (3.15) are

(3.19) $\displaystyle r_{1}(t)=\frac{r(0)e^{2ct}}{r(0)+1-r(0)e^{2ct}},$ (3.20)
$\displaystyle r_{2}(t)=\frac{s^{2}+s\pm
2s\sqrt{s}}{s^{2}-2s+1},\quad\text{where}\quad
s=s(t)=\frac{r(0)}{(\sqrt{r(0)}+1)^{2}}e^{2ct}.$

respectively. Again, the constants depends only on $\|\mathbf{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

$\displaystyle
T_{1}<\frac{1}{2c}\ln\bigg{(}1+\frac{1}{\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}}\bigg{)}$

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,\omega,T)$ to (1.1)–(1.5) as compared to $r_{1}$
since

$\displaystyle\frac{1}{2c}\ln\bigg{(}1+\frac{1}{\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}}\bigg{)}\leq\frac{1}{2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})}.$

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},\omega^{1},T)$ and $(b^{2},\omega^{2},T)$ be
two solutions of (1.1)–(1.5) sharing the same data $\mathbf{u}_{h}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$, $f\in
W^{2,2}(\mathbb{T}^{2})$ and $(b_{0},\omega_{0})$. So in particular,

(3.21)
$\displaystyle\sup_{t\in[0,T)}\|b^{i}(t,\cdot)\|_{3,2}<\infty,\qquad\sup_{t\in[0,T)}\|\omega^{i}(t,\cdot)\|_{2,2}<\infty$

hold for each $i=1,2$. We now set $b^{12}:=b^{1}-b^{2}$,
$\omega^{12}:=\omega^{1}-\omega^{2}$ and
$\mathbf{u}^{12}:=\mathbf{u}^{1}-\mathbf{u}^{2}$ so that
$(b^{12},\mathbf{u}^{12},\omega^{12})$ satisfies

(3.22) $\displaystyle\frac{\partial}{\partial
t}b^{12}+\mathbf{u}^{12}\cdot\nabla b^{1}+\mathbf{u}^{2}\cdot\nabla b^{12}=0,$
(3.23) $\displaystyle\frac{\partial}{\partial
t}\omega^{12}+\mathbf{u}^{12}\cdot\nabla(\omega^{1}-b^{1})+\mathbf{u}^{2}\cdot\nabla(\omega^{12}-b^{12})=-\mathbf{u}_{h}\cdot\nabla
b^{12},$

where

(3.24)
$\displaystyle\mathbf{u}^{12}=\nabla^{\perp}\psi^{12},\qquad\mathbf{u}_{h}=\frac{1}{2}\nabla^{\perp}h,\qquad\omega^{12}=(\Delta-1)\psi^{12}.$

We now show the following stability and uniqueness results in the space
$W^{2,2}(\mathbb{T}^{2})\times W^{1,2}(\mathbb{T}^{2})$ which is larger than
the space $W^{3,2}(\mathbb{T}^{2})\times W^{2,2}(\mathbb{T}^{2})$ of existence
of the individual solutions. Indeed, it suffices to show uniqueness in the
even larger space $W^{1,2}(\mathbb{T}^{2})\times L^{2}(\mathbb{T}^{2})$ but we
avoid doing so since we wish to concurrently show the proof of Theorem 2.11 as
well.

###### Proof of Theorem 2.11 and Corollary 2.13.

If we apply $\partial^{\beta}$ to (3.22) with $|\beta|\leq 2$, we obtain

(3.25) $\displaystyle\frac{\partial}{\partial
t}\partial^{\beta}b^{12}+\mathbf{u}^{2}\cdot\nabla\partial^{\beta}b^{12}=S_{1}+S_{2}-\mathbf{u}^{12}\cdot\nabla\partial^{\beta}b^{1},$

where

$\displaystyle S_{1}$
$\displaystyle:=\mathbf{u}^{12}\cdot\partial^{\beta}\nabla
b^{1}-\partial^{\beta}(\mathbf{u}^{12}\cdot\nabla b^{1}),$ $\displaystyle
S_{2}$ $\displaystyle:=\mathbf{u}^{2}\cdot\partial^{\beta}\nabla
b^{12}-\partial^{\beta}(\mathbf{u}^{2}\cdot\nabla b^{12})$

are such that

(3.26) $\displaystyle\|S_{1}\|_{2}$
$\displaystyle\lesssim\|\nabla\mathbf{u}^{12}\|_{4}\|b^{1}\|_{2,4}+\|\nabla
b^{1}\|_{\infty}\|\mathbf{u}^{12}\|_{2,2}\lesssim\|\mathbf{u}^{12}\|_{2,2}\|b^{1}\|_{3,2},$
(3.27) $\displaystyle\|S_{2}\|_{2}$
$\displaystyle\lesssim\|\nabla\mathbf{u}^{2}\|_{\infty}\|b^{12}\|_{2,2}+\|\nabla
b^{12}\|_{4}\|\mathbf{u}^{2}\|_{2,4}\lesssim\|\mathbf{u}^{2}\|_{3,2}\|b^{12}\|_{2,2}.$

Additionally, the following estimate for the $L^{2}$ inner product of the last
term (3.25) with $\partial^{\beta}b^{12}$ holds

(3.28)
$\displaystyle\Big{|}\big{\langle}\mathbf{u}^{12}\cdot\nabla\partial^{\beta}b^{1}\,,\,\partial^{\beta}b^{12}\big{\rangle}\Big{|}\lesssim\|b^{1}\|_{3,2}\|\mathbf{u}^{12}\|_{\infty}\|b^{12}\|_{2,2}.$

If we test (3.25) with $\partial^{\beta}b^{12}$ and sum over $|\beta|\leq 2$,
we obtain from the above estimates together with (2.7) and (3.21),

(3.29) $\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|b^{12}\|_{2,2}^{2}$
$\displaystyle\lesssim\|\omega^{12}\|_{1,2}^{2}+\|b^{12}\|_{2,2}^{2}.$

Next, we apply $\partial^{\beta}$ to (3.23) with $|\beta|\leq 1$ to get

(3.30) $\displaystyle\frac{\partial}{\partial
t}\partial^{\beta}\omega^{12}+\mathbf{u}^{2}\cdot\nabla\partial^{\beta}\omega^{12}=$
$\displaystyle-\mathbf{u}^{12}\cdot\nabla\partial^{\beta}(\omega^{1}-b^{1})+(\mathbf{u}^{2}-\mathbf{u}_{h})\cdot\nabla\partial^{\beta}b^{12}$
$\displaystyle+S_{3}+\ldots+S_{7}$

where

$\displaystyle S_{3}$
$\displaystyle:=\mathbf{u}^{2}\cdot\partial^{\beta}\nabla\omega^{12}-\partial^{\beta}(\mathbf{u}^{2}\cdot\nabla\omega^{12}),$
$\displaystyle S_{4}$
$\displaystyle:=\mathbf{u}^{12}\cdot\partial^{\beta}\nabla\omega^{1}-\partial^{\beta}(\mathbf{u}^{12}\cdot\nabla\omega^{1}),$
$\displaystyle S_{5}$
$\displaystyle:=-\mathbf{u}^{12}\cdot\partial^{\beta}\nabla
b^{1}+\partial^{\beta}(\mathbf{u}^{12}\cdot\nabla b^{1}),$ $\displaystyle
S_{6}$ $\displaystyle:=-\mathbf{u}^{2}\cdot\partial^{\beta}\nabla
b^{12}+\partial^{\beta}(\mathbf{u}^{2}\cdot\nabla b^{12}),$ $\displaystyle
S_{7}$ $\displaystyle:=\mathbf{u}_{h}\cdot\partial^{\beta}\nabla
b^{12}-\partial^{\beta}(\mathbf{u}_{h}\cdot\nabla b^{12})$

are such that

(3.31) $\displaystyle\|S_{3}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{2}\|_{1,\infty}\|\omega^{12}\|_{1,2}\lesssim\|\omega^{12}\|_{1,2},$
(3.32) $\displaystyle\|S_{4}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{1,4}\|\omega^{1}\|_{1,4}\lesssim\|\mathbf{u}^{12}\|_{2,2}\lesssim\|\omega^{12}\|_{1,2},$
(3.33) $\displaystyle\|S_{5}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{1,4}\|b^{1}\|_{1,4}\lesssim\|\mathbf{u}^{12}\|_{2,2}\lesssim\|\omega^{12}\|_{1,2},$
(3.34) $\displaystyle\|S_{6}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{2}\|_{1,\infty}\|b^{12}\|_{1,2}\lesssim\|b^{12}\|_{1,2}\lesssim\|b^{12}\|_{2,2},$
(3.35) $\displaystyle\|S_{6}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}_{h}\|_{1,\infty}\|b^{12}\|_{1,2}\lesssim\|b^{12}\|_{1,2}\lesssim\|b^{12}\|_{2,2}.$

Also, since $\mathrm{div}(\mathbf{u}^{2})=0$

(3.36)
$\displaystyle\Big{|}\big{\langle}\mathbf{u}^{2}\cdot\nabla\partial^{\beta}\omega^{12}\,,\,\partial^{\beta}\omega^{12}\big{\rangle}\Big{|}=0,$
$\displaystyle\Big{|}\big{\langle}\mathbf{u}^{12}\cdot\nabla\partial^{\beta}(\omega^{1}-b^{1})\,,\,\partial^{\beta}\omega^{12}\big{\rangle}\Big{|}\lesssim\|\omega^{1}-b^{1}\|_{2,2}\|\mathbf{u}^{12}\|_{\infty}\|\omega^{12}\|_{1,2},$
$\displaystyle\Big{|}\big{\langle}(\mathbf{u}^{2}-\mathbf{u}_{h})\cdot\nabla\partial^{\beta}b^{12}\,,\,\partial^{\beta}\omega^{12}\big{\rangle}\Big{|}\lesssim\|\mathbf{u}^{2}-\mathbf{u}_{h}\|_{\infty}\|b^{12}\|_{2,2}\|\omega^{12}\|_{1,2}.$

Therefore,

(3.37)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\omega^{12}\|_{1,2}^{2}$
$\displaystyle\lesssim\ \|\omega^{12}\|_{1,2}^{2}+\|b^{12}\|_{2,2}^{2}.$

From (3.29) and (3.37), we can conclude from Grönwall’s lemma that

(3.38)
$\displaystyle\|b^{12}(t,\cdot)\|_{2,2}^{2}+\|\omega^{12}(t,\cdot)\|_{1,2}^{2}\leq\exp(ct)\big{(}\|b^{12}_{0}\|_{2,2}^{2}+\|\omega^{12}_{0}\|_{1,2}^{2}\big{)}$

for all $t\geq 0$ where the constant $c>0$ depends on
$(\|b^{1}\|_{3,2},\|\omega^{1}\|_{2,2},\|\omega^{2}\|_{2,2},\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2})$.
This finishes the proof of Theorem 2.11. To obtain uniqueness, i.e. Corollary
2.13, we use the fact that $b^{12}_{0}=0$ and $\omega^{12}_{0}=0$. ∎

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

$\displaystyle b\in C([0,T];W^{3,2}(\mathbb{T}^{2})),\qquad\omega\in
C([0,T];W^{2,2}(\mathbb{T}^{2})).$

###### 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,\omega,T_{0})$ of (1.1)–(1.2) where $T_{0}\in[0,T]$ is fixed but
arbitrary. Then by (3.6), this solution will satisfy the bound

(3.39)
$\displaystyle\|(b,\omega)(T_{0})\|^{2}_{\mathcal{M}}\leq\frac{\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}}+2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})T_{0}}{1-2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})T_{0}}$

$c=c(\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2})$. We can now construct a new
solution $(b^{1},\omega^{1},T_{0}+T_{1})$ by taking $(b,\omega)(T_{0})$ as an
initial condition. A repetition of the argument leading to (3.39) means that
this new pair $(b^{1},\omega^{1})$ solving (1.1)–(1.2) will satisfy the
inequality

(3.40)
$\displaystyle\|(b^{1},\omega^{1})(t)\|^{2}_{\mathcal{M}}\leq\frac{\|(b,\omega)(T_{0})\|^{2}_{\mathcal{M}}+2c(1+\|(b,\omega)(T_{0})\|^{2}_{\mathcal{M}})t}{1-2c(1+\|(b,\omega)(T_{0})\|^{2}_{\mathcal{M}})t}$

on $[T_{0},T_{0}+T_{1}]$ where

$\displaystyle T_{1}<\frac{1}{2c(1+\|(b,\omega)(T_{0})\|^{2}_{\mathcal{M}})}.$

Furthermore, $(b^{1},\omega^{1})$ must coincide with $(b,\omega)$ on
$[T_{0},T_{0}+T_{1}]\cap[0,T]$ by uniqueness and the fact that they agree at
$T_{0}\in[0,T]$. Subsequently, we obtain from (3.40), strong right-continuity
of $\|(b^{1},\omega^{1})\|_{\mathcal{M}}$ at time $t=T_{0}$ just as was done
for $t=0$ in Remark 3.2. Again, by uniqueness, this implies that
$\|(b,\omega)\|_{\mathcal{M}}$ is also strongly right-continuous at time
$t=T_{0}$. Since $T_{0}$ was chosen arbitrarily, this means that $(b,\omega)$
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

(3.41) $\displaystyle
T<\frac{1}{2c(1+\|(b_{0},\omega_{0})\|^{2}_{\mathcal{M}})}$

with $c=c(\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2})$ such that $(b,\omega,T)$ is a
unique solution of (1.1)–(1.5). At time $T>0$ given above, we now choose
$(b_{T},\omega_{T})\in\mathcal{M}$, $\mathbf{u}_{h}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ and $f\in
W^{2,2}(\mathbb{T}^{2})$ as our new data where

$\displaystyle b_{T}:=b(T,\cdot),\qquad\omega_{T}:=\omega(T,\cdot).$

By repeating the argument, we can also find

(3.42)
$\displaystyle\tilde{T}<\frac{1}{2c(1+\|(b_{T},\omega_{T})\|^{2}_{\mathcal{M}})}$

with $c=c(\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2})$, such that with the new
initial conditions $(b_{T},\omega_{T})$, $(b,\omega)$ uniquely solves
(1.1)–(1.5), in the sense of Definition 2.5, on $[T,T+\tilde{T}]$. By a gluing
argument, we obtain a solution $(b,\omega,T_{1})$ of (1.1)–(1.5) with
$T_{1}=T+\tilde{T}$. By this iterative procedure:

  * •

we obtain an increasing family $(T_{n})){n\in\mathbb{N}}$ of time steps whose
limit is $T_{\max}\in(0,\infty]$;

  * •

for each $n\in\mathbb{N}$, $(b,\omega,T_{n})$ is a solution of (1.1)–(1.5);

  * •

if $T_{\max}<\infty$, then

$\displaystyle\limsup_{T_{n}\rightarrow
T_{\max}}\|(b,\omega)(T_{n})\|_{\mathcal{M}}^{2}=\infty$

since otherwise, we can repeat the procedure above to obtain
$\tilde{\tilde{T}}>0$ satisfying

$\displaystyle\tilde{\tilde{T}}<\frac{1}{2c(1+\|(b,\omega)(T_{\max})\|^{2}_{\mathcal{M}})}$

with $c=c(\|\mathbf{u}_{h}\|_{3,2},\|f\|_{2,2})$, such that
$(b,\omega,T_{\max}^{\max})$ is a solution of (1.1)–(1.5) with
$T_{\max}^{\max}=T_{\max}+\tilde{\tilde{T}}$. This will contradict the fact
that $T_{\max}$ is the maximal time. ∎

##  4\. Convergence of $\alpha$-TQG to TQG

In this section, for a given $\mathbf{u}_{h}\in
W^{3,2}_{\mathrm{div}}(\mathbb{T}^{2};\mathbb{R}^{2})$ and $f\in
W^{2,2}(\mathbb{T}^{2})$ serving as data for both the TQG (1.1)–(1.5) and the
$\alpha$-TQG (2.12)–(2.14), we aim to show that any family
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})_{\alpha>0}$ of unique maximal
solutions to the $\alpha$-TQG model (2.12)–(2.14) converges strongly to the
unique maximal solution $(b,\omega,T_{\max})$ of the TQG model (1.1)–(1.5) as
$\alpha\rightarrow 0$ provided that both system share the same initial
conditions and are defined on a common time interval $[0,T]$. Our main result
is the following.

######  Proposition 4.1.

Under Assumption 2.4, let
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})_{\alpha>0}$ be a family of
unique maximal solutions of (2.12)–(2.14) and let $(b,\omega,T_{\max})$ be the
unique maximal solution of (1.1)–(1.5). Then

(4.1) $\displaystyle\sup_{t\in[0,T]}\big{(}\|({b}^{\alpha}-b)(t)\|_{1,2}^{2}$
$\displaystyle+\|({\omega}^{\alpha}-\omega)(t)\|_{2}^{2}\big{)}\lesssim\alpha^{2}T\exp(cT)\big{[}1+\exp(cT)\big{]}$

for any $T<\min\\{T_{\max},T_{\max}^{\alpha}\\}$ and

(4.2) $\displaystyle\sup_{t\in[0,T]}\|({b}^{\alpha}-b)(t)\|_{1,2}\rightarrow
0,\quad\qquad\sup_{t\in[0,T]}\|({\omega}^{\alpha}-\omega)(t)\|_{2}\rightarrow
0$

as $\alpha\rightarrow 0$.

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
$(\overline{b}^{\alpha},\overline{\omega}^{\alpha},\overline{T}^{\alpha}_{\max})$
with initial conditions $b_{0}\in W^{3,2}(\mathbb{T}^{2})$ and $\omega_{0}\in
W^{2,2}(\mathbb{T}^{2})$ of the following “intermediate $\alpha$-TQG” equation
given by

(4.3) $\displaystyle\frac{\partial}{\partial
t}\overline{b}^{\alpha}+\overline{\mathbf{u}}^{\alpha}\cdot\nabla\overline{b}^{\alpha}=0,$
(4.4) $\displaystyle\frac{\partial}{\partial
t}\overline{\omega}^{\alpha}+\overline{\mathbf{u}}^{\alpha}\cdot\nabla(\overline{\omega}^{\alpha}-\overline{b}^{\alpha})=-\mathbf{u}_{h}\cdot\nabla\overline{b}^{\alpha},$

where like (2.15),

(4.5)
$\displaystyle\overline{\mathbf{u}}^{\alpha}=\nabla^{\perp}(1-\alpha\Delta)^{-1}(\Delta-1)^{-1}(\omega-f),\quad\mathbf{u}_{h}=\frac{1}{2}\nabla^{\perp}h,$

and where $\omega$(and not $\omega^{\alpha}$) 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
$\bf{i}\alpha$-TQG for short, follow from the proof of existence for the
original TQG and it particular, for a uniform-in-$\alpha$ set of data
$(\mathbf{u}_{h},f,b_{0},\omega_{0})$, we have the following uniform-
in-$\alpha$ estimates

(4.6)
$\displaystyle\sup_{\alpha>0}\sup_{t\in[0,\overline{T}^{\alpha}_{\max})}\|\overline{b}^{\alpha}\|_{3,2}\lesssim
1,\quad\sup_{\alpha>0}\sup_{t\in[0,\overline{T}^{\alpha}_{\max})}\|\overline{\omega}^{\alpha}\|_{2,2}\lesssim
1$

for the maximal solution of (4.3)–(4.5). Indeed, exactly as in Theorem 2.15,
we also have the following result.

######  Lemma 4.2.

Fix $\alpha>0$. There exist a unique maximal solution
$(\overline{b}^{\alpha},\overline{\omega}^{\alpha},\overline{T}^{\alpha}_{\max})$
of (4.3)–(4.5) under Assumption 2.4. In particular, uniformly of $\alpha>0$,
the inequality

(4.7)
$\displaystyle\|(\overline{b}^{\alpha},\overline{\omega}^{\alpha})(t)\|_{\mathcal{M}}\lesssim
1$

holds for $t<\overline{T}^{\alpha}_{\max}$.

As a next step, we show that any family
$(\overline{b}^{\alpha},\overline{\omega}^{\alpha},\overline{T}^{\alpha}_{\max})_{\alpha>0}$
of maximal solutions to the intermediate $\mathbf{i}\alpha$-TQG model
(4.3)–(4.5) (rather than of the $\alpha$-TQG model (2.12)–(2.14)) converges
strongly to the unique maximal solution of the TQG model (1.1)–(1.5) as
$\alpha\rightarrow 0$ on the time interval $[0,T]$ where
$T<\min\\{\overline{T}^{\alpha}_{\max},T_{\max}\\}$. In order to achieve this
goal, we replicate the uniqueness argument in Section 3.3 by setting
$b^{12}:=\overline{b}^{\alpha}-b$,
$\omega^{12}:=\overline{\omega}^{\alpha}-\omega$ and
$\mathbf{u}^{12}:=\overline{\mathbf{u}}^{\alpha}-\mathbf{u}$ so that
$(b^{12},\mathbf{u}^{12},\omega^{12})$ satisfies

(4.8) $\displaystyle\frac{\partial}{\partial
t}b^{12}+\mathbf{u}^{12}\cdot\nabla\overline{b}^{\alpha}+\mathbf{u}\cdot\nabla
b^{12}=0,$ (4.9) $\displaystyle\frac{\partial}{\partial
t}\omega^{12}+\mathbf{u}^{12}\cdot\nabla(\overline{\omega}^{\alpha}-\overline{b}^{\alpha})+\mathbf{u}\cdot\nabla(\omega^{12}-b^{12})=-\mathbf{u}_{h}\cdot\nabla
b^{12}.$

Similar to the uniqueness argument, if we apply $\partial^{\beta}$ to the
equation for $b^{12}$ above where now $|\beta|\leq 1$, we obtain

(4.10) $\displaystyle\frac{\partial}{\partial
t}\partial^{\beta}b^{12}+\mathbf{u}\cdot\nabla\partial^{\beta}b^{12}=S_{1}+S_{2}-\mathbf{u}^{12}\cdot\nabla\partial^{\beta}\overline{b}^{\alpha},$

where

$\displaystyle S_{1}$
$\displaystyle:=\mathbf{u}^{12}\cdot\partial^{\beta}\nabla\overline{b}^{\alpha}-\partial^{\beta}(\mathbf{u}^{12}\cdot\nabla\overline{b}^{\alpha})=-\partial\mathbf{u}^{12}\cdot\nabla\overline{b}^{\alpha},$
$\displaystyle S_{2}$ $\displaystyle:=\mathbf{u}\cdot\partial^{\beta}\nabla
b^{12}-\partial^{\beta}(\mathbf{u}\cdot\nabla
b^{12})=-\partial\mathbf{u}\cdot\nabla b^{12}$

are such that

(4.11)
$\displaystyle\|S_{1}\|_{2}\lesssim\|\mathbf{u}^{12}\|_{1,2}\|\overline{b}^{\alpha}\|_{3,2},\qquad\qquad\|S_{2}\|_{2}\lesssim\|\mathbf{u}\|_{3,2}\|b^{12}\|_{1,2}.$

Also,

(4.12)
$\displaystyle\Big{|}\big{\langle}\mathbf{u}^{12}\cdot\nabla\partial^{\beta}\overline{b}^{\alpha}\,,\,\partial^{\beta}b^{12}\big{\rangle}\Big{|}\lesssim\|\mathbf{u}^{12}\|_{1,2}\|\overline{b}^{\alpha}\|_{3,2}\|b^{12}\|_{1,2}.$

Collecting the information above yields the estimate

(4.13) $\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|b^{12}\|_{1,2}^{2}$
$\displaystyle\lesssim\big{(}\|\mathbf{u}^{12}\|_{1,2}\|\overline{b}^{\alpha}\|_{3,2},+\|\mathbf{u}\|_{3,2}\|b^{12}\|_{1,2}\big{)}\|b^{12}\|_{1,2}.$

For the equation for $\omega^{12}$, we aim to derive a square-integrable
estimate in space. For this, we first note that

(4.14)
$\displaystyle\Big{|}\big{\langle}\mathbf{u}\cdot\nabla\omega^{12}\,,\,\omega^{12}\big{\rangle}\Big{|}=0.$

Next, we have that

(4.15)
$\displaystyle\Big{|}\big{\langle}(\mathbf{u}-\mathbf{u}_{h})\cdot\nabla
b^{12}\,,\,\omega^{12}\big{\rangle}\Big{|}\lesssim\|\mathbf{u}-\mathbf{u}_{h}\|_{\infty}\|b^{12}\|_{1,2}\|\omega^{12}\|_{2}.$

Finally, we also have that

(4.16)
$\displaystyle\Big{|}\big{\langle}\mathbf{u}^{12}\cdot\nabla(\overline{\omega}^{\alpha}-\overline{b}^{\alpha})\,,\,\omega^{12}\big{\rangle}\Big{|}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{4}\|\nabla(\overline{\omega}^{\alpha}-\overline{b}^{\alpha})\|_{4}\|\omega^{12}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{1,2}\|\overline{\omega}^{\alpha}-\overline{b}^{\alpha}\|_{2,2}\|\omega^{12}\|_{2}.$

Therefore,

(4.17)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\omega^{12}\|_{2}^{2}$
$\displaystyle\lesssim\big{(}\|\mathbf{u}-\mathbf{u}_{h}\|_{\infty}\|b^{12}\|_{1,2}+\|\mathbf{u}^{12}\|_{1,2}\|\overline{\omega}^{\alpha}-\overline{b}^{\alpha}\|_{2,2}\big{)}\|\omega^{12}\|_{2}.$

Summing up with the estimate for $b^{12}$ then yield

(4.18)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{(}\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}\big{)}$
$\displaystyle\lesssim\big{(}\|\mathbf{u}^{12}\|_{1,2}\|\overline{b}^{\alpha}\|_{3,2},+\|\mathbf{u}\|_{3,2}\|b^{12}\|_{1,2}\big{)}\|b^{12}\|_{1,2}$
$\displaystyle+\big{(}\|\mathbf{u}-\mathbf{u}_{h}\|_{\infty}\|b^{12}\|_{1,2}+\|\mathbf{u}^{12}\|_{1,2}\|\overline{\omega}^{\alpha}-\overline{b}^{\alpha}\|_{2,2}\big{)}\|\omega^{12}\|_{2}.$

If we use (2.7) and (4.7), we obtain

(4.19)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{(}\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}\big{)}$
$\displaystyle\lesssim\big{(}\|\mathbf{u}^{12}\|_{1,2}+\|b^{12}\|_{1,2}\big{)}\|b^{12}\|_{1,2}$
$\displaystyle+\big{(}\|b^{12}\|_{1,2}+\|\mathbf{u}^{12}\|_{1,2}\big{)}\|\omega^{12}\|_{2}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{1,2}^{2}+\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}.$

Now recall that the velocity fields for the $\alpha$-TQG, the
$\bf{i}\alpha$-TQG, and the TQG are given by

(4.20)
$\displaystyle{\mathbf{u}}^{\alpha}=\nabla^{\perp}(1-\alpha\Delta)^{-1}(\Delta-1)^{-1}(\omega^{\alpha}-f),$
(4.21)
$\displaystyle\overline{\mathbf{u}}^{\alpha}=\nabla^{\perp}(1-\alpha\Delta)^{-1}(\Delta-1)^{-1}(\omega-f),$
(4.22) $\displaystyle\mathbf{u}=\nabla^{\perp}(\Delta-1)^{-1}(\omega-f)$

respectively. As a result, in particular, the difference

$\displaystyle\overline{\mathbf{u}}^{\alpha}-\mathbf{u}=\alpha\Delta\nabla^{\perp}(1-\alpha\Delta)^{-1}(\Delta-1)^{-1}(\omega-f)$

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)
$\displaystyle\|\overline{\mathbf{u}}^{\alpha}-\mathbf{u}\|_{1,2}^{2}\lesssim\alpha^{2}\|\mathbf{u}\|_{3,2}^{2}\lesssim\alpha^{2}\|\omega-f\|_{2,2}^{2}\lesssim\alpha^{2}\big{(}\|\omega\|_{2,2}^{2}+\|f\|_{2,2}^{2}\big{)}\lesssim\alpha^{2}.$

This means that for the comparison of the TQG and $\bf{i}\alpha$-TQG, we can
conclude from (4.19) that

(4.24)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{(}\|\overline{b}^{\alpha}-b\|_{1,2}^{2}$
$\displaystyle+\|\overline{\omega}^{\alpha}-\omega\|_{2}^{2}\big{)}\lesssim\alpha^{2}+\|\overline{b}^{\alpha}-b\|_{1,2}^{2}+\|\overline{\omega}^{\alpha}-\omega\|_{2}^{2}$

and by Grönwall’s lemma, (note that $b^{12}_{0}=0$ and $\omega^{12}_{0}=0$)

(4.25) $\displaystyle\|(\overline{b}^{\alpha}-b)(t)\|_{1,2}^{2}$
$\displaystyle+\|(\overline{\omega}^{\alpha}-\omega)(t)\|_{2}^{2}\lesssim\alpha^{2}T\exp(cT),\qquad
t\in[0,T].$

The above gives the decay rate of the difference of the solution to the
$\mathbf{i}\alpha$-TQG and TQG equations. Our next goal is to obtain a decay
rate for the difference between the $\mathbf{i}\alpha$-TQG and the
$\alpha$-TQG equations. Recall that they share the same initial data. For
this, we use the estimate

$\displaystyle\|\overline{\mathbf{u}}^{\alpha}-{\mathbf{u}}^{\alpha}\|_{1,2}^{2}\lesssim\|\omega-\omega^{\alpha}\|_{2}^{2}\lesssim\|\overline{\omega}^{\alpha}-\omega\|_{2}^{2}+\|\overline{\omega}^{\alpha}-\omega^{\alpha}\|_{2}^{2}$

so that by setting $b^{12}:=\overline{b}^{\alpha}-b^{\alpha}$,
$\omega^{12}:=\overline{\omega}^{\alpha}-\omega^{\alpha}$ and
$\mathbf{u}^{12}:=\overline{\mathbf{u}}^{\alpha}-{\mathbf{u}}^{\alpha}$, we
obtain from (4.19),

(4.26)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{(}\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}\big{)}$
$\displaystyle\lesssim\|\mathbf{u}^{12}\|_{1,2}^{2}+\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}$
$\displaystyle\lesssim\|\overline{\omega}^{\alpha}-\omega\|_{2}^{2}+\|b^{12}\|_{1,2}^{2}+\|\omega^{12}\|_{2}^{2}.$

Since both systems share the same initial conditions, by Grönwall’s lemma and
(4.25), it follows that

(4.27)
$\displaystyle\|(\overline{b}^{\alpha}-b^{\alpha})(t)\|_{1,2}^{2}+\|(\overline{\omega}^{\alpha}-\omega^{\alpha})(t)\|_{2}^{2}$
$\displaystyle\leq\exp(cT)\int_{0}^{T}c\,\|(\overline{\omega}^{\alpha}-\omega)(t)\|_{2}^{2}\,\mathrm{d}t$
$\displaystyle\lesssim\alpha^{2}T\big{[}\exp(cT)\big{]}^{2}.$

By the triangle inequality, it follows from (4.25) and (4.27) that

(4.28) $\displaystyle\|({b}^{\alpha}-b)(t)\|_{1,2}^{2}$
$\displaystyle+\|({\omega}^{\alpha}-\omega)(t)\|_{2}^{2}\lesssim\alpha^{2}T\exp(cT)\big{[}1+\exp(cT)\big{]}$

so that

$\displaystyle\sup_{t\in[0,T]}\|({b}^{\alpha}-b)(t)\|_{1,2}\rightarrow
0,\quad\qquad\sup_{t\in[0,T]}\|({\omega}^{\alpha}-\omega)(t)\|_{2}\rightarrow
0$

as $\alpha\rightarrow 0$. This ends the proof.

######  Remark 4.3.

We conjecture that Proposition 4.1 may be extended to the stronger space
$W^{2,2}(\mathbb{T}^{2})\times W^{1,2}(\mathbb{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.

##  5\. Conditions for blowup of the $\alpha$-TQG model

In the following, we give Beale–Kato–Majda [2] type conditions under which we
expect the strong solution of the $\alpha$-TQG to blowup.

######  Theorem 5.1.

Fix $\alpha>0$. Under Assumption 2.4, let
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ be a maximal solution of
(2.12)–(2.14). If $T_{\max}^{\alpha}<\infty$, then

$\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla
b^{\alpha}\|_{\infty}\,\mathrm{d}t=\infty.$

###### Proof.

First of all, fix $\alpha>0$ and let $T_{\max}^{\alpha}>0$ be the maximal time
so that

(5.1) $\displaystyle\limsup_{T_{n}\rightarrow
T_{\max}^{\alpha}}\|(b^{\alpha},\omega^{\alpha})(T_{n})\|_{\mathcal{M}}^{2}=\infty.$

We now suppose that

(5.2) $\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla
b^{\alpha}\|_{\infty}\,\mathrm{d}t=K<\infty$

and show that

(5.3)
$\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}\lesssim
1,\qquad t<T_{\max}^{\alpha}$

holds yielding a contradiction to (5.1).
Before we show the estimate (5.3), we first require preliminary estimates for
$\|\nabla\mathbf{u}^{\alpha}\|_{\infty}$ and $\|\omega^{\alpha}\|_{2}$. To
obtain these estimates, we note that the Fourier multiplier $m^{\alpha}(\xi)$
defined below satisfies the bound

(5.4) $\displaystyle
m^{\alpha}(\xi):=\frac{(1+|\xi|^{2})|\xi|^{2}}{(1+\alpha|\xi|^{2})^{2}}\lesssim_{\alpha}1$

for all fixed $\alpha>0$ and in particular, the bound may be taken uniformly
of all $\alpha\geq 1$. Due to (5.4) and the continuous embedding
$W^{3,2}(\mathbb{T}^{2})\hookrightarrow W^{1,\infty}(\mathbb{T}^{2})$, we can
conclude that

(5.5)
$\displaystyle\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\lesssim\|\mathbf{u}^{\alpha}\|_{3,2}\lesssim_{\alpha}\|\omega^{\alpha}-f\|_{2}\lesssim\|\omega^{\alpha}-f\|_{2,2}.$

On the other hand, if we test (2.13) with $\omega^{\alpha}$ and use (2.7) for
$w=(1-\alpha\Delta)^{-1}(\omega^{\alpha}-f)$ and $k=0$, we obtain

(5.6)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\omega^{\alpha}\|_{2}^{2}$
$\displaystyle\lesssim\Big{(}\|\mathbf{u}^{\alpha}\|_{2}\|\omega^{\alpha}\|_{2}+\|\mathbf{u}_{h}\|_{2}\|\omega^{\alpha}\|_{2}\Big{)}\|\nabla
b^{\alpha}\|_{\infty}\lesssim\Big{(}1+\|\omega^{\alpha}\|_{2}^{2}\Big{)}\|\nabla
b^{\alpha}\|_{\infty}$

for a constant depending only on $\|\mathbf{u}_{h}\|_{2}$ and $\|f\|_{2}$. It
therefore follow from (5.6) and (5.2) that

(5.7)
$\displaystyle\|\omega^{\alpha}(t)\|_{2}^{2}\lesssim\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK),\qquad
t<T_{\max}^{\alpha}$

which when combined with the first estimate in (5.5) yields

(5.8)
$\displaystyle\|\nabla\mathbf{u}^{\alpha}(t)\|_{\infty}\lesssim_{\alpha}1+\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK)\lesssim_{\alpha}\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK),\qquad
t<T_{\max}^{\alpha}.$

Recall that $f\in W^{2,2}(\mathbb{T}^{2})$ by assumption.
With these preliminary estimates in hand, we can proceed to derive estimates
for $(b^{\alpha},\omega^{\alpha})$ in the space $\mathcal{M}$ of existence of
a solution. Since the space of smooth functions is dense in $\mathcal{M}$, it
suffices to show our result for a smooth solution pair
$(b^{\alpha},\omega^{\alpha})$.
To achieve our goal, we apply $\partial^{\beta}$ to (2.12) for $|\beta|\leq 3$
to obtain

(5.9) $\displaystyle\frac{\partial}{\partial
t}\partial^{\beta}b^{\alpha}+\mathbf{u}^{\alpha}\cdot\nabla\partial^{\beta}b^{\alpha}=R_{1}$

where

$\displaystyle R_{1}:=\mathbf{u}^{\alpha}\cdot\partial^{\beta}\nabla
b^{\alpha}-\partial^{\beta}(\mathbf{u}^{\alpha}\cdot\nabla b^{\alpha}).$

Now since $\mathrm{div}\mathbf{u}^{\alpha}=0$, if we multiply (5.9) by
$\partial^{\beta}b^{\alpha}$ and sum over the multiindex $\beta$ so that
$|\beta|\leq 3$, we obtain

(5.10)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|b^{\alpha}\|_{3,2}^{2}$
$\displaystyle\lesssim\Big{(}\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\|b^{\alpha}\|_{3,2}+\|\nabla
b^{\alpha}\|_{\infty}\|\mathbf{u}^{\alpha}\|_{3,2}\Big{)}\|b^{\alpha}\|_{3,2}$
$\displaystyle\lesssim(\|\nabla\mathbf{u}^{\alpha}\|_{\infty}+\|\nabla
b^{\alpha}\|_{\infty})(1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2})$

where we have used (2.7) for $w=(1-\alpha\Delta)^{-1}(\omega^{\alpha}-f)$ and
$k=2$.
Next, we find a bound for $\|\omega^{\alpha}\|^{2}_{2,2}$. For this, we apply
$\partial^{\beta}$ to (2.13) for $|\beta|\leq 2$ and we obtain

(5.11) $\displaystyle\frac{\partial}{\partial
t}\partial^{\beta}\omega^{\alpha}+\mathbf{u}^{\alpha}\cdot\nabla\partial^{\beta}(\omega^{\alpha}-b^{\alpha})+\mathbf{u}_{h}\cdot\nabla\partial^{\beta}b^{\alpha}=R_{2}+R_{3}+R_{4}$

where

$\displaystyle R_{2}$
$\displaystyle:=\mathbf{u}^{\alpha}\cdot\partial^{\beta}\nabla\omega^{\alpha}-\partial^{\beta}(\mathbf{u}^{\alpha}\cdot\nabla\omega^{\alpha}),$
$\displaystyle R_{3}$
$\displaystyle:=-\mathbf{u}^{\alpha}\cdot\partial^{\beta}\nabla
b^{\alpha}+\partial^{\beta}(\mathbf{u}^{\alpha}\cdot\nabla b^{\alpha}),$
$\displaystyle R_{4}$
$\displaystyle:=\mathbf{u}_{h}\cdot\partial^{\beta}\nabla
b^{\alpha}-\partial^{\beta}(\mathbf{u}_{h}\cdot\nabla b^{\alpha})$

are such that

(5.12)
$\displaystyle\big{|}\big{\langle}R_{2}\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}\big{|}\lesssim\|R_{2}\|_{2}\|\omega^{\alpha}\|_{2,2}$
$\displaystyle\lesssim\Big{(}\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\|\omega^{\alpha}\|_{2,2}+\|\omega^{\alpha}\|_{2,2}(1+\|\omega^{\alpha}\|_{2})\Big{)}\|\omega^{\alpha}\|_{2,2},$
(5.13)
$\displaystyle\big{|}\big{\langle}R_{3}\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}\big{|}\lesssim\|R_{3}\|_{2}\|\omega^{\alpha}\|_{2,2}$
$\displaystyle\lesssim\Big{(}\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\|b^{\alpha}\|_{3,2}+\|\nabla
b^{\alpha}\|_{\infty}(1+\|\omega^{\alpha}\|_{2,2})\Big{)}\|\omega^{\alpha}\|_{2,2},$
(5.14)
$\displaystyle\big{|}\big{\langle}R_{4}\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}\big{|}\lesssim\|R_{4}\|_{2}\|\omega^{\alpha}\|_{2,2}$
$\displaystyle\lesssim\Big{(}\|b^{\alpha}\|_{3,2}+\|\nabla
b^{\alpha}\|_{\infty}\Big{)}\|\omega^{\alpha}\|_{2,2}$

holds true. Recall that $\mathbf{u}_{h}\in W^{3,2}(\mathbb{T}^{2})$ by
assumption and note that we have used the continuous embedding
$W^{3,2}(\mathbb{T}^{2})\hookrightarrow W^{2,4}(\mathbb{T}^{2})$ and the
second inequality in (5.5) in order to obtain the term
$1+\|\omega^{\alpha}\|_{2}$ for the estimates (5.12). Next, by using
$\mathrm{div}\mathbf{u}^{\alpha}=0$, we have the following identity

(5.15)
$\displaystyle\big{\langle}(\mathbf{u}^{\alpha}\cdot\nabla\partial^{\beta}\omega^{\alpha})\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}=\frac{1}{2}\int_{\mathbb{T}^{2}}\mathrm{div}(\mathbf{u}^{\alpha}|\partial^{\beta}\omega^{\alpha}|^{2})\,\mathrm{d}x=0.$

Additionally, the following estimates for $L^{2}$ inner products holds true

(5.16)
$\displaystyle\Big{|}\big{\langle}(\mathbf{u}^{\alpha}\cdot\nabla\partial^{\beta}b^{\alpha})\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}\Big{|}$
$\displaystyle\lesssim\|\mathbf{u}^{\alpha}\|_{\infty}\|b^{\alpha}\|_{3,2}^{2}+\|\mathbf{u}^{\alpha}\|_{\infty}\|\omega^{\alpha}\|_{2,2}^{2},$
(5.17)
$\displaystyle\Big{|}\big{\langle}(\mathbf{u}_{h}\cdot\nabla\partial^{\beta}b^{\alpha})\,,\,\partial^{\beta}\omega^{\alpha}\big{\rangle}\Big{|}$
$\displaystyle\lesssim\|b^{\alpha}\|_{3,2}^{2}+\|\omega^{\alpha}\|_{2,2}^{2}.$

since $\mathbf{u}_{h}\in W^{3,2}(\mathbb{T}^{2})$. If we now collect the
estimates above (keeping in mind that $f\in W^{2,2}(\mathbb{T}^{2})$ and
$\mathbf{u}_{h}\in W^{3,2}(\mathbb{T}^{2})$), we obtain by multiplying (5.11)
by $\partial^{\beta}\omega^{\alpha}$ with $|\beta|\leq 2$ and then summing
over $|\beta|\leq 2$, the following

(5.18)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\omega^{\alpha}\|_{2,2}^{2}$
$\displaystyle\lesssim\big{(}1+\|\mathbf{u}^{\alpha}\|_{1,\infty}+\|\nabla
b^{\alpha}\|_{\infty}+\|\omega^{\alpha}\|_{2}\big{)}\big{(}1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2}\big{)}.$

Summing up (5.10) and (5.18) and using (5.7)–(5.5) yields

(5.19)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}(1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2})$
$\displaystyle\lesssim_{\alpha}\Big{(}\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK)+\|\nabla
b^{\alpha}\|_{\infty}\Big{)}(1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2})$

so that by Grönwall’s lemma and (5.2), we obtain

(5.20) $\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}$
$\displaystyle\leq\big{[}1+\|(b_{0},\omega_{0})\|_{\mathcal{M}}^{2}\big{]}\exp\Big{(}c(\alpha)T_{\max}\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK)+c(\alpha)K\Big{)}$
$\displaystyle\lesssim_{T_{\max}^{\alpha},K,f,\mathbf{u}_{h},\omega_{0},b_{0},\alpha}1$

for all $t<{T_{\max}^{\alpha}}$ contradicting (5.1). ∎

######  Remark 5.2 (Global existence for constant-in-space buoyancy).

We note that for the $\alpha$-TQG model, when the buoyancy is constant in
space so that $\nabla b^{\alpha}=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
$\nabla b^{\alpha}=0$.

######  Remark 5.3 (Global existence for super diffusive buoyancy).

Fix $\alpha>0$. Note that by adding any super diffusive term
$\Lambda^{\beta}b^{\alpha}$, $\beta\geq 1$ to the right-hand side of (2.12) so
that

$\displaystyle\int_{0}^{t}\|\nabla
b^{\alpha}(s)\|_{\infty}\,\mathrm{d}s\lesssim\int_{0}^{t}\|\Lambda^{\beta/2}b^{\alpha}(s)\|_{2}^{2}\,\mathrm{d}s$

holds for all $t<T_{\max}^{\alpha}$, then we obtain a global solution since in
this case, we obtain the energy estimate

$\displaystyle\|b^{\alpha}(t)\|_{2}^{2}+\int_{0}^{t}\|\Lambda^{\beta/2}b^{\alpha}(s)\|_{2}^{2}\,\mathrm{d}s\leq\|b_{0}\|_{2}\leq\|b_{0}\|_{3,2}$

for all $t<T_{\max}^{\alpha}$.

Next, we also give an alternative to the blowup condition in Theorem 5.1
above.

######  Theorem 5.4.

Fix $\alpha>0$. Under Assumption 2.4, let
$(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ be a maximal solution of
(2.12)–(2.14). If $T_{\max}^{\alpha}<\infty$, then

$\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\,\mathrm{d}t=\infty.$

###### Proof.

In analogy with the proof of Theorem 5.1, we fix $\alpha>0$, let
$T_{\max}^{\alpha}>0$ be the maximal time so that

(5.21) $\displaystyle\limsup_{T_{n}\rightarrow
T_{\max}^{\alpha}}\|(b^{\alpha},\omega^{\alpha})(T_{n})\|_{\mathcal{M}}^{2}=\infty.$

We now suppose that

(5.22)
$\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\,\mathrm{d}t=K<\infty$

and show that

(5.23)
$\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}\lesssim
1,\qquad t<T_{\max}^{\alpha}$

holds yielding a contradiction to (5.21).
To show this, we first fix $\alpha>0$. If we differentiate (2.12) in space,
multiply the resulting equation by $p|\nabla b^{\alpha}|^{p-2}\nabla
b^{\alpha}$ where $p>1$ is fixed and finite in this instant, and then
integrate over $\mathbb{T}^{2}$, we obtain

(5.24) $\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\nabla
b^{\alpha}\|_{p}^{p}\leq p\|\nabla\mathbf{u}^{\alpha}\|_{\infty}\|\nabla
b^{\alpha}\|_{p}^{p}.$

Rather than use Grönwall’s lemma at this point, we use the analogous
separation of variable technique to obtain

$\displaystyle\frac{1}{p}\ln\|\nabla
b^{\alpha}(t)\|_{p}^{p}\leq\frac{1}{p}\ln\|\nabla
b_{0}\|_{p}^{p}+\int_{0}^{t}\|\nabla\mathbf{u}^{\alpha}(s)\|_{\infty}\,\mathrm{d}s$

so that

$\displaystyle\|\nabla b^{\alpha}(t)\|_{p}\leq\|\nabla
b_{0}\|_{p}\exp\bigg{(}\int_{0}^{t}\|\nabla\mathbf{u}^{\alpha}(s)\|_{\infty}\,\mathrm{d}s\bigg{)}$

holds for all $t<T_{\max}^{\alpha}$ uniformly in $p>1$. Since $b_{0}\in
W^{3,2}(\mathbb{T}^{2})$, by Sobolev embedding and (5.22), we obtain

(5.25) $\displaystyle\|\nabla
b^{\alpha}(t)\|_{\infty}\lesssim\|b_{0}\|_{3,2}\exp(K)$

for all $t<T_{\max}^{\alpha}$ so that

$\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla
b^{\alpha}(t)\|_{\infty}\,\mathrm{d}t\leq
c\,T_{\max}^{\alpha}\|b_{0}\|_{3,2}\exp(K).$

We can therefore deduce from the estimate (5.6) that

(5.26)
$\displaystyle\|\omega^{\alpha}(t)\|_{2}^{2}\lesssim\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp\big{(}cT_{\max}^{\alpha}\|b_{0}\|_{3,2}\exp(K)\big{)},\qquad
t<T_{\max}^{\alpha}.$

Summing up (5.10) and (5.18) and using (5.25)–(5.36) yields

(5.27)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}(1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2})$
$\displaystyle\lesssim\big{(}1+\|\mathbf{u}^{\alpha}\|_{1,\infty}+\mathfrak{A}\big{)}(1+\|(b^{\alpha},\omega^{\alpha})\|_{\mathcal{M}}^{2})$

where

$\displaystyle\mathfrak{A}:=\|b_{0}\|_{3,2}\exp(K)+\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp\big{(}cT_{\max}^{\alpha}\|b_{0}\|_{3,2}\exp(K)\big{)}.$

We can therefore conclude from Grönwall’s lemma and (5.22) that

(5.28) $\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}$
$\displaystyle\leq\big{[}1+\|(b_{0},\omega_{0})\|_{\mathcal{M}}^{2}\big{]}\exp\Big{(}cT_{\max}^{\alpha}+cK+cT_{\max}^{\alpha}\mathfrak{A}\Big{)}$
$\displaystyle\lesssim_{T_{\max}^{\alpha},K,f,\mathbf{u}_{h},\omega_{0},b_{0}}1$

for all $t<{T_{\max}^{\alpha}}$ contradicting (5.21). ∎

In order the give the final blowup criterion, let us now recall the following
endpoint Sobolev inequality by Brézis and Gallouet [8].

######  Lemma 5.5.

If $f\in W^{2,2}(\mathbb{T}^{2})$ then

$\displaystyle\|f\|_{\infty}\lesssim(1+\|f\|_{1,2})\sqrt{\ln(\mathrm{e}+\|f\|_{2,2})}.$

######  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.

######  Theorem 5.7.

Fix $\alpha>0$. Let $(b^{\alpha},\omega^{\alpha},T_{\max}^{\alpha})$ be a
maximal solution of (2.12)–(2.14). If $T_{\max}^{\alpha}<\infty$, then

$\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla
b^{\alpha}\|_{1,2}\,\mathrm{d}t=\infty.$

###### Proof.

Fix $\alpha>0$ and let $T_{\max}^{\alpha}>0$ be the maximal time so that

(5.29) $\displaystyle\limsup_{T_{n}\rightarrow
T_{\max}^{\alpha}}\|(b^{\alpha},\omega^{\alpha})(T_{n})\|_{\mathcal{M}}^{2}=\infty.$

We now suppose that

(5.30) $\displaystyle\int_{0}^{T_{\max}^{\alpha}}\|\nabla
b^{\alpha}\|_{1,2}\,\mathrm{d}t=K<\infty$

and show that

(5.31)
$\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}\lesssim
1,\qquad t<T_{\max}^{\alpha}$

holds and thereby yields a contradiction to (5.29).
To show this, we first fix $\alpha>0$ and define

$\displaystyle
g(t):=\mathrm{e}+\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2},\qquad
t>0.$

With this definition in hand, it follows from estimates (5.10) and (5.18) that

(5.32)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}g\lesssim\big{(}1+\|\mathbf{u}^{\alpha}\|_{1,\infty}+\|\omega^{\alpha}\|_{2}+\|\nabla
b^{\alpha}\|_{\infty}\big{)}g.$

By using Lemma 5.5 and the monotonic property of logarithms however, we can
deduce that

(5.33) $\displaystyle\|\nabla b^{\alpha}\|_{\infty}\lesssim\big{(}1+\|\nabla
b^{\alpha}\|_{1,2}\big{)}\ln\,g.$

Next, we observe that the second estimate in (5.5) yields

(5.34) $\displaystyle
1+\|\mathbf{u}^{\alpha}\|_{1,\infty}+\|\omega^{\alpha}\|_{2}\lesssim_{\alpha}1+\|\omega^{\alpha}\|_{2}$

for $f\in L^{2}(\mathbb{T}^{2})$. However, if we test (2.13) with
$\omega^{\alpha}$, we also obtain

(5.35)
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\omega^{\alpha}\|_{2}^{2}$
$\displaystyle\lesssim\Big{(}\|\mathbf{u}^{\alpha}\|_{\infty}\|\omega^{\alpha}\|_{2}+\|\mathbf{u}_{h}\|_{\infty}\|\omega^{\alpha}\|_{2}\Big{)}\|\nabla
b^{\alpha}\|_{2}\lesssim\Big{(}1+\|\omega^{\alpha}\|_{2}^{2}\Big{)}\|\nabla
b^{\alpha}\|_{1,2}$

for a constant depending only on $\|\mathbf{u}_{h}\|_{\infty}$ and
$\|f\|_{2}$. It therefore follow from (5.30) that

(5.36)
$\displaystyle\|\omega^{\alpha}(t)\|_{2}^{2}\lesssim\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK),\qquad
t<T_{\max}^{\alpha}.$

Using the fact that $\ln\,g\geq 1$, we can conclude from (5.34) and (5.36)
that

(5.37) $\displaystyle
1+\|\mathbf{u}^{\alpha}\|_{1,\infty}+\|\omega^{\alpha}\|_{2}\lesssim_{\alpha}\big{(}1+\|\omega_{0}\|_{2}^{2}\big{)}\exp(cK)\ln\,g.$

If we now combine this estimate with (5.33), we can conclude from (5.32) that

$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}g\lesssim\big{(}1+\|\omega_{0}\|_{2}+\|\nabla
b^{\alpha}\|_{1,2}\big{)}\exp(cK)\,g\ln\,g$

so that

$\displaystyle g(t)\leq g(0)^{\exp\int_{0}^{t}c(1+\|\omega_{0}\|_{2}+\|\nabla
b^{\alpha}(s)\|_{1,2})\exp(cK)\,\mathrm{d}s}$

holds for $t<T_{\max}^{\alpha}$. In particular, given (5.30), we have shown
that

(5.38) $\displaystyle\|(b^{\alpha},\omega^{\alpha})(t)\|_{\mathcal{M}}^{2}$
$\displaystyle\leq\big{[}\mathrm{e}+\|(b^{\alpha}_{0},\omega^{\alpha}_{0})\|_{\mathcal{M}}^{2}\big{]}^{\exp(c(T_{\max}^{\alpha}+K)\exp(cK))}$
$\displaystyle\lesssim_{T_{\max}^{\alpha},K,f,\mathbf{u}_{h},\omega_{0},b_{0}}1$

for all $t<{T_{\max}^{\alpha}}$ contradicting (5.29). ∎

##  6\. 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 $\alpha$–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 $\alpha$-TQG system. By setting $\alpha$ to
zero we obtain our numerical setup for the TQG system.

Consider a bounded domain ${\mathcal{D}}$. Let $\partial{\mathcal{D}}$ denote
the boundary. We impose Dirichlet boundary conditions

(6.1) $\psi^{\alpha}=0,\ \Delta\psi^{\alpha}=0,\qquad\text{on
}\partial{\mathcal{D}}.$

For ${\mathcal{D}}=\mathbb{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}\left({\mathcal{D}}\right)$ denote the Sobolev
$W^{1,2}\left({\mathcal{D}}\right)$ space and let
$\left\|.\right\|_{\partial{\mathcal{D}}}$ denote the
$L^{2}\left(\partial{\mathcal{D}}\right)$ norm. Define the space

(6.2) $W^{1}\left({\mathcal{D}}\right):=\left\\{\nu\in
H^{1}\left({\mathcal{D}}\right)\left|\left\|\nu\right\|_{\partial{\mathcal{D}}}=0\right.\right\\}.$

We express (2.14) as two inhomogeneous Helmholtz equations

(6.3) $\displaystyle\omega^{\alpha}-f$ $\displaystyle=(\Delta-1)\tilde{\psi}$
(6.4) $\displaystyle-\tilde{\psi}$
$\displaystyle=(\alpha\Delta-1)\psi^{\alpha}.$

We take $\tilde{\psi},\psi^{\alpha}\in W^{1}({\mathcal{D}})$. Since the two
equations are of the same form, let us first consider (6.3). Using an
arbitrary test function $\phi\in W^{1}({\mathcal{D}})$, we obtain the
following weak form of (6.3),

(6.5)
$\displaystyle\langle\nabla\tilde{\psi},\nabla\phi\rangle_{\mathcal{D}}+\langle\tilde{\psi},\phi\rangle_{\mathcal{D}}$
$\displaystyle=-\langle\omega^{\alpha}-f,\phi\rangle_{\mathcal{D}}.$

Define the functionals

(6.6) $\displaystyle L_{\alpha}(v,\phi)$ $\displaystyle:=\alpha\langle\nabla
v,\nabla\phi\rangle_{\mathcal{D}}+\langle v,\phi\rangle_{\mathcal{D}}$ (6.7)
$\displaystyle F_{\cdot}(\phi)$
$\displaystyle:=-\langle\cdot,\phi\rangle_{\mathcal{D}}$

for $v,\phi\in W^{1}({\mathcal{D}})$, then (6.5) can be written as

(6.8) $L_{1}(\tilde{\psi},\phi)=F_{\omega^{\alpha}-f^{\alpha}}(\phi).$

And similarly for (6.4), we have

(6.9) $L_{\alpha}(\psi^{\alpha},\phi)=F_{-\tilde{\psi}}(\phi).$

We discretise equations (6.8) and (6.9) using a continuous Galerkin (CG)
discretisation scheme.

Let $\delta$ be the discretisation parameter, and let ${\mathcal{D}}_{\delta}$
denote a space filling triangulation of the domain, that consists of geometry-
conforming nonoverlapping elements. Define the approximation space

(6.10) $W_{\delta}^{k}({\mathcal{D}}):=\left\\{\phi_{\delta}\in
W^{1}\left({\mathcal{D}}\right)\ :\ \phi_{\delta}\in
C\left({\mathcal{D}}\right),\left.\phi_{\delta}\right|_{K}\in\Pi^{k}\left(K\right)\text{
each }K\in{\mathcal{D}}_{\delta}\right\\}.$

in which $C({\mathcal{D}})$ is the space of continuous functions on
${\mathcal{D}}$, and $\Pi^{k}\left(K\right)$ denotes the space of polynomials
of degree at most $k$ on element $K\in{\mathcal{D}}_{\delta}$.

For (6.8), given $f_{\delta}\in W_{\delta}^{k}({\mathcal{D}})$ and
$\omega_{\delta}\in V_{\delta}^{k}({\mathcal{D}})$ (see (6.13) for the
definition of $V_{\delta}^{k}({\mathcal{D}})$), our numerical approximation is
the solution $\tilde{\psi}_{\delta}\in W_{\delta}^{k}({\mathcal{D}})$ that
satisfies

(6.11)
$L_{1}(\tilde{\psi}_{\delta},\phi_{\delta})=F_{\omega_{\delta}-f_{\delta}}(\phi_{\delta})$

for all test functions $\phi_{\delta}\in W_{\delta}^{k}({\mathcal{D}})$. Then,
using $\tilde{\psi}_{\delta}$, our numerical approximation of $\psi^{\alpha}$
is the solution $\psi_{\delta}^{\alpha}\in W_{\delta}^{k}({\mathcal{D}})$ that
satisfies

(6.12)
$L_{\alpha}(\psi_{\delta}^{\alpha},\phi_{\delta})=F_{-\tilde{\psi}_{\delta}}(\phi_{\delta})$

for all test functions $\phi_{\delta}\in W_{\delta}^{k}({\mathcal{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_{\delta}^{k}({\mathcal{D}})$,
to be the element-wise polynomial space,

(6.13) $V_{\delta}^{k}({\mathcal{D}})=\left\\{\left.v_{\delta}\in
L^{2}({\mathcal{D}})\right|\forall K\in{\mathcal{D}}_{\delta},\
\exists\phi_{\delta}\in\Pi^{k}(K):\
\left.v_{\delta}\right|_{K}=\left.\phi_{\delta}\right|_{K}\right\\}.$

We look to approximate $b^{\alpha}$ and $\omega^{\alpha}$ in the space
$V_{\delta}^{k}({\mathcal{D}})$. Essentially, this means our approximations of
$b^{\alpha}$ and $\omega^{\alpha}$ are each an direct sum over the elements in
${\mathcal{D}}_{\delta}$. Additional constraints on the numerical fluxes
across shared element boundaries are needed to ensure conservation properties
and stability. Further, note that $W^{k}_{\delta}({\mathcal{D}})\subset
V^{k}_{\delta}({\mathcal{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

(6.14) $\displaystyle\langle\partial_{t}b^{\alpha},\nu_{\delta}\rangle_{K}$
$\displaystyle=\langle b^{\alpha}{\bf
u}^{\alpha},\nabla\nu_{\delta}\rangle_{K}-\langle b^{\alpha}{\bf
u}^{\alpha}\cdot{\hat{\mathbf{n}}},\nu_{\delta}\rangle_{\partial K},\quad
K\in{\mathcal{D}}_{\delta}$

where $\nu_{\delta}\in V^{k}_{\delta}({\mathcal{D}})$ is any test function,
$\partial K$ denotes the boundary of $K$, and ${\hat{\mathbf{n}}}$ denotes the
unit normal vector to $\partial K$. Let $b^{\alpha}_{\delta}$ be the
approximation of $b^{\alpha}$ in $V_{\delta}^{k}$, and let
$\mathbf{u}^{\alpha}_{\delta}=\nabla^{\perp}\psi^{\alpha}_{\delta}$ for
$\psi^{\alpha}_{\delta}\in W^{k}_{\delta}$. Our discretised buoyancy equation
over each element is given by

(6.15) $\langle\partial_{t}b^{\alpha}_{\delta},\nu_{\delta}\rangle_{K}=\langle
b^{\alpha}_{\delta}{\bf
u}^{\alpha}_{\delta},\nabla\nu_{\delta}\rangle_{K}-\langle
b^{\alpha}_{\delta}{\bf
u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}},\nu_{\delta}\rangle_{\partial
K},\quad K\in{\mathcal{D}}_{\delta}.$

Similarly, let $\omega^{\alpha}_{\delta}\in V_{\delta}^{k}({\mathcal{D}})$ be
the approximation of $\omega^{\alpha}$, and let $h_{\delta}\in
W_{\delta}^{k}$. We obtain the following discretised variational formulation
that corresponds to (2.13),

$\displaystyle\langle\partial_{t}\omega^{\alpha}_{\delta},\nu_{\delta}\rangle_{K}$
$\displaystyle=\langle(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta}){\bf
u}^{\alpha}_{\delta},\nabla\nu_{\delta}\rangle_{K}-\langle(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta}){\bf
u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}},\nu_{\delta}\rangle_{\partial K}$
(6.16)
$\displaystyle\qquad-\frac{1}{2}\langle\nabla\cdot(b^{\alpha}_{\delta}\nabla^{\perp}h_{\delta}),\nu_{\delta}\rangle_{K},\qquad
K\in{\mathcal{D}}_{\delta},$

for test function $\nu_{\delta}\in V_{\delta}^{k}({\mathcal{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
${\mathcal{D}}_{\delta}$. In doing so, the $\partial K$ terms in (6.15) and
(6.16) must be treated carefully. Let $\partial K_{\text{ext}}$ denote the
part of cell boundary that is contained in $\partial{\mathcal{D}}$. Let
$\partial K_{\text{int}}$ denote the part of the cell boundary that is
contained in the interior of the domain
${\mathcal{D}}\backslash\partial{\mathcal{D}}$. On $\partial K_{\text{ext}}$
we simply impose the PDE boundary conditions. However, on $\partial
K_{\text{int}}$ we need to consider the contribution from each of the
neighbouring elements. By choice, the approximant $\psi^{\alpha}_{\delta}\in
W_{\delta}^{k}$ is continuous on $\partial K$. And since

(6.17)
$\mathbf{u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}=\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}=-\nabla\psi^{\alpha}_{\delta}\cdot\hat{\tau}=-\frac{\,\mathrm{d}\psi^{\alpha}_{\delta}}{\,\mathrm{d}\hat{\tau}},$

where $\hat{\tau}$ denotes the unit tangential vector to $\partial K$,
$\mathbf{u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}$ is also continuous. This
means $\mathbf{u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{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}^{\alpha}({\mathcal{D}})$,
$b^{\alpha}_{\delta}$ and $\omega^{\alpha}_{\delta}$ are multi-valued on
$\partial K_{\text{int}}$. Thus, as our approximation of $\omega^{\alpha}$ and
$b^{\alpha}$ over the whole domain is the sum over
$K\in{\mathcal{D}}_{\delta}$, we have to constrain the flux on the set
$\left(\bigcup_{K\in{\mathcal{D}}_{\delta}}\partial
K\right)\setminus\partial{\mathcal{D}}$. This is done using appropriately
chosen numerical flux fields in the boundary terms of (6.15) and (6.16).

Let $\nu^{-}:=\lim_{\epsilon\uparrow 0}\nu({\bf
x}+\epsilon{\hat{\mathbf{n}}})$ and $\nu^{+}:=\lim_{\epsilon\downarrow
0}\nu({\bf x}+\epsilon{\hat{\mathbf{n}}})$, for ${\bf x}\in\partial K$, be the
inside and outside (with respect to a fixed element $K$) values respectively,
of a function $\nu$ on the boundary. Let $\hat{f}$ be a _numerical flux_
function that satisfies the following properties:

  1. (i)

consistency

(6.18)
$\hat{f}(\nu,\nu,\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}})=\nu\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}}$

  2. (ii)

conservative

(6.19)
$\hat{f}(\nu^{+},\nu^{-},\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}})=-\hat{f}(\nu^{-},\nu^{+},-\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}})$

  3. (iii)

$L^{2}$ stable in the enstrophy norm with respect to the buoyancy equation,
see [4, Section 6].

With such an $\hat{f}$, we replace
$b^{\alpha}_{\delta}\mathbf{u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}$ by
the numerical flux $\hat{f}(b_{\delta}^{\alpha,+},b_{\delta}^{\alpha,-},{\bf
u}^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}})$ in (6.15). Similarly, in (6.16),
we replace
$(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})\mathbf{u}_{\delta}^{\alpha}\cdot{\hat{\mathbf{n}}}$
by
$\hat{f}((\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{+},(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{-},\mathbf{u}_{\delta}^{\alpha}\cdot{\hat{\mathbf{n}}})$.

######  Remark 6.1.

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,

(6.20)
$\hat{f}(\nu^{+},\nu^{-},\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}})=\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}}\\{\\{\nu\\}\\}-\frac{|\mathbf{u}_{\delta}\cdot{\hat{\mathbf{n}}}|}{2}\llbracket\nu\rrbracket$

where

(6.21)
$\\{\\{\nu\\}\\}:=\frac{1}{2}(\nu^{-}+\nu^{+}),\qquad\llbracket\nu\rrbracket:={\hat{\mathbf{n}}}^{-}\nu^{-}+{\hat{\mathbf{n}}}^{+}\nu^{+}.$

Finally, our goal is to find $b^{\alpha}_{\delta},\omega^{\alpha}_{\delta}\in
V_{\delta}^{k}({\mathcal{D}})$ such that for all $\nu_{\delta}\in
V_{\delta}^{k}({\mathcal{D}})$ we have

(6.22)
$\displaystyle\sum_{K\in{\mathcal{D}}_{\delta}}\langle\partial_{t}b^{\alpha}_{\delta},\nu_{\delta}\rangle_{K}$
$\displaystyle=\sum_{K\in{\mathcal{D}}_{\delta}}\big{\\{}\langle
b^{\alpha}_{\delta}\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla\nu_{\delta}\rangle_{K}-\langle\hat{f}_{b^{\alpha}}(b_{\delta}^{\alpha,+},b_{\delta}^{\alpha,-},\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}),\nu_{\delta}^{-}\rangle_{\partial
K}\big{\\}},$ (6.23)
$\displaystyle\sum_{K\in{\mathcal{D}}_{\delta}}\langle\partial_{t}\omega^{\alpha}_{\delta},\nu_{\delta}\rangle_{K}$
$\displaystyle=\sum_{K\in{\mathcal{D}}_{\delta}}\big{\\{}\langle(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla\nu_{\delta}\rangle_{K}-\frac{1}{2}\langle\nabla\cdot(b^{\alpha}_{\delta}\nabla^{\perp}h_{\delta}),\nu_{\delta}\rangle_{K}$
$\displaystyle\qquad-\langle\hat{f}_{\omega^{\alpha}}((\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{+},(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{-},\nabla^{\perp}\psi_{\delta}^{\alpha}\cdot{\hat{\mathbf{n}}}),\nu_{\delta}^{-}\rangle_{\partial
K}\big{\\}}$

with $\psi_{\delta}^{\alpha}\in W_{\delta}^{k}({\mathcal{D}})$ being the
numerical approximation to the stream function.

######  Remark 6.2.

In (6.22) and (6.23) we do not explicitly distinguish $\partial
K_{\text{ext}}$ and $\partial K_{\text{int}}$ because for the boundary
conditions (6.1), the $\partial K_{\text{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 $\alpha$–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 $\mathbb{T}^{2}$ .

######  Proposition 6.3 ($\alpha$–TQG conserved quantities).

On a bounded domain ${\mathcal{D}}$ with boundary $\partial{\mathcal{D}}$,
consider the $\alpha$-TQG system (2.12) – (2.14) with boundary conditions

(6.24) ${\hat{\mathbf{n}}}\cdot\mathbf{u}^{\alpha}=0,\
{\hat{\mathbf{n}}}\times\nabla b^{\alpha}=0,\
\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{\partial{\mathcal{D}}}\nabla\psi^{\alpha}\cdot{\hat{\mathbf{n}}}=0,\
\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{\partial{\mathcal{D}}}\nabla\Delta\psi^{\alpha}\cdot{\hat{\mathbf{n}}}=0,\
{\hat{\mathbf{n}}}\times\nabla\Delta\psi^{\alpha}=0$

We have

(6.25)
$E^{\alpha}(t):=-\frac{1}{2}\int_{\mathcal{D}}\\{(\omega^{\alpha}-f)\psi^{\alpha}+{h}b^{\alpha}\\}\,\mathrm{d}\mathbf{x},$

and

(6.26)
$C^{\alpha}_{\Psi,\Phi}(t):=\int_{\mathcal{D}}\\{\Psi(b^{\alpha})+\omega^{\alpha}\Phi(b^{\alpha})\\}\,\mathrm{d}\mathbf{x},\qquad\forall\Psi,\Phi\in
C^{\infty}$

are conserved, i.e. $\,\mathrm{d}_{t}E^{\alpha}(t)=0$ and
$\,\mathrm{d}_{t}C^{\alpha}_{\Psi,\Phi}(t)=0$.

###### Proof.

For the energy $E^{\alpha}(t)$, we obtain

$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}E^{\alpha}(t)$
$\displaystyle=-\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{\mathcal{D}}\\{\psi^{\alpha}(\Delta-1)(1-\alpha\Delta)\psi^{\alpha}+{h}b^{\alpha}\\}\,\mathrm{d}\mathbf{x}$
$\displaystyle=\int_{\mathcal{D}}\\{-\psi^{\alpha}\
\partial_{t}\omega^{\alpha}-\frac{1}{2}h\partial_{t}b^{\alpha}\\}\,\mathrm{d}\mathbf{x}.$

Substitute $\partial_{t}\omega^{\alpha}$ and $\partial_{t}b^{\alpha}$ using
(2.12) and (2.13), the result then follows from direct calculations.

Similarly for the casimirs, we obtain the result by directly evaluating
$\,\mathrm{d}_{t}C^{\alpha}_{\Psi,\Phi}(t)$. ∎

######  Remark 6.4.

In (6.24), except for the integral Neumann boundary conditions, we effectively
have Dirichlet boundary conditions for $\psi^{\alpha}$, $b^{\alpha}$ and
$\Delta\psi^{\alpha}$. The integral Neumann boundary conditions can be viewed
as imposing the Kelvin theorem on the boundary – consider

(6.27)
$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int_{\partial{\mathcal{D}}}\nabla(\psi^{\alpha}-\alpha\Delta\psi^{\alpha})\cdot{\hat{\mathbf{n}}}=0,$

and apply the divergence theorem.

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 $\psi^{\alpha}$ and $\Delta\psi^{\alpha}$
on $\partial{\mathcal{D}}$ to be _zero_.

We now analyse the conservation properties of our spatially discretised
$\alpha$-TQG system. Let $\langle\cdot,\cdot\rangle_{H^{1}({\mathcal{D}})}$ be
the $H^{1}({\mathcal{D}})$ inner product. Define the numerical total energy

(6.28)
$E^{\alpha}_{\delta}(t)=\frac{1}{2}\langle\psi^{\alpha}_{\delta},\tilde{\psi}_{\delta}\rangle_{H^{1}({\mathcal{D}})}-\frac{1}{2}\langle{h}_{\delta},b^{\alpha}_{\delta}\rangle_{\mathcal{D}},$

in which $\psi_{\delta}^{\alpha},\tilde{\psi}_{\delta}\in
W_{\delta}^{\alpha}({\mathcal{D}})$ and $b_{\delta}^{\alpha}\in
V_{\delta}^{k}({\mathcal{D}})$ are the numerical approximations of
$\psi^{\alpha}$, $\tilde{\psi}$ and $b^{\alpha}$ respectively. Define the
numerical casimir functional

(6.29)
$C_{\delta}^{\alpha}(t;\Psi,\Phi):=\langle\Psi(b_{\delta}^{\alpha}),1\rangle_{\mathcal{D}}+\langle\omega_{\delta}^{\alpha}\Phi(b_{\delta}^{\alpha})\rangle_{\mathcal{D}},\quad\Psi,\Phi\in
C^{\infty}.$

######  Lemma 6.5.

With boundary conditions (6.1)

(6.30) $\frac{\,\mathrm{d}}{\,\mathrm{d}t}E^{\alpha}_{\delta}(t)=0$

for $\alpha=0$.

In other words, the semi-discrete discretisation conserves energy in the TQG
case.

###### Proof.

We consider $\,\mathrm{d}/\,\mathrm{d}_{t}$ evaluated at an arbitrary fixed
value $t_{0}$. First note that when $\alpha=0$, we have
$\psi^{\alpha}_{\delta}=\tilde{\psi}_{\delta}$. Then, from (6.6) we have

(6.31) $\displaystyle 2E_{\delta}^{\alpha}(t)$
$\displaystyle=L_{1}(\psi_{\delta}^{\alpha},{\psi}^{\alpha}_{\delta})-\langle{h}_{\delta},b_{\delta}^{\alpha}\rangle_{\Omega}$

For the first term in (6.31), following (6.11) we have

(6.32)
$\displaystyle\left.L_{1}(\partial_{t}\psi_{\delta},\phi_{\delta})\right|_{t=t_{0}}=-\sum_{K\in\mathcal{{\mathcal{D}}}_{\delta}}\left.\langle\partial_{t}\omega^{\alpha}_{\delta},\phi_{\delta}\rangle_{K}\right|_{t=t_{0}},\quad\forall\phi_{\delta}\in
W_{\delta}^{k}({\mathcal{D}}).$

Thus, substituting in the discretised equation (6.23) for $\omega^{\alpha}$,
we obtain

$\displaystyle-
L_{1}(\left.\partial_{t}\psi_{\delta},\phi_{\delta})\right|_{t=t_{0}}$
$\displaystyle=\sum_{K\in{\mathcal{D}}_{\delta}}\left[\big{\\{}\langle(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla\phi_{\delta}\rangle_{K}-\frac{1}{2}\langle\nabla\cdot(b^{\alpha}_{\delta}\nabla^{\perp}h_{\delta}),\phi_{\delta}\rangle_{K}\right.$
(6.33)
$\displaystyle\qquad-\left.\langle\hat{f}_{\omega^{\alpha}}((\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{+},(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{-},\nabla^{\perp}\psi_{\delta}^{\alpha}\cdot{\hat{\mathbf{n}}}),\phi_{\delta}^{-}\rangle_{\partial
K}\big{\\}}\right]_{t=t_{0}},$

in which we can choose $\phi_{\delta}=[\psi^{\alpha}_{\delta}]_{t=t_{0}}$.
This choice is consistent with the discretisation space
$V_{\delta}^{k}({\mathcal{D}})$ for $\omega^{\alpha}$, since
$W_{\delta}^{k}({\mathcal{D}})\subset V_{\delta}^{k}({\mathcal{D}})$.

Similarly, for the second term in (6.31), using (6.22) we obtain

(6.34) $\left.\frac{\,\mathrm{d}}{\,\mathrm{d}t}\right|_{t=t_{0}}\langle
b^{\alpha}_{\delta},\nu_{\delta}\rangle_{\mathcal{D}}=\sum_{K\in{\mathcal{D}}_{\delta}}\left[\big{\\{}\langle
b^{\alpha}_{\delta}\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla\nu_{\delta}\rangle_{K}-\langle\hat{f}_{b^{\alpha}}(b_{\delta}^{\alpha,+},b_{\delta}^{\alpha,-},\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}),\nu_{\delta}^{-}\rangle_{\partial
K}\big{\\}}\right]_{t=t_{0}},$

where we can choose $\nu_{\delta}={h}_{\delta}$.

Putting together (6.31), (6.33) and (6.34), and substituting in ${h}_{\delta}$
and $[\psi^{\alpha}_{\delta}]_{t=t_{0}}$ for the arbitrary choices, we obtain

$\displaystyle
2\frac{\,\mathrm{d}}{\,\mathrm{d}t}\Big{|}_{t=t_{0}}E^{\alpha}_{\delta}(t)$
$\displaystyle=-\sum_{K\in{\mathcal{D}}_{\delta}}\left[2\Big{\\{}\langle(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})\underbrace{\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla\psi^{\alpha}_{\delta}}_{=0}\rangle_{K}-\frac{1}{2}\langle\nabla\cdot(b^{\alpha}_{\delta}\nabla^{\perp}h_{\delta}),\psi^{\alpha}_{\delta}\rangle_{K}\right.$
$\displaystyle\qquad-\left.\langle\hat{f}_{\omega^{\alpha}}((\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{+},(\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{-},\nabla^{\perp}\psi_{\delta}^{\alpha}\cdot{\hat{\mathbf{n}}}),{\psi^{\alpha}_{\delta}}^{-}\rangle_{\partial
K}\Big{\\}}\right.$ (6.35) $\displaystyle+\left.\Big{\\{}\underbrace{\langle
b^{\alpha}_{\delta}\nabla^{\perp}\psi^{\alpha}_{\delta},\nabla{h}_{\delta}\rangle_{K}}_{=-\langle
b^{\alpha}_{\delta}\nabla^{\perp}{h}_{\delta},\nabla\psi^{\alpha}_{\delta}\rangle_{K}}-\langle\hat{f}_{b^{\alpha}}(b_{\delta}^{\alpha,+},b_{\delta}^{\alpha,-},\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}),{h}_{\delta}^{-}\rangle_{\partial
K}\Big{\\}}\right]_{t=t_{0}}$
$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}E^{\alpha}_{\delta}(t)$
$\displaystyle=\sum_{K\in\mathcal{{\mathcal{D}}}_{\delta}}\Big{[}-\langle\hat{f}_{\omega^{\alpha}}((\omega^{\alpha}_{\delta}-b^{\alpha}_{\delta})^{+},(\omega_{\delta}^{\alpha}-b^{\alpha}_{\delta})^{-},\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}),{\psi^{\alpha}_{\delta}}^{-}\rangle_{\partial
K}$ (6.36)
$\displaystyle\qquad-\frac{1}{2}\langle\hat{f}_{b^{\alpha}}({b^{\alpha}_{\delta}}^{+},{b^{\alpha}_{\delta}}^{-},\nabla^{\perp}\psi^{\alpha}_{\delta}\cdot{\hat{\mathbf{n}}}),{h}_{\delta}^{-}\rangle_{\partial
K}\Big{]}_{t=t_{0}}$

Since $\psi_{\delta}$ and ${h}_{\delta}$ are continuous across element
boundaries, we have $\psi_{\delta}^{+}=\psi_{\delta}^{-}$ and
${h}_{\delta}^{+}={h}_{\delta}^{-}$. Thus, if the numerical fluxes
$\hat{f}_{\omega^{\alpha}}$ and $\hat{f}_{b^{\alpha}}$ satisfy the
conservative property (6.19), the sum of flux terms is zero. ∎

######  Remark 6.6.

In order for our semi-discrete discretisation to conserve numerical energy
when $\alpha>0$, more regularity is required for the approximation space for
$\psi^{\alpha}_{\delta}$. Additionally, the current scheme coupled with the
time stepping algorithm do not conserve casimirs. For future work, we look to
explore casimir conserving schemes for the model.

####  6.1.4. Time stepping

To discretise the time derivative, we use the strong stability preserving
Runge-Kutta of order 3 (SSPRK3) scheme, see [14]. Writing the finite element
spatial discretisation formally as

(6.37) $\partial_{t}b^{\alpha}_{\delta}={\rm f}_{\delta}(b^{\alpha}_{\delta})$

where ${\rm f}_{\delta}$ is the discretisation operator that follows from
(6.22), and

(6.38) $\partial_{t}\omega={\rm g}_{\delta}(\omega,b)$

where ${\rm g}_{\delta}$ is the discretisation operator that follows from
(6.23). Let $b^{\alpha,n}_{\delta}$ and $\omega^{\alpha,n}_{\delta}$ denote
the approximation of $b^{\alpha}_{\delta}$ and $\omega^{\alpha}_{\delta}$ at
time step $t_{n}$. The SSPRK3 time discretisation is as follows

(6.39a) $\displaystyle b^{(1)}$ $\displaystyle=b^{\alpha,n}_{\delta}+\Delta t\
{\rm f}_{\delta}(b^{\alpha,n}_{\delta})$ (6.39b) $\displaystyle\omega^{(1)}$
$\displaystyle=\omega^{\alpha,n}_{\delta}+\Delta t\ {\rm
g}_{\delta}\left(\omega^{\alpha,n}_{\delta},b^{\alpha,n}_{\delta}\right)$
(6.39c) $\displaystyle b^{\left(2\right)}$
$\displaystyle=\frac{3}{4}b^{\alpha,n}_{\delta}+\frac{1}{4}\left(b^{\left(1\right)}+\Delta
t\ {\rm f}_{\delta}\left(b^{\left(1\right)}\right)\right)$ (6.39d)
$\displaystyle\omega^{\left(2\right)}$
$\displaystyle=\frac{3}{4}\omega^{\alpha,n}_{\delta}+\frac{1}{4}\left(\omega^{\left(1\right)}+\Delta
t\ {\rm g}_{\delta}\left(\omega^{\left(1\right)},b^{(1)}\right)\right)$
(6.39e) $\displaystyle b^{\alpha,n+1}_{\delta}$
$\displaystyle=\frac{1}{3}b^{\alpha,n}_{\delta}+\frac{2}{3}\left(b^{\left(2\right)}+\Delta
t\ {\rm f}_{\delta}\left(b^{\left(2\right)}\right)\right)$ (6.39f)
$\displaystyle\omega^{\alpha,n+1}_{\delta}$
$\displaystyle=\frac{1}{3}\omega^{\alpha,n}_{\delta}+\frac{2}{3}\left(\omega^{\left(2\right)}+\Delta
t\ {\rm g}_{\delta}\left(\omega^{\left(2\right)},b^{(2)}\right)\right)$

where $\Delta t=t_{n+1}-t_{n}$ each $n$.

###  6.2. $\alpha$-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 $|{\bf
k}|^{-1}$, in the limit that $|{\bf k}|\to\infty$. Consequently, the TQG
dynamics is linearly well-posed. 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 $\alpha$-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 $|{\bf k}|^{-2}$, in the limit that $|{\bf k}|\to\infty$.
In the limit that $\alpha$ tends to zero, one recovers the TQG dispersion
relation derived in [17].

For the reader’s convenience, we repeat the equations solved by the pair of
variables $(b^{\alpha},\omega^{\alpha})$ in the $\alpha$-TQG system in
(2.12)–(2.13), as formulated now in vorticity-streamfunction form with fluid
velocity given by ${\mathbf{u}}^{\alpha}=\nabla^{\perp}\psi^{\alpha}$. This
formulation is given by

(6.40) $\displaystyle\frac{\partial}{\partial
t}b^{\alpha}+J(\psi^{\alpha},b^{\alpha})$ $\displaystyle=0,$ (6.41)
$\displaystyle\frac{\partial}{\partial
t}\omega^{\alpha}+J(\psi^{\alpha},\omega^{\alpha}-b^{\alpha})$
$\displaystyle=-\frac{1}{2}J(h,b^{\alpha}),$ (6.42)
$\displaystyle\omega^{\alpha}$
$\displaystyle=(\Delta-1)(1-\alpha\Delta)\psi^{\alpha}+f.$

Here, as before, $J(a,b)=\nabla^{\perp}a\cdot\nabla b=a_{x}b_{y}-b_{x}a_{y}$
denotes the Jacobian of any smooth functions $a$ and $b$ defined on the
$(x,y)$ plane.

Equilibrium states of the $\alpha$-TQG system in the equations (6.40)-(6.42)
satisfy
$J(\psi^{\alpha}_{e},\omega^{\alpha}_{e})=J(\psi^{\alpha}_{e},b^{\alpha}_{e})=J(h,b^{\alpha}_{e})=0.$
The equilibrium TQG state considered in [17] was given via the specification
of the following gradient fields

$\begin{array}[]{c}\nabla\psi^{\alpha}_{e}=-U{\hat{\mathbf{y}}},\
\nabla\omega^{\alpha}_{e}=(U-\beta){\hat{\mathbf{y}}},\ \nabla
f=-\beta{\hat{\mathbf{y}}}\\\ \\\ \nabla b^{\alpha}_{e}=-B{\hat{\mathbf{y}}},\
\nabla{h}=-H{\hat{\mathbf{y}}}\end{array}$

where we have taken $f=1-\beta y$ and the equilibrium parameters
$U,\beta,B,H\in{\mathbb{R}}$ are all constants. Linearising the $\alpha$-TQG
system around the steady state produces the following evolution equations for
the perturbations $\omega^{\prime}$, $b^{\prime}$ and ${\psi}^{\prime}$

(6.43a)
$\displaystyle\omega_{t}^{\prime}+U\omega^{\prime}_{x}+(U+B-\beta){\psi}_{x}^{\prime}$
$\displaystyle=(U-H/2)b^{\prime}_{x}$ (6.43b) $\displaystyle
b^{\prime}_{t}+Ub_{x}^{\prime}-B{\psi}_{x}^{\prime}$ $\displaystyle=0$ (6.43c)
$\displaystyle\Big{(}(1-\alpha\Delta)(1-\Delta)\Big{)}{\psi}^{\prime}$
$\displaystyle=-\,\omega^{\prime}.$

Since these equations are linear with constant coefficients, one has the plane
wave solutions

$\omega^{\prime}=e^{i({\bf k}\cdot\mathbf{x}-\nu t)}\hat{\omega},\quad
b^{\prime}=e^{i({\bf k}\cdot\mathbf{x}-\nu
t)}\hat{b},\quad{\psi}^{\prime}=e^{i({\bf k}\cdot\mathbf{x}-\nu
t)}\hat{{\psi}}.$

Further, from (6.43c) we obtain

$-\,\hat{\omega}=(|{\bf k}|^{2}+1)(\alpha|{\bf k}|^{2}+1)\hat{{\psi}}.$

Substituting these solutions into the linearised equations we have

$\displaystyle(\nu-kU)\hat{\omega}-k(U+B-\beta)\hat{{\psi}}$
$\displaystyle=-k(U-H/2)\hat{b}$ $\displaystyle\left(\nu-kU\right)\hat{b}$
$\displaystyle=-kB\hat{{\psi}}$

where $k$ is the first component of ${\bf k}$. From the above we obtain a
quadratic formula for the Doppler-shifted phase speed
$C=C(\alpha):=\left(\nu(\alpha)-kU\right)/k$

$C^{2}(|{\bf k}|^{2}+1)(\alpha|{\bf k}|^{2}+1)+CX+Y=0.$

where $X:=U+B-\beta$ and $Y:=(U-H/2)B$. Thus, the dispersion relation for the
thermal Rossby wave for TQG possesses two branches, corresponding to two
different phase velocities,

(6.44) $C(\alpha)=\frac{-X\pm\sqrt{X^{2}-4Y(|{\bf k}|^{2}+1)(\alpha|{\bf
k}|^{2}+1)}}{2(|{\bf k}|^{2}+1)(\alpha|{\bf k}|^{2}+1)}.$

Upon setting $\alpha=0$ in (6.44), we recover the Doppler-shifted phase speed
of thermal Rossby waves for the TQG system in [17]. When $\alpha=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 $|{\bf k}|\gg 1$; namely for $(|{\bf k}|^{2}+1)(\alpha|{\bf
k}|^{2}+1)\geq X^{2}/(4Y)$, for both branches of the dispersion relation.
However, the growth-rate of the instability found from the imaginary part of
$C(\alpha)$ in (6.44) decays to zero as $O(|{\bf k}|^{-2})$, for $|{\bf k}|\gg
1$. Indeed, the linear stability analysis leading to $C(\alpha)$ predicts that
the maximum growth rate occurs at a finite wavenumber $|{\bf k}|_{\max}$ which
depends on the value of $\alpha$, beyond which the growth rate of the
linearised TQG wave amplitude falls rapidly to zero. One would expect that
simulated solutions of $\alpha$-TQG would be most active at the length-scale
corresponding to $|{\bf 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.

Figure 3. An example of the growth-rate of linearised TQG waves as determined
from the imaginary part of $C(\alpha)$ in equation (6.44), plotted for
different values of $\alpha$. We observe that increasing $\alpha$ shifts the
wavenumber at the maximum growth-rate to lower values.

The magnitude of $\alpha$ controls the unstable growth rate of linearised TQG
waves at asymptotically high wavenumbers. That is, the presence of $\alpha$
regularises the wave activity at high wavenumbers, see Figure 3. Thus, from
the perspective of numerics, the $\alpha$ 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 ${\mathcal{D}}$ is an unit square with doubly periodic boundaries. We
chose to discretise ${\mathcal{D}}$ using a grid that consists of $256\times
256$ cells, i.e.
$\operatorname*{cardinality}({\mathcal{D}}_{\delta})=256\times 256$. This was
the maximum resolution we could computationally afford, to obtain results over
a reasonable amount of time.

We computed $\alpha$-TQG solutions for the following values of $\alpha$ –
$\frac{1}{16^{2}},\frac{1}{32^{2}},\frac{1}{64^{2}},\frac{1}{128^{2}},\frac{1}{180^{2}},\frac{1}{220^{2}},\frac{1}{256^{2}}$
and $0$. Note when $\alpha=0$ we get the TQG system. For all cases, the
following initial conditions were used,

(6.45) $\displaystyle\omega(0,x,y)$ $\displaystyle=\sin(8\pi x)\sin(8\pi