Surface Quasi-Geostrophy

Oceanic and atmospheric dynamics are often interpreted through potential vorticity, as this quantity is conserved along the geostrophic flow. However, in addition to potential vorticity, surface buoyancy is a conserved quantity, and this also affects the dynamics. Buoyancy at the ocean surface or at the atmospheric tropopause plays the same role of an active tracer as potential vorticity does since the velocity field can be deduced from these quantities. The surface quasi-geostrophic model has been proposed to explain the dynamics associated with surface buoyancy conservation and seems appealing for both the ocean and the atmosphere. In this review, we present its main characteristics in terms of coherent structures, instabilities and turbulent cascades. Furthermore, this model is mathematically studied for the possible formation of singularities, as it presents some analogies with three-dimensional Euler equations. Finally, we discuss its relevance for the ocean and the atmosphere.


Introduction
A large fraction of kinetic energy of oceanic flows is concentrated in quasi-horizontal motions of horizontal scales between 10 and 500 km and covering the upper 1500 m (the thermocline).The theory describing these motions mostly relies on the works of Charney [1] and is based on the conservation of quasi-geostrophic potential vorticity along the geostrophic flow.This theory applies to the ocean interior dynamics and has had numerous successes in understanding the energy exchanges between scales and the interaction between oceanic eddies (see, e.g., [2,3]).In the atmosphere, another picture emerges for which mid-latitude perturbations develop through the interactions of frontal structures at the tropopause and at the surface (the so-called Eady [4] model of baroclinic instability), which had some successes, for instance, in explaining the evolution of wave packets in storm tracks (e.g., [5]).Since the last 20 years, numerical simulations of the ocean at mid-latitudes unveiled the significant role of frontal structures at the ocean surface as was already remarked for the atmosphere.This was not foreseen because the scales of these oceanic fronts are of a few tens of kilometers (much below the Rossby radius of deformation) while they are of the order of 500 km in the atmosphere.Some studies [6,7] pointed out the possible use of the Surface Quasi-Geostrophy (SQG) model that was introduced by Blumen [8] to interpret ocean dynamics at these scales.This model has also attracted the community studying quasi-geostrophic turbulence, as it is a model that lays between purely two-dimensional (2D) barotropic and three-dimensional (3D) baroclinic dynamics (see Held et al. [9]).Its mathematical properties in terms of possible non-regularity of the solutions starting from smooth initial conditions were also put into evidence (see Constantin et al. [10]), in addition to its formal similarities with three-dimensional Euler equations.
In this review, we will discuss the characteristics of surface quasi-geostrophic flows.First, we will recall the general formulation of SQG, by introducing it in the quasi-geostrophic setting.This model is based on the conservation of an active scalar (surface buoyancy) along the horizontal geostrophic flow and on a particular relation between velocity and buoyancy.We will then discuss the existence of quasi-steady coherent structures and their related instabilities.Turbulent cascades in SQG will then be discussed, and the mathematical aspects of SQG will be covered in relation with the possible appearance of a singularity in finite time.The characteristics of passive tracers advected by SQG flows are finally discussed, as well as the connections between SQG and the ocean and atmosphere dynamics.

General Formulation of Surface Quasi-Geostrophy (SQG)
2.1.General Quasi-Geostrophic (QG) Theory The SQG model stems from a special case of the Quasi-Geostrophic (QG) equations.Quasi-geostrophy is the approximation of the equations of motion for a stratified fluid in rotation, with a small Rossby number (rotation of the Earth much larger than the rotation of the motion) and a small Froude number (high stratification).General theories for the QG system can be found in Pedlosky [2] and Vallis [3].We recall here the main points that are necessary for the presentation.
The salient feature of the QG system is the conservation of Potential Vorticity (PV) along the geostrophic flow trajectories.Potential vorticity is expressed as: in a Cartesian coordinate system (x, y, z) which is centered at a latitude φ 0 .Here, f 0 = 2Ω sin φ 0 is the Coriolis parameter at latitude φ 0 and Ω the rate of rotation of the Earth.The quantity ∂v/∂x − ∂u/∂y is the vertical component of the relative vorticity of the horizontal flow u = (u, v).In the ocean, buoyancy b is related to density ρ through: with g the gravity constant and ρ 0 a reference density.In the atmosphere, b is related to potential temperature θ through: with θ 0 a reference potential temperature.In both cases, N is the Brunt-Väisälä frequency (N 2 = ∂ z b with b a mean vertical buoyancy profile).
In the absence of forcing and dissipation, PV is conserved along the geostrophic flow: with ∇ the horizontal gradient.
Since the geostrophic velocity is divergent free, a streamfunction ψ can be introduced, so that: Relation (6) derives from the hydrostatic and geostrophic approximations.Relation (7) involves three terms: planetary and relative vorticities and a vortex stretching term (compression of column of fluid).In the following, we will solely examine PV anomalies, which consist of omitting the planetary vorticity term f 0 in the PV equation.
In addition to PV conservation, conditions on the horizontal and vertical boundaries of the domain are to be defined.In particular, surface buoyancy is conserved along the surface flow: where b s = b(x, y, z = 0) and u s = u(x, y, z = 0).Equation ( 8) is valid in the case of a rigid lid, that is to say, when vertical velocity w = 0 at z = 0.

Inversion of Potential Vorticity (PV) Equation
From Equation ( 7), it appears possible to determine the streamfunction field and other quantities (e.g., velocity, buoyancy, pressure fields) knowing the three-dimensional distribution of PV (at a given time).This is called the principle of PV inversion [11].Equation ( 7) is an elliptic equation if N 2 > 0 everywhere in the domain (condition of stable static stability).Therefore, solving this equation requires additional equations on the (lateral and vertical) boundaries of the domain.Setting aside lateral conditions (for instance, for a doubly-periodic domain) and imagining a semi-infinite domain in the vertical, the remaining condition is at the surface z = 0, i.e., Equation (9).Therefore, both PV in the fluid interior and surface buoyancy are needed in determining the three-dimensional streamfunction field.
Due to the particular form of the surface condition and the elliptic operator, Bretherton [12] observed that, in the QG approximation, buoyancy at the surface plays the same role as potential vorticity in the interior of the fluid.This can be analytically established (see [7,12]), and the surface condition: can be replaced by: provided that Equation (7) becomes: Hence, the surface buoyancy plays the role of a Dirac function for the PV field (i.e., it can be considered as a "PV sheet").Note that this idea has been extended in the context of the primitive equations by Schneider et al. [13].The analogy between interior PV and surface buoyancy is also apparent when considering baroclinic instability criteria [14], as necessary conditions of instability involve the change in sign of either the interior potential vorticity gradient or the opposite signs between the interior PV gradient and the surface buoyancy gradient.
From Equation (12), the contributions of the surface buoyancy and the interior PV to the total flow can be separated leading to interior and a surface-induced dynamics [7,11].An interior-induced dynamics corresponds to the situation of conservation of the interior PV along the geostrophic flow with uniform buoyancy at the surface: This is the classical model of the quasi-geostrophic theory [1], and its characteristics are relatively well known (see, e.g., [2,3,15,16]).
A surface-induced dynamics corresponds to the situation of conservation of the surface buoyancy along the surface geostrophic flow with uniform PV in the interior: This model has been considered, among others, by Blumen [8] and Held et al. [9], who first studied its characteristics.We will now present the different results that are known about this system.

Surface Quasi-Geostrophy (SQG) Formulation
The surface quasi-geostrophic equations in their non-dimensional form can be expressed as: To obtain these non-dimensional equations, a constant Brunt-Väisälä frequency N 2 has been assumed, and the vertical coordinate has been transformed using zN/ f 0 as a new vertical coordinate.The non-dimensional variable θ corresponds to the surface buoyancy b s .Here, we choose to place ourselves in the atmospheric case (considering a semi-infinite domain with z > 0).
From Equations ( 19)-( 21), the following relation is obtained in the horizontal Fourier domain: with K = | k| the wavenumber modulus.ψ( k, z) is the two-dimensional Fourier transform of the streamfunction ψ at the altitude z, and θ( k) is the Fourier transform of the surface buoyancy.Equation (24) shows that, for each Fourier component, the streamfunction decreases exponentially with z.Observe in particular that the decrease is even faster when small horizontal scales are considered (large K).

Invariants
There are two invariants in this system: the surface potential energy (often called generalized enstrophy by analogy to relative enstrophy for barotropic flow [17]), and the total vertically-integrated (kinetic and potential) energy (often called generalized energy [17]), with ψ s = ψ(x, y, z = 0).The last equality in Equation ( 26) comes from the uniformity of PV, i.e., Equation (19).From Equation ( 24), we deduce that surface kinetic energy is proportional to surface potential energy: with k = (k, l) the horizontal wavenumber.

Relation between the Active Tracer and Streamfunction
The buoyancy θ is an active tracer since it is advected by the flow it creates.This is similar to the nature of relative vorticity of a two-dimensional barotropic flow.The expression relating buoyancy to the streamfunction is: while the relation relating relative vorticity ζ = ∂v/∂x − ∂u/∂y = ∂ 2 ψ/∂x 2 + ∂ 2 ψ/∂y 2 to streamfunction is: Such analogies persist in the physical space, taking θ as a Dirac density function θ( x) = δ( x − x 0 ): In the SQG case, the velocity modulus decreases as 1/| x − x 0 | 2 , much sharper than 1/| x − x 0 | for the barotropic case [9].This means that barotropic vortices have a farther range influence compared to their SQG counterparts.

Exact Surface Quasi-Geostrophy (SQG) Solutions
Several exact solutions of the SQG equations that are steady in some Galilean reference frame have been found.A first class of solutions is shear lines, i.e., buoyancy being a function of a single dimension (e.g., θ = f (y)).Rectilinear strips with varying or uniform buoyancy over a finite band in the transverse direction belong to this class [18].
Concerning mono-polar vortices, circular vortices are evidently steady solutions of the SQG equations.Dritschel [19] presented an SQG solution corresponding to a steadily rotating ellipse with nonuniform buoyancy in the interior of the ellipse.The possibility of the existence of a rotating ellipse with uniform buoyancy was treated by Castro et al. [20], who showed that such an ellipse would not persist in time while keeping its shape.Nonetheless, Castro et al. [21] proved the existence of convex global rotating vortices with C ∞ regularity for their boundary curves and with uniform buoyancy in their interior, but their proof does not allow for any explicit solutions.
The case of two vortices was treated by several authors.A first simple model consists of two point vortices of the same or different buoyancies.Lim and Majda [22] showed that, in this case, the distance between the two partners remains fixed in time, and they form curvilinear trajectories.In the case of the same buoyancy anomaly, they rotate around each other, while they steadily translate in a straight line when their buoyancy is opposite.For continuous buoyancy profiles, the analytical solution of such steadily-translating vortex dipoles was described by Muraki and Snyder [23].Carton et al. [24] presented the solution of two co-rotating vortices with uniform buoyancy in the presence or absence of a large-scale strain field.This last study also examined analytically and numerically the merger of these vortices in the presence or absence of an external deformation field.It was found that the distance of vortex merging was generally smaller than the two-dimensional barotropic case [24].This can be explained by the rapid decrease of the velocity as a function of distance (see Section 2.5), which prevents the vortices from exerting an influence on each other at a large distance.

Contour Dynamics
Vortices formed by piecewise-constant buoyancy within finite-area regions may be seen as building block models to study vortex dynamics since only the time evolution of vortex boundaries needs to be known.This is the essence of the approach of contour dynamics [25].However, contrary to the barotropic case, the equation of evolution of the interface for SQG flows requires a careful treatment.Indeed, the tangential velocity along the contour is divergent, while the local contribution to the normal velocity (displacing the material curve) vanishes at leading order, allowing some regularity in the contour equations [9,18].
In this context, Cordoba et al. [26] determined the evolution equation for the interface of an "almost sharp front", that is to say, a strip of finite width δ separating two regions of constant buoyancy, δ being a small parameter.They showed that the equation with δ tending to zero corresponds to a weak solution of the equation of motion, and Fefferman and Rodrigo [27] were able to construct analytic solutions for this problem.Rodrigo [28] examined the case of sharp fronts (discontinuities in buoyancy) and showed that, if θ(x, y, t) is a weak solution of the equations, then it is possible to obtain an equation for the interface.This result was extended by Gancedo [29], and local-in-time existence and uniqueness for smooth contours were obtained.Cordoba et al. [30] examined numerical simulations of contour dynamics and showed an example of the collapse of two iso-buoyancy contours at a single point.In the case of the instability of an ellipse with uniform buoyancy, Scott and Dritschel [31] also observed the collapse of a filament width in finite time.However, as shown by Gancedo and Strain [32], such a collapse can only occur if the curvature of the contour tends to infinity in finite time, which was the case in the studies of Cordoba et al. [30] and Scott and Dritschel [31].These different examples do not contradict the mathematical results of Rodrigo [28] and Gancedo [29] because their theorems rely on the smoothness of the vortex boundary, which is not the case here.

Instabilities
The instability of strips of uniform buoyancy was studied by Juckes [18], who showed that the number of the most unstable waves, as well as their rate of growth are inversely proportional to the width of the strip.This is due to the conservation of buoyancy that implies that the relative vorticity ζ increases when the filament becomes thinner since, by Equation ( 28), the vorticity scales as θ/L, L being the width of the filament.This leads to the typical roll-up of the filament, as in the case of barotropic flows, but with a much smaller timescale [9].The general mechanism to inhibit the instability of a filament is the presence of a background strain, for instance due to a distant eddy [33].However, due to the fast decay of the velocity field from a point source, which follows from Equation (30), this mechanism will be less effective than in the barotropic case.Indeed, Harvey and Ambaum [34] showed that, contrary to the barotropic case, the large-scale deformation or adverse shear is not capable of stabilizing these filaments because the thinning of the filament due to these processes increases the growth rate of the instability.
The fact that the instability of a filament is more rapid as its width decreases with the subsequent formation of new vortices of smaller scales may lead to a cascade of instabilities, as these vortices will also stir new filaments at a thinner and thinner scale.This cascade of filament instabilities may lead to a finite time singularity as discussed by Scott [35] and Scott and Dritschel [31].This is well illustrated in Figure 1, which represents the buoyancy contours at a time close to the formation of the singularity.The scenario of the instability can be described as follows: the first filament developed an instability with a time scale that is inversely proportional to its width.Vortices at the scale of the instability (i.e., the filament width) are then formed and stir a filament between them.As explained above, the shear induced by the vortices is not sufficient to inhibit the instability of the new filament.A new instability takes place, which is more rapid and at a smaller scale, since the filament is made thinner by the straining of the nearby vortices.This process continues more and more rapidly, leading to a self-similar cascade of instabilities.The final evolution consists of the collapse of the filament width to zero and to its infinite curvature.Concerning coherent vortices, Carton [36] numerically examined the instability of circular vortices with different buoyancy profiles.For moderate buoyancy gradients at the edge of the vortex, the linear instability is similar to the classical barotropic instability, while it is more unstable for stronger gradients.In addition, the nonlinear equilibration proceeds as in the barotropic case with the formation of multipoles, such as tripoles [36].Moreover, Harvey and Ambaum [37] examined the instability of a circular vortex with uniform buoyancy (a Rankine vortex) and showed that it was also similar to the instability of the barotropic circular vortex.Other classes of axisymmetric vortices were also studied, such as shielded vortices (vortices with zero net buoyancy) with a constant-buoyancy core surrounded by a constant-buoyancy annulus [38].This last study confirmed the numerical results of Carton [36].In addition to axisymmetric vortices, one can examine the instability of the steadily-rotating ellipse documented by Dritschel [19].It was found to be stable to perturbations for a smaller range of aspect ratios compared to the two-dimensional barotropic case.The instability leads to the ejection of filaments as in the barotropic case, but these filaments develop shear line instabilities, rolling up into small vortices, which was not observed in the barotropic case.More general studies on the properties of instabilities of SQG flows have been done.Friedlander and Shvydkoy [39] examined the unstable spectrum of solutions (eigenvectors of the linear problem) for a particular spatially oscillating shear flow.In such a case, the general solutions (the unstable spectrum) were also found to share the same properties as the standard 2D barotropic case.These different examples illustrate that barotropic and SQG flows bear some common properties in terms of linear instability.
Nonlinear instability was treated by Friedlander et al. [40] for the SQG equations in the presence of a constant forcing in time and with a dissipation operator.It was found that, like barotropic flows, the linear instability of stationary solutions implies nonlinear instability (à la Lyapunov).A similar problem about the nonlinear growth of small perturbations in the context of atmospheric predictability was considered by Rotunno and Snyder [41].It was found that the SQG model exhibits limited predictability, unlike the barotropic model, as the initial small-scale errors rapidly grow and saturate, and a slower growth is then observed when these perturbations reach larger scales.On the contrary, in the barotropic case, the errors at small scales cascade towards errors at larger scales, but still having small amplitudes, and then, saturate relatively slowly in time.A last result that can be discussed concerns the number of degrees of freedom as defined as the number of positive local Lyapunov exponents (corresponding to an orthogonal basis of perturbations with an instantaneous positive growth rate).Tran et al. [42] showed that this number grows at most in Re 3/2 with Re the Reynolds number, a result that is in agreement with the phenomenology of turbulent cascades.

Turbulent Cascades
Surface quasi-geostrophy possesses particular properties in term of turbulent cascades as it conjugates two-dimensionality through the advection of buoyancy and possible transfers between potential and kinetic energies.

Buoyancy Variance Spectra
Blumen [8] showed that, because of the existence of two invariants (total energy E and potential energy P), two turbulent regimes can be identified: an inverse cascade of total energy associated with a spectrum of buoyancy variance in k −1 ; a direct cascade of potential energy with a spectrum of buoyancy variance in k −5/3 .This has to be contrasted with the barotropic QG system for which the direct cascade of enstrophy gives an enstrophy spectrum steeper than k −1 , and in the inverse kinetic energy, the enstrophy spectrum is in k 1/3 , corresponding to a kinetic energy spectrum in k −5/3 [43].Using an Eddy-Damped Quasi-Normal Markovian (EDQNM) model, Hoyer and Sadourny [44] confirmed the scenario of the direct and inverse cascades and the corresponding spectra.(see also [45,46] for a discussion of the inverse cascade for a particular large-scale dissipation operator.)Different numerical simulations were done, in the presence or absence of forcing or large-scale dissipation.Pierrehumbert et al. [47] and Celani et al. [48] examined the direct cascade of buoyancy variance and found that the spectral slope was slightly steeper than k −5/3 in the case of large-scale forcing.Sukhatme and Pierrehumbert [49] and Capet et al. [50] also found steeper slopes in freely-decaying simulations during the time when an inertial range can be observed.These studies indicate that a correction to the standard similarity argument is needed.Watanabe and Iwayama [51] reconsidered the standard similarity hypothesis of Kraichnan [43] and obtained a prediction for the buoyancy variance slope shallower than k −5/3 , which they observe in their simulation where they force and dissipate the energy at large scales.As discussed by Tran and Bowman [45], Tran [46] and Constantin [52], a particular form of the dissipation operator can lead to different upper bounds for the spectra.The discrepancy of [51] with other results may thus be due to the types of forcing and dissipation used in each numerical setting (see Burgess et al. [53] for a discussion of such an effect).It may also depend on the presence of vortices that generally steepen the spectra.
For the inverse cascade, Smith et al. [17] and Tobias and Cattaneo [54] observed a spectrum consistent with k −1 for the buoyancy variance in the presence of forcing and dissipation.Burgess et al. [53] observed a steeper spectrum when dissipation is absent, which they attribute to the formation of coherent vortices.

Physical Properties of the Cascade
As can be seen in Figure 2a, SQG turbulence is characterized by the presence of vortices at all scales because instabilities of fine-scale filaments are weakly inhibited by larger-scale vortices, leading to new instabilities [31].This suggests a possible multifractal nature of the turbulent field, as noted by Sukhatme and Pierrehumbert [49].For a simulation of decaying turbulence (hence, in the direct cascade of buoyancy variance), they observed that the buoyancy field becomes rough in a particular range of scales.Moreover, in the SQG inverse cascade, Bernard et al. [55] showed that the fractal dimension of the buoyancy field is d ≈ 1.5, as for the relative vorticity field in the direct cascade of 2D barotropic turbulence, This study also showed that the zero-buoyancy isolines are conformally invariant as in the 2D inverse cascade of barotropic turbulence.This result indicates that buoyancy isolines are equivalent to curves that can be mapped into a one-dimensional Brownian walk.Concerning now the properties of the buoyancy field at large time, Venaille et al. [56] showed that the formation of a sharp interface between two regions of homogenized buoyancy will occur as time tends to infinity due to a coarsening process.The statistical equilibrium generally corresponds to a one-dimensional flow field and not a large-scale dipole that would occur for the inverse cascade of barotropic turbulence.The evolution to statistical equilibria can also be examined in terms of spectral properties.In this context, Teitelbaum and Mininni [57] showed that generalized enstrophy thermalizes (the spectra of buoyancy variance behaving as k for large k), while generalized energy condenses (piling up at the largest scale).

Role of Meridional Gradients and Linear Damping
One can examine the more general problem of an SQG system with an additional meridional gradient of buoyancy and a linear damping term: The last term on the left-hand side of Equation ( 32) is the meridional advection of a large-scale buoyancy θ = βy, while the term on the right-hand side exerts a "thermal" damping on buoyancy.When β is larger than the advective terms, the solutions of this equation are the equivalent to Rossby waves of barotropic turbulence.If we continue with the analogy with barotropic turbulence [58], when β is as large as the nonlinear advective terms, these waves will interact with the nonlinearities of the flow, and zonal jets will develop.This is illustrated by Figure 3a, which shows the total buoyancy field θ + βy in a numerical simulation resolving Equation (32).These jets lie between regions of homogenized buoyancy.In between, small-scale vortices interact between each other without cascading to larger scales (Figure 3b).In the inverse cascade regime, the kinetic energy spectrum should follow a k −3 dependency along the meridional direction [17].For comparison, in the inverse energy cascade of barotropic turbulence, one would have instead an anisotropic spectrum in k −5 .Sukhatme and Smith [59] were not able to see the k −3 regime in their numerical simulation resolving the inverse SQG cascade.They instead noted shallower spectra for increasing β in the inertial range domain.In the direct cascade, the buoyancy variance spectra were observed to be mostly insensitive to the value of β.
The effect of the linear damping term in absence of β was examined by Smith et al. [17].This study showed that the inverse cascade should be halted at a particular wavenumber k r ∼ rg −1/3 with g the energy injection rate to the system, and the generalized energy spectrum should behave as 2 .Numerical simulations were not very conclusive about the shape of the energy spectrum, but the arrest scale was well predicted.

Energy Transfers
The prediction of the buoyancy variance spectrum in k −5/3 is obtained if one assumes that only local scales can affect the buoyancy variance spectral transfers between scales.In barotropic turbulence, this is not the case, and a steeper spectrum of enstrophy than predicted was observed by numerous studies (see, e.g., Vallis [3]).This can be measured through the study of the contributions of local and nonlocal scale interactions in the buoyancy variance flux at a given scale.Watanabe and Iwayama [60] showed that both local and nonlocal triad interactions contribute to the total transfer.This is contrary to the standard expectation based on a strain-rate variance σ 2 (k) = k 0 K 2 Q(K)dK scaling as k 4/3 for a buoyancy variance spectrum Q(K) in k −5/3 [47].Such a scaling would indicate that local scales dominate the time scale of the spectral fluxes.Watanabe and Iwayama [51] presented a correction of the k −5/3 spectrum that was in fair agreement with the result of their numerical simulation.The issue of locality or nonlocality in the spectral transfers for the direct cascade is not yet totally solved, as the flatter buoyancy spectrum observed by Watanabe and Iwayama [51] is in disagreement with other studies, which show steeper spectra.
The SQG system possesses at the same time surface potential energy P = θ 2 /2 and surface kinetic energy E = | u s | 2 /2, and one might expect some transfers between the two as is the case for three-dimensional QG turbulent flows.The conservation equation of buoyancy variance gives: where () * stands for the complex conjugate and [] the real part.In the absence of dissipation, the interaction among triads transfers buoyancy variance to small scales [60].By manipulating the momentum equations, the conservation of kinetic energy in spectral space can now be written as: where w s z = ∂w/∂z at z = 0 is the opposite of the surface divergence [50].The first term on the right-hand side of this equation corresponds to the transfer of kinetic energy among wavenumbers.The second term corresponds to the effect of the ageostrophic pressure that allows the transfer of potential energy into kinetic energy within the fluid.Introducing the spectral transfer functions: where is an average over time and spectral shells between K and K + dK and using (28), which implies that the buoyancy variance is equal to the kinetic energy for each wavenumber, we obtain: in absence of source or dissipation.In the case of the direct cascade to small scales, there exists an inertial range over which Π θ is constant and positive.As shown by Capet et al. [50], in general, a transfer of potential to kinetic energy can be expected at all scales, so that Π a should be positive.This implies that Π u < 0, which means that in the direct cascade of buoyancy variance, the nonlinear interactions promote an inverse cascade of surface kinetic energy.The different fluxes are shown in Figure 2b for a simulation in free-decay in a regime where the buoyancy variance flux Π θ shows the presence of an inertial range.The signs of Π a and Π u are in agreement with the argument explained above [50].It is interesting to note that the inverse kinetic energy transfer was observed in the ocean using satellite observations [61], suggesting that the SQG turbulence could be a possible interpretation to it.

Surface Quasi-Geostrophy (SQG) from a Mathematical Point of View
The SQG system has attracted the mathematician community due to its formal proximity with the three-dimensional Euler equations.It thus appeared that if some results could be obtained about the existence of finite-time singularity from smoothed initial conditions in SQG, the methods employed might be a good start to tackle the 3D problem, as well [10,62].

Analogy with 3D Euler Equations
Constantin et al. [10] pointed out that SQG equations and the three-dimensional Euler equations share some common properties.Indeed, the time evolution of the horizontal buoyancy gradient in SQG is expressed as: with ∇ ⊥ θ = (−∂θ/∂y, ∂θ/∂x), D/Dt the Lagrangian derivative and ∇ the horizontal gradient.Moreover, the velocity field can be expressed as: using Equation (30).In this case, the velocity gradient tensor can be expressed as a function of buoyancy: in the horizontal Fourier space with k the two-dimensional horizontal wavenumber.The operator ⊗ is the matrix formed by the tensor product of two vectors [62].
On the other hand, for the 3D Euler equations, the time evolution equation for the vorticity vector ω is expressed as: with ∇ 3 the three-dimensional gradient.The 3D velocity v field takes the form: and the velocity gradient tensor can be expressed as a function of the vorticity vector: in spectral three-dimensional space with k the three-dimensional wavenumber.
In both systems, the vectors representing the dynamics (∇ ⊥ θ and ω) are stretched by a velocity gradient that linearly depends on them.Hence, the time evolution induced by this nonlinearity may lead to a rapid increase of the 3D vorticity vector and the buoyancy gradient norms with the possible blow-up of the solutions [10,62].

Regularity of Surface Quasi-Geostrophy (SQG) Solutions
For two-dimensional Euler equations, no singularity can occur in finite time in the case of smooth initial conditions.Only the vorticity gradient can increase without bound as time goes to infinity (see, e.g., [63][64][65] for a demonstration and related review on earlier work).This result was then extended by Dutton [66], Bennett and Kloeden [67], Chae [68] for three-dimensional quasi-geostrophic flows, but restricted to the case with no surface buoyancy.The SQG case was not investigated until the work of Constantin et al. [10].
For the three-dimensional Euler equations, necessary conditions for finite-time singularity can be derived [69] and provide some help in understanding the regularity of SQG solutions.Consider an initial condition v 0 belonging to the Sobolev space H s (R 3 ) whose distributional derivatives up to order s (with s ≥ 3) are in L 2 ,s (R 3 ).A smooth solution of the 3D Euler equations v cannot be continued to t = T if and only if the vorticity vector satisfies: where The regularity in the velocity field results in the finite growth of the vorticity vector.The analogy between the vorticity vector and the buoyancy gradient discussed in the previous section still goes on with a similar criterion for SQG [10].In this case, the singularity can occur at t = T if and only if: where The singularity is thus associated with an infinite buoyancy gradient in finite time.Constantin et al. [10] proposed an initial condition (a zonal shear oscillating in the x direction) that was a candidate to singularity development (see also [62]).However, studies using high-resolution numerical models [70] showed that this was not the case.Instead, the buoyancy gradient norm |∇θ| exhibits a growth rate in exp(exp(σt)).Such a double exponential can be understood when examining Equations ( 39) and ( 41) which indicate that the strain tensor that applies to the buoyancy gradient is proportional to this gradient in return.A mathematical result by Denisov [71] shows that such a double exponential law can also occur for the vorticity gradient in standard 2D Euler equations.The absence of the finite time singularity for the initial solution considered in Constantin et al. [10] was confirmed by Cordoba [72,73] and Constantin et al. [74], who derived necessary conditions based on the curvature of the saddle-point of the flow.The question of the finite-time singularity raised by Constantin et al. [10] motivated different studies (e.g., [72,75,76], among others) who examined the regularity of solutions of the SQG system in the more general setting of forced-dissipative systems in a periodic domain: with f a stationary or time-varying forcing, ∆ the 2D Laplacian operator and with a viscosity parameter κ > 0. The exponent γ determines different kinds of dissipation, the most interesting one being γ = 1, which is similar to an Ekman friction.Constantin et al. [10] showed that, in absence of forcing and for κ = 0, smooth solutions (starting from smooth initial conditions) exist over finite time, while Resnick [77] showed that weak solutions exist globally, but might not be unique.This was further studied by Constantin and Wu [78], who showed that, for 1 < γ ≤ 2, smooth solutions exist for all time in the presence of a smooth and time-varying forcing.They also showed that weak solutions must coincide with strong solutions as long as these latter exist.For 0 < γ < 1, the existence and uniqueness of weak solutions in the absence of forcing was obtained by Ju [76,79] and Wu [80], while Chae and Lee [81] and Cordoba and Cordoba [82] considered the case of smooth solutions.The case γ = 1 was first considered by Constantin et al. [83], who proved the global existence of solutions for small initial solutions.Later, Kiselev et al. [84], Caffarelli and Vasseur [85], Constantin et al. [86] extended this result, removing the condition of smallness and in the presence of a sufficiently smooth forcing.

Passive Tracers
The behavior of a tracer passively advected by the flow was examined by different authors.Consider a tracer C that obeys the conservation equation: where u g is the horizontal velocity field at some depth z that can be obtained through Equation (24).At the surface z = 0, the velocity field is related to a turbulent regime with a well-defined inertial range (either in the direct cascade of generalized enstrophy or the inverse cascade of generalized energy).Using the standard assumption of self-similar cascades, the spectrum of the tracer variance can be predicted to follow a law in k −2 in the inverse cascade and in k −5/3 in the direct cascade [17,87].Smith et al. [17] confirmed in their numerical simulation the k −2 regime in the inverse cascade for the tracer variance.However, simulations in the direct cascade regime show that the tracer is shallower than k −5/3 [48,88].A tentative explanation would rely on the steeper spectrum of kinetic energy (in k −2 ) compared to the theory (in k −5/3 ) so that both large-scale and small-scale components of the flow contribute to the stirring of the tracer.One can imagine that the large-scale component represents an important contribution to this stirring for a steeper spectral slope of kinetic energy, decoupling the stirring scale to the tracer scale.In that situation, the theory of Batchelor [89] would explain a spectrum in k −1 .
If, now, we consider an altitude z > 0, from (24), we observe that the velocity field decays almost exponentially for sufficiently large wavenumbers.Introducing a transitional wavenumber k c ∼ 1/2z, the kinetic energy spectrum decays as E(k, z) = E(k, z = 0) exp(−k/k c ) with E(k, z = 0) its surface value.Two limiting cases are obtained: for k k c , the kinetic energy spectrum will resemble E(k, z = 0), while for k k c , it decays exponentially with k.Using this argument, Scott [88] showed that, in the direct cascade, for k k c , the spectrum should follow a k −5/3 power law, while for k k c , the spectrum should be in k −1 .In this latter case, only the large-scale component of the strain field stirs the tracer following the argument of Batchelor [89].Scott [88] observed a good agreement between these predictions and a numerical simulation of SQG turbulence in the direct cascade of surface buoyancy variance.This is also confirmed by the study of Wirth et al. [90], who examined the advection of a passive tracer at different altitudes z in the case of the merging of two surface vortices.They observed that the tracer in the interior of the fluid was stirred in numerous filaments as would occur for a scale separation between the stirring and the tracer fields [89].
A final remark concerns the passive tracer spectra at a particular altitude z in the inverse cascade.If we use the same argument as Scott [88], we obtain a tracer variance spectrum in k −2 for k k c and a dependence in k −1 for k k c .While the regime k k c corresponds to what would occur at the surface, the regime k k c still represents a scale separation between the stirring and the passive tracer fields.Such a result was not investigated to our knowledge.

Horizontal Motions from Surface Quasi-Geostrophy (SQG) Equilibrium
The ocean surface is characterized by mesoscale structures (eddies and fronts of scales of 50-300 km, corresponding to the Rossby radius of deformation) and sub-mesoscale structures (fronts and filaments of 10 km in width and hundreds of kilometers in length).These filaments are associated with horizontal gradients of buoyancy, temperature, salinity and other tracers.Since the early 2000s, these sub-mesoscales were assumed to play a passive role in the ocean dynamics, as they were supposed as weakly energetic.More recently, numerical high-resolution simulations resolving these structures unveiled the importance of these small scales due to the high vertical velocities associated with them, their large values of relative vorticity and their strong contribution to the vertical heat flux [91,92].This was further confirmed during the "Scalable Lateral Mixing and Coherent Turbulence" (LatMix) experiment that took place in the Sargasso Sea [93].However, the standard QG theory of Charney [1] could not apply to understanding this dynamics due to different reasons: it does not take into account the surface buoyancy gradient; it considers a small Rossby number and fixed static stability.As remarked by Lapeyre and Klein [7], oceanic surface buoyancy plays an analogous role to the potential vorticity within the ocean (see Section 2.2).Because of this, the upper layers of the ocean (the first 500 m) have a very different dynamics than the interior layers.They are governed by the dynamics related to these small scales associated with strong buoyancy gradients, and the SQG model was revealed to be a good starting point to understand their properties.The potential of this model for the surface ocean dynamics has probably been first put forward by Lesieur and Sadourny [87], who proposed to explain the observed phytoplankton spectra by SQG theory.One can also mention the work of Johnson [94], who proposed SQG solutions for deep vortices trapped above topographic bumps.Furthermore, Held et al. [9], in their fairly detailed study of the SQG model, remarked at the end of their conclusion that this model could be useful for the ocean, but without specifying in what aspect.A clear link between SQG and the upper ocean dynamics was made by LaCasce and Mahadevan [6] and independently by Lapeyre and Klein [7].
Lapeyre and Klein [7] and Klein et al. [91] examined whether the relation between streamfunction and buoyancy Equation ( 24) was observed in a primitive equation simulation of an idealized baroclinically unstable front.More precisely, in dimensional units, the balance ( 24) is simply: with b buoyancy defined by Equation (2).Comparing relative vorticity ζ = ∂v/∂x − ∂u/∂y with its prediction using ζ = ∇ 2 ψ and relation (49), they confirmed that the surface of the ocean was in an SQG balance.Other situations were examined, such as a realistic simulation of the North Atlantic [95] or the North Pacific [96].In Isern-Fontanet et al. [95], it was shown that the buoyancy at the base of the oceanic mixed layer was a better proxy for b(x, y, z = 0).Moreover, using instead sea surface temperature showed good results in reconstructing relative vorticity from buoyancy.These different studies observed that the SQG relation (49) was satisfied down to typical depths of 500 m and horizontal scales between 300 and 50 km [95,97].This situation occurs only if the contribution of the surface signal on the dynamics is dominant, which may not be the case in all regions of the ocean [98].The use of the SQG equilibrium together with satellite datasets opens the possibility to access the sub-mesoscale features from sea surface temperature images or sea surface heights obtained by altimetry.Isern-Fontanet et al. [99] indeed confirmed that geostrophic velocities obtained through sea surface temperature using the SQG balance compare well with the geostrophic velocities from altimetry in the Gulf Stream region (see Figure 4).This was also confirmed on a more global scale by Gonzalez-Haro and Isern-Fontanet [100].Note the good agreement between the two velocity fields.From [99], with permission of John Wiley and Sons.

Relation between Surface and Interior Dynamics
Usually, the decomposition of oceanic motions involves barotropic and baroclinic modes [1,101].These modes are orthogonal for the total (potential and kinetic) energy norm and are the solution of an eigenvalue problem involving the diagonalization of potential vorticity.However, Lapeyre [98] showed that this problem does not include the surface oceanic condition (surface buoyancy); hence, the decomposition in baroclinic modes is not complete, as a mode corresponding to uniform potential vorticity (satisfying the SQG equilibrium) is needed.As a result, the new basis is not orthogonal anymore, and the barotropic and the first baroclinic modes have a strong projection on the "SQG" mode.A consequence of this is that the sea surface height measured by altimetry is not the signal of the first baroclinic mode (as suggested by Wunsch [101], Stammer [102], among others), but an SQG mode trapped in the top 500 m of the ocean [98].Another consequence is that the isotropization of energy as described by Charney [1] may not take place if surface buoyancy is non-uniform.
To circumvent the non-completeness of standard baroclinic modes, Smith and Vanneste [103] introduced a new "surface-aware" basis of eigenmodes taking into account the surface boundary condition.Their new modes are now dependent on the wavenumber, contrary to the standard baroclinic modes, which only depend on the stratification and the Coriolis parameter.A drawback of the new basis introduced by Smith and Vanneste [103] is that it necessitates the choice of a nondimensional parameter that seems to depend on the particular situation under study.
Roullet et al. [104] confirmed the paradigm of Lapeyre and Klein [7] as they put in evidence the separation between a surface dynamics controlled by surface buoyancy anomalies and an interior dynamics controlled by interior PV.A question then arises about a possible correlation between the contributions of the surface buoyancy and the interior potential vorticity to the dynamics.Lapeyre and Klein [7] and Klein et al. [97] showed that due to the meridional distribution of the density field, the large-scale potential vorticity gradients are often similar to the large-scale surface buoyancy gradients.As a result, it is possible to show that PV anomalies at mesoscales are proportional to surface buoyancy anomalies [7,105].The contribution due to the interior potential vorticity to the dynamics is then highly correlated with the surface contribution, which may explain why the SQG equilibrium seems to be robust in the ocean.
New models including the interactions of a few baroclinic modes (representing PV anomalies) and the SQG modes were developed to better understand how interior and surface dynamics interact [106,107].These ideas served in developing techniques to reconstruct the three-dimensional velocity field from the sea surface temperature and sea surface height [105,108,109].One can also mention that SQG models served to develop idealized models for the oceanic mixed layer dynamics [110,111].These last studies used a system of uniform potential vorticity within two layers of different depth.For a top layer with small depth, one recovers some basic properties of the mixed layer instabilities observed in more realistic flows by Boccaletti et al. [112].

Vertical Motions from Surface Quasi-Geostrophy (SQG) Equilibrium
From the general quasi-geostrophic equations, it is possible to compute vertical velocities by solving the so-called omega-equation: where [∇ u] is the horizontal velocity gradient tensor and [ ] T the matrix transpose operator [113].
To obtain w from this equation, an elliptic problem needs to be solved with w = 0 at the surface.For the SQG model, a simpler form exists in spectral space: where the subscript s is for surface quantities.This last relation in addition to Equations ( 6) and (49) shows that vertical velocities are completely determined by surface quantities.Such an equation was used to diagnose vertical velocities from surface buoyancy in different numerical simulations with some success [7,96,114].Ponte et al. [115] improved this estimation when taking into account wind-driven Ekman motions and a diabatic contribution due to surface mixing as parametrized by Garrett and Loder [116].

Frontogenesis
An important and relatively little considered feature of SQG (except [117] and the semi-geostrophic studies that considered cases of uniform PV) is its direct link with frontogenesis, i.e., the production process of horizontal gradients of surface buoyancy by stretching, rotation and deformation [118,119].As was shown by Sawyer [120], Eliassen [121] and Hoskins and Bretherton [122], this frontogenesis has significant impact on the three-dimensional dynamics, in particular vertical velocities as illustrated by Equation (50).The equation of horizontal buoyancy gradient at some depth z is: with: The vector Q was introduced by Hoskins and Bretherton [122] and quantifies the effect of deformation and stretching of the buoyancy contours, which enhances buoyancy gradients [118,123].The growth of buoyancy gradients by straining processes is then compensated by the setting up of an ageostrophic circulation and related vertical velocities through Equation ( 50).This mechanism does not only apply to SQG flows, but also to any class of quasi-geostrophic ones (even with no surface buoyancy anomalies).However, it is relatively efficient in the SQG case, since at the surface, there is no vertical velocity to oppose the formation of buoyancy gradients, so that there is a permanent frontogenesis that is able to occur.An important point is that there will be large values of vertical velocities within the smallest scales that are affected by this process since they concentrate the buoyancy gradients.
A final remark is that the frontogenesis process can be related to the transfer of potential energy to kinetic energy as discussed in Section 4.4.Using the omega equation and after some manipulation, one can obtain a relation between the spectral flux Π a and frontogenesis: in non-dimensional form and using the integration in the Fourier space [50].This equation shows that Π a is directly connected to the production of buoyancy gradients.Due to the cascade of surface buoyancy to small scales, this term will be positive and explains the transfer of potential energy to kinetic energy at these scales.

Impact on Tracers
As discussed in Section 6, tracer properties differ when comparing standard barotropic and SQG dynamics.At first approximation, phytoplankton can be considered as a passive tracer being stirred by the eddy field once it has sufficiently grown.Phytoplankton spectra at the ocean surface can be measured using satellite images of ocean color.In particular, Gower et al. [124] observed a typical spectrum in k −2 for this tracer and proposed to relate it to the theory of Charney [1].However, such a theory would argue towards a much flatter spectrum in k −1 , and Lesieur and Sadourny [87] proposed an alternative interpretation based on the SQG model.Indeed, the spectrum of a tracer horizontally advected by an SQG flow would be in k −5/3 in the direct cascade of buoyancy variance [17].Horizontal transport is not the sole process to explain phytoplankton variability.Phytoplankton concentrates in the surface oceanic layers, but needs the nutrient uptake from deeper layers to grow.The coupling between horizontal transport (responsible for stirring and for the formation of small horizontal scales) and vertical transport (responsible for the injection of nutrients at the surface) is thus critical in this process.To understand this coupling, one can examine the time evolution and spatial organization of a tracer C at the ocean surface following this equation: where we consider the tracer constant equal to one below the depth z = −H.There are two concomitant effects: the effect of vertical transport (related to w) that injects the tracer in the surface layers and the effect of horizontal transport that horizontally redistributes the tracer concentration.Since the horizontal deformation leads to small-scale formation and is responsible for vertical velocities through Equation ( 50), vertical velocities will occur in filamentary regions between large-scale vortices.Once the tracer is advected from the deeper layers inside these filaments, it will be rapidly spread horizontally over a large distance.Simulations to study this process were examined by Lapeyre and Klein [125], who highlighted the role of eddy-eddy interactions in making this process more efficient.They also showed that half of the tracer can be trapped in these sub-mesoscale filaments.In order to see the role of vertical injection for the phytoplankton dynamics, a biological model was coupled with an SQG model [126].More specifically, the competition between two species of phytoplankton was examined, and Perruche et al. [126] showed that different species of phytoplankton can dominate in different oceanic regions, i.e., either inside vortices or inside filaments, which create ecological niches for plankton growth.These studies are a first step to understand how sub-mesoscales can affect biological dynamics, but their potential role at a more global scale still needs to be quantified.

Surface Quasi-Geostrophy (SQG) and the Atmosphere Dynamics
Oceanic dynamics is quite different from atmosphere processes, in the sense that the ocean is generally thought to be driven both by interior PV dynamics for its interior part and surface buoyancy dynamics for its surface part.On the contrary, the atmosphere is in general interpreted through the interaction of potential anomalies at the tropopause and at the surface [4].Juckes [127] considered the case of a fluid composed of a troposphere and a stratosphere with uniform potential vorticities.Using the SQG approximation, he showed that for an anomaly of potential temperature at the tropopause θ tp , there is a net vertical displacement of the tropopause δz, such that: where N t and N s are respectively the tropospheric and stratospheric Brunt-Väisälä frequencies and θ 0 a reference potential temperature.As illustrated in Figure 5, this relation is revealed to be well satisfied in atmospheric general circulation models [127].In addition, Juckes [127] showed that the streamfunction ψ at the tropopause can be related to θ tp in the Fourier space, such that: More general solutions for a continuous profile of N across the tropopause were discussed by Smith and Bernard [128] and Plougonven and Vanneste [129].Equation (57) shows that the tropopause acts as a material surface because of the discontinuity in stratification, and therefore, the tropopause temperature serves as an active tracer.Moreover, it has much similarity with the oceanic case for which: with ρ the density anomaly.Indeed, the two results are the same if N t tends to infinity.Rivest et al. [130] obtained similar equations when examining the linear solutions of an instability of a shear flow.
The applicability of such solutions was investigated by Tomikawa et al. [131], who confirmed its relevance for explaining the structure of tropopause perturbations close to the sub-polar jet in the Southern Hemisphere.Different applications of SQG for the understanding of the atmospheric dynamics have been done.With a similar approach to Lapeyre and Klein [125], Wirth et al. [132] proposed an interpretation for the roll-up of stratospheric intrusion that is observed on satellite images of water vapor.The destabilization of a potential temperature strip gives rise to the formation of horizontal vortices and, at the same time, to vertical advection of water vapor, which is in agreement with observed water vapor images.Another problem where SQG dynamics has been invoked concerns the horizontal kinetic energy and potential temperature variance spectra near the tropopause that were observed by Nastrom and Gage [133].This study shows a transition from k −3 for large scales to k −5/3 for scales below 500 km.However, the classical quasi-geostrophic theory would indicate a spectrum in k −5/3 at scales larger than the deformation radius (3000 km in the atmosphere) and k −3 for smaller ones.Different theories [134,135] have been proposed to explain the discrepancy with the Charney [1] theory.In this context, Tulloch and Smith [136] offered an interpretation based on a depth-limited version of SQG.Instead of choosing a condition of the semi-infinite domain in z, one can use a domain with a fixed depth H with no buoyancy anomaly at the surface z = 0 and an SQG dynamics at the tropopause z = H.In this case, relation (24) becomes: in dimensional units and with K the wavenumber modulus.For small values of K, potential temperature is proportional to relative vorticity ( θ ∝ −K 2 ψ( k, H)), while for large values of K, we recover the standard SQG model ( θ ∝ −K ψ( k, H)).Tulloch and Smith [136] confirmed that this model could explain the Nastrom and Gage spectrum, and they later proposed a more elaborate version including the role of potential vorticity anomalies [107].However, Asselin et al. [137] showed that it is not possible to rely only on SQG dynamics, as this model leads to further static instabilities near the tropopause if the horizontal speed at the tropopause is large enough.This prevents the application of the SQG theory to scales below 100 km.The more general question of the vertical structure of the troposphere can be examined through an approach of coupling surface-induced and interior-induced dynamics.Morss et al. [138] studied the tropospheric dynamics forced by a baroclinic jet in a quasi-geostrophic model resolving both surface and interior dynamics.They observed the instability of filaments at the tropopause similar to the one for SQG flows [9].They noted a transition between an interior and a surface dynamical regime, which suggests that both interior potential vorticity and tropopause potential temperature control the flow characteristics.Greenslade and Haynes [139] used a similar model to study the transport barrier across a baroclinic jet.Within the tropopause, the jet serves as a barrier to transport, while near the surface, the jet is responsible for strong horizontal mixing, which confirms a vertical transition between the two dynamics.For this situation, the surface buoyancy is the active tracer.This separation between two regimes depending on the flow depth is similar to the oceanic case discussed by Roullet et al. [104].

Conclusions
This review concerned surface quasi-geostrophic flows, a special class of quasi-geostrophic flows, which is based on the conservation of surface buoyancy with uniform potential vorticity within the fluid.The SQG dynamics allows a better understanding of the dynamics of the upper 500 m of the ocean between 10 and 300 km.For the atmosphere, several studies showed that surface quasi-geostrophy could provide some understanding of the tropopause dynamics, but work is still needed to further address its usefulness.Surface buoyancy serves as an active tracer, as is the case for relative vorticity for two-dimensional flows.However, it is associated with a more local dynamics, in the sense that the velocity field associated with a buoyancy anomaly decreases more rapidly as a function of distance compared to a vorticity anomaly for the 2D classical case.This allows small scales to be more energetic as they weakly perceive the effects of other dynamical structures.Concerning SQG turbulence, there exists a direct cascade of surface buoyancy variance toward small scales and an inverse cascade of total energy to large scales.In contrast, for barotropic flows, there is a direct cascade of enstrophy to small scales and an inverse cascade of kinetic energy to large scales.However, while in 2D flows, nonlocal interactions between scales are important in the direct cascade of enstrophy, the contribution of these nonlocal interactions in the direct cascade of buoyancy variance in SQG is not yet clear.An important aspect of SQG flows is that the large-scale component of the flow is relatively inefficient to inhibit filament instability.This may lead to the cascade of instabilities with explosive thinning and roll-up of filaments at smaller and smaller scales.The regularity of SQG solutions in that they do not develop infinite buoyancy gradients in finite time is still a subject of study, in particular because SQG buoyancy gradients bear some formal analogies with the vorticity vector in 3D flows.Still, some results were obtained about the regularity in several cases of dissipation or about the existence of weak solutions.Finally, the SQG equations were used in different (more exotic) contexts, such as the kinetic dynamo [54] or oceanic wind-driven gyres [140].
Several extensions to the SQG model exist.First, one can think of models with finite depth and uniform PV.In this case, two buoyancy equations (at the top and bottom of the fluid) drive the dynamics and the interactions are through the advection of the top (respectively bottom) buoyancy by the velocity induced by the bottom (respectively top) buoyancy, as given by relation (59).This is the essence of the Eady [4] model of baroclinic instability.Such a model was used by Tulloch and Smith [136] to interpret the energy spectra observed at the tropopause.Similar models (but using two layers with homogeneous PV in each one) were developed for understanding of oceanic mixed layer instabilities [110,111].A more general class of QG models consists of coupling the surface buoyancy with interior potential vorticity dynamics [141].The linear instability of a vertical shear flow in the presence of a large-scale potential vorticity gradient (such as the β effect) was addressed by different authors [141][142][143].These studies pointed out that one can understand the instability of such a flow as the interaction of two Rossby waves, one at the surface and one at some critical level.A relevant parameter to determine the typical height of the critical level is the Charney scale: where ∂ z U is the large-scale mean vertical shear and ∂ y PV is the large-scale mean meridional PV gradient [144].The corresponding horizontal scale is λ = Nh/ f 0 .For z h, the dynamics will be surface buoyancy-driven, while for z h, it will be PV-driven.This is consistent with the results of Tulloch and Smith [107], who considered a coupled turbulent model and were able to find this transition in the energy spectra (see also [104,138,139]).Another consequence in the coupling of the surface and the interior dynamics is the barotropization of the flow [145].However, energy transfers between SQG, barotropic and baroclinic modes and between scales are far from obvious.
In particular, one could ask about the possibility of energy isotropization between horizontal and vertical wavenumbers as in the case of Charney [1].Concerning coherent structures, Perrot et al. [146] and Reinaud et al. [147] examined the interactions between surface and interior vortices and filaments using contour dynamics, while Lim and Majda [22] examined the dynamics of point vortices.This latter study proved that clusters of point vortices coupling surface vortices and vortices as a particular depth can exist, such that the cluster cloud propagates as a long-lived and non-dispersive structure.This behavior also occurs for vortex clusters at two different depths [22].In conclusion, these different studies point out that coupling the surface and the interior leads to potentially new phenomena that were not predicted from the classical theory of Charney [1].
As was shown by Williams [148], frontogenesis in a primitive equation gives rise to a rapid contraction of a baroclinic unstable front and to the formation of a discontinuity.This is not the case for the quasi-geostrophic frontogenesis, for which the discontinuity may appear in infinite time, as shown by Williams and Plotkin [119] and as was discussed in this review.Hoskins and Bretherton [122] indicate that, for a baroclinic unstable flow, the formation of the discontinuity occurs after a time τ = N/(ζ ∂ z U) with ζ relative vorticity and ∂ z U vertical shear.Taking typical values for the ocean and the atmosphere leads to τ ≈ 50 days in the ocean and τ ≈ 12 h in the atmosphere.Therefore, the oceanic fronts are in general much weaker than their atmospheric counterparts, and the surface quasi-geostrophic approximation remains a good approximation for the ocean.To include the ageostrophic mechanisms related to atmospheric frontogenesis, the SQG model was extended for flows with finite Rossby numbers: these correspond to semi-geostrophic flows [149], which also assume uniform PV, but which also include ageostrophic effects in the advection of surface buoyancy.In this case, there is a possible formation of singularity in finite time [122].Badin [150] examined the properties of the semi-geostrophy solution in an oceanic context, while Ragone and Badin [151] studied the turbulence in this system and the development of vortex asymmetries.Hakim et al. [117] also examined the same problem using an ageostrophic model different from semigeostrophy (the SQG +1 model).These different studies put in evidence the asymmetry towards more intense cyclones and a restratification of the flow near the surface (an effect obtained in primitive equations models by [152]).Comparisons between SQG and primitive equation solutions were examined in different contexts.Juckes [153] compared both models for the evolution of a shear line with uniform buoyancy at the tropopause.He found that the SQG solution gives a good approximation for small displacements of the tropopause and within the limit of small Rossby numbers.Bembenek et al. [154] examined the time-evolution of an elliptical surface-intensified vortex in primitive equations compared to its SQG counterpart.They showed that both models agree for a small Rossby number, and ageostrophy leads to the inhibition of filament instabilities and the radiation of gravity waves, as well as gravitational and inertial instabilities.The time evolution of a vortex dipole in initial SQG equilibrium was considered by Snyder et al. [155], who showed the generation of inertia-gravity waves with an expansion of the anticyclone and a contraction of the cyclone.To conclude on finite Rossby numbers, even if SQG is appealing for understanding the basic aspects of frontogenesis, it lacks the formation of true discontinuity in the buoyancy field and other ageostrophic mechanisms that affect cyclone-anticyclone asymmetries, non-geostrophic instabilities or inertial gravity wave emission.
As a final comment, surface quasi-geostrophy has been the subject of numerous works in the last twenty years because of the simplicity of its equations bearing analogies with both 2D and 3D Euler equations.More importantly, it attracted the community of geophysical fluid dynamics, as it was found that surface buoyancy gradients are an important ingredient to understand both the upper oceanic layers and the tropopause dynamics.New directions go towards the interaction of interior potential vorticity and surface buoyancy dynamics and towards a better understanding of oceanic surface dynamics using more elaborate models or new observations.

Figure 1 .
Figure 1.Cascade of filament instabilities near the time of singularity (infinite curvature).On the bottom left (box with white background), two elliptical vortices are elongating a thin filament.Each panel in clockwise order shows a close-up of the filament.Finer and finer scale structures can be observed with the rolling-up of the filament.Reprinted with permission from [31].Copyright (2014) by the American Physical Society .

Figure 2 .
Figure 2. (a) Surface buoyancy in a freely-decaying simulation at a resolution of 1024 × 1024.Note that the vortices at all scales develop from filament instability.(b) Spectral fluxes Π θ (continuous curve), Π u (dashed curve) and Π a (dash and dotted curve).See the text for the definition.Panel (b) is from [50], with permission of Cambridge University Press.

Figure 3 .
Figure 3. (a) Surface buoyancy θ + βy in a forced Surface Quasi-Geostrophy (SQG) simulation in the presence of a β effect; (b) corresponding relative vorticity.Only a subdomain is shown.The model is forced at wavenumber k f = 40 for a domain size of [0, 2π] × [0, 2π] and a resolution of 512 × 512.

Figure 4 .
Figure 4. Validation of the SQG balance in the Gulf Stream region.(a) Sea surface temperature (shadings) and sea surface height (white contours); (b) corresponding velocity field, deduced from altimetry (blue arrows) and from SQG balance (red arrows).Note the good agreement between the two velocity fields.From [99], with permission of John Wiley and Sons.

Figure 5 .
Figure 5. Validation of the SQG balance at the tropopause.(a) Geopotential height contours of the Potential Vorticity (PV) = −2 surface (500-m intervals, heavy contour 10 km).The shading is the region where potential temperature lies between 310 and 330 K. Simulations taken from a General Circulation Model.(b) Contoured scatterplot of potential temperature anomalies (from zonal mean) versus geopotential height anomalies on the tropopause (PV = −2 PV units).In panel (a), there is a clear correspondence between geopotential height contours and potential temperature shadings.In panel (b), we clearly see the linearity between geopotential height and potential temperature at the tropopause.From [127].c American Meteorological Society.Used with permission.