2.1 General equations

We start from the SG equations in the $f$ -plane in Boussinesq form. For a complete analysis of the properties of this system of equations and a general overview on SG theory, the reader can refer to Cullen (Reference Cullen2006). Introducing the velocity field $(u,v,w)$ , the streamfunction (rescaled pressure) ${\it\phi}$ , the potential temperature anomaly ${\it\theta}$ and background ${\it\theta}_{0}$ , and the gravitational acceleration $g$ , the SG equations in dimensional form are

where

and the geostrophic velocity field $(u_{g},v_{g})$ is defined by the streamfunction as

The model conserves the potential temperature ${\it\theta}$ and the SG potential vorticity

where ${\bf\zeta}_{sg}$ is the SG absolute vorticity

with $J_{ab}$ the Jacobian operator with respect to the variables $a$ and $b$ . The definition of SG vorticity differs from the QG case for the presence of the nonlinear terms in (2.5). The model conserves also the energy integral over the entire domain,

where the horizontal kinetic energy involves only the geostrophic velocities and is identical to the QG case.

The basic difference between QG and SG is that the materially conserved quantities in QG are advected by the geostrophic velocity field $(u_{g},v_{g})$ , while in SG they are advected by the full (geostrophic and ageostrophic) horizontal velocity field $(u,v)$ . However, $u$ and $v$ are implicit in the SG equations; therefore (2.1) as it is cannot be written as a system of conservation equations with an inversion equation for the streamfunction. Hoskins & Bretherton (Reference Hoskins and Bretherton1972) and Hoskins (Reference Hoskins1975), inspired by previous works by Eliassen (Reference Eliassen1948) and Fjortoft (Reference Fjortoft1962), have shown that, introducing the geostrophic coordinates

the horizontal advection operator expressed in the new coordinates includes explicitly only the geostrophic velocities,

The Bernoulli function

then acts as a streamfunction in the new coordinate system,

Equations (2.7) and (2.9) define a Legendre or contact transformation, and their mathematical properties have been studied by Blumen (Reference Blumen1981) and Purser (Reference Purser1993, Reference Purser1999). In the geostrophic space $(X,Y,Z)$ the dynamics can be expressed in terms of conservation equations for the potential temperature ${\it\theta}$ and the potential vorticity $Q_{sg}$ , plus the inversion equation for the Bernoulli function,

The SG problem in this form can in principle be approached as the QG problem, with two crucial differences: (1) the inversion equation for the streamfunction is a nonlinear Monge–Ampère type equation instead of a linear elliptic equation; and (2) the model is formulated in geostrophic coordinates, which include implicitly the advection by the geostrophic velocity field.

2.2 Boundary conditions and finite-depth surface semi-geostrophic equations

The finite-depth SSG approximation is obtained by defining the domain and the boundary conditions as for finite-depth SQG, and setting constant potential vorticity in the interior of the domain (Badin Reference Badin2013). The final equations are consistent with, for example, Hoskins & West (Reference Hoskins and West1979) and Snyder et al. (Reference Snyder, Skamarock and Rotunno1991), with a few differences. We have chosen the finite-depth version of SSG rather than the version defined on a semi-infinite domain because in the former case the solution procedure is much easier to treat numerically, as discussed in the following. We consider a domain doubly periodic in the horizontal directions and bounded in the vertical by two horizontal surfaces at $Z=0$ and $Z=H$ . At the boundaries we impose the rigid lid condition $w=0$ . Setting constant stratification and potential vorticity $Q_{sg}$ uniform in the interior of the domain, the time-dependent problem occurs only at the boundaries as horizontal advection of ${\it\theta}$ . Fixing ${\it\theta}=0$ at $Z=H$ (Tulloch & Smith Reference Tulloch and Smith2006), the state of the system is determined by the evolution of the 2D dynamics of ${\it\theta}$ at $Z=0$ , which provides the boundary condition to (2.11) to determine the full 3D structure of the flow. As in SQG, in SSG the surface potential temperature plays the same role taken by vorticity in the Euler equations, as the active tracer advected by the dynamics (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). Note that here we have chosen the active boundary at the bottom of the domain, consistent with the formalism of the $\text{SQG}^{+1}$ model of Hakim et al. (Reference Hakim, Snyder and Muraki2002) (which is defined in a semi-infinite domain, as classic SQG) and the finite-depth SQG model of Tulloch & Smith (Reference Tulloch and Smith2006). This choice is appropriate for atmospheric applications. If one wants to consider an oceanic application, the active boundary has to be taken at the top of the domain, which results in an inversion of the vertical coordinates. From the structure of the equations, it can be immediately seen that our results hold also in this case, simply changing the sign of the potential temperature.

Note that, in order to compare the results of SSG and SQG turbulence, as well as with the results from $\text{SQG}^{+1}$ turbulence that are present in the literature, we have adopted here the formulation of SG in geostrophic and height coordinates. Although this is the most commonly used form of the equations, the natural formulation of SG in terms of potential vorticity advection and inversion is, instead, in geostrophic and isentropic coordinates. The proofs that SG is well posed construct a mapping from geostrophic and isentropic coordinates to physical space, but cannot require that the boundary of the region where $Q_{sg}$ is constant in geostrophic and isentropic space maps to the physical boundary (Cullen Reference Cullen2006). As a result, the boundary conditions (2.13) could result in an overdetermined problem and there is the possibility that the results are affected by a lack of well-posedness in the problem. Future work will have to clarify these issues.

In order to formulate the model in non-dimensional form, we define as typical vertical length scale the depth of the domain $H$ . The typical horizontal length scale in geostrophic space is taken as the Rossby deformation radius $L=L_{R}=NH/f$ , where $N$ is the Brunt–Väisälä frequency, here taken as constant. Note that, introducing the length scale in geostrophic space, $L$ represents the typical distance between lines of absolute momentum (Craig Reference Craig1993). Denoting the horizontal velocity scale as $U$ , we introduce the geostrophic space Rossby number ${\it\epsilon}=U/Lf\ll 1$ , and we consider a time scale $T=1/{\it\epsilon}f$ larger than the inertial time scale. Setting $Q_{sg}=N^{2}>0$ and with a suitable rescaling of streamfunction and potential temperature (Hoskins & West Reference Hoskins and West1979), the inversion equation in non-dimensional form can be written as

with boundary conditions

Note that the rescaling transforms (2.13) into a homogeneous equation. The time evolution at $Z=0$ is given by

with geostrophic velocities

The non-dimensional geostrophic coordinates and the Bernoulli function are connected to the non-dimensional physical coordinates and to the streamfunction by Equations (2.12)–(2.14) are the same as in Hoskins & West (Reference Hoskins and West1979), with two differences: (1) we keep the full nonlinear form of equation (2.12) instead of neglecting the nonlinear terms; and (2) we take homogenous boundary conditions at the bottom as in the finite-depth SQG model (Tulloch & Smith Reference Tulloch and Smith2006).

Finite-depth SQG features a transition scale corresponding to a critical wavenumber $k_{t}=f/NH$ (Tulloch & Smith Reference Tulloch and Smith2006), such that the kinetic energy spectrum follows a $k^{-3}$ law for $k\ll k_{t}$ and a $k^{-5/3}$ law for $k\gg k_{t}$ . When the length scale $H$ is taken as the whole depth of the domain (as in our case), the transition scale corresponds to the Rossby radius of deformation, that is, the horizontal unit length. In this case the finite-depth SQG spectrum is expected to follow essentially a $k^{-5/3}$ law, since the $k^{-3}$ regime is not represented in the power-law regime of the spectrum (the power-law scaling emerges for scales smaller than the length scale of the vortices dominating the dynamics, that is, the Rossby radius of deformation). In the current work we set the system in these conditions, as we are interested in the dynamics emerging in the frontal regime. It must be noted that, when dealing with finite-depth SSG, we do not know in general what will be the effect on the transition scale due to (1) the presence of the nonlinear term and (2) the coordinate transformation. We may explore this point in future studies.

The invertibility of the coordinate transformation requires the Jacobian of (2.16) to be positive. Following Hoskins (Reference Hoskins1975), the Jacobian of (2.16) is equal to the vertical component of the absolute SG vorticity (2.5). In non-dimensional geostrophic coordinates one obtains that

When the invertibility condition breaks down, the model produces singular solutions. From a physical point of view, the invertibility condition limits the values of the relative SG vorticity to be smaller than $f$ , thus filtering out inertial instabilities.

In general, if $Q_{sg}$ changes sign within the domain, (2.11) is a mixed type non-homogeneous partial differential equation (Tricomi Reference Tricomi1923) of Monge–Ampère type (elliptic for $Q_{sg}>0$ , hyperbolic for $Q_{sg}<0$ ). For example, in the ocean the potential vorticity can assume positive and negative values, e.g. due to the action of down-front winds which act to destroy potential vorticity (Thomas Reference Thomas2005; D’Asaro et al. Reference D’Asaro, Lee, Rainville, Harcourt and Thomas2011). However, for $Q_{sg}<0$ the assumptions behind the SG equations break down; therefore, given the qualitative, theoretical nature of our study, we limit our analysis to the case $Q_{sg}>0$ .

The SSG inversion equation (2.12) with $Q_{sg}>0$ is elliptic and homogeneous, and substantially easier to treat than the more general form (2.11). However, solving (2.12) is still much less straightforward than solving the Laplace equation arising in SQG (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995; Tulloch & Smith Reference Tulloch and Smith2006). Although other methods of solution of the SG equations have been developed in the past (see Cullen Reference Cullen2006, and references therein), in the vast majority of studies present in the literature numerical solutions of the SG equations have been obtained by approximating equations (2.12)–(2.14). With one notable exception (Snyder et al. Reference Snyder, Skamarock and Rotunno1991), when treating the 3D problem, the nonlinear terms of (2.11) and (2.12) have always been neglected, advocating the fact that they were small compared to the linear terms (Hoskins Reference Hoskins1975, Reference Hoskins1976; Hoskins & West Reference Hoskins and West1979). Alternatively, the dynamics has been limited to a 2D vertical plane, thus automatically eliminating the nonlinearity (Hoskins & Bretherton Reference Hoskins and Bretherton1972; Badin Reference Badin2013). In both cases, the difference between QG and SG dynamics was therefore limited to the use of geostrophic coordinates in the latter; see e.g. the early comment in Hoskins (Reference Hoskins1975) about SG theory providing merely a distortion of the QG solutions, its utility being thus essentially limited to justify the use of the formally and practically much simpler QG theory outside its strict range of applicability. However, Snyder et al. (Reference Snyder, Skamarock and Rotunno1991) showed that including the nonlinear terms of (2.12) does have a substantial effect on the development of cyclones. It is therefore arguable that these terms could play a non-negligible role in determining the properties of SSG turbulence.

4.1 Model description

The numerical integration of the SSG model is performed in geostrophic space and it involves two steps. First, starting from an initial condition in geostrophic coordinates $(X,Y,Z)$ satisfying (2.12), the potential temperature ${\it\theta}$ at $Z=0$ is advected with velocities given by the streamfunction ${\it\Phi}$ at $Z=0$ . Then, the streamfunction ${\it\Phi}$ is computed in the whole domain solving (2.12) with the new potential temperature field as boundary condition at $Z=0$ . The solution of (2.12) is given by (3.11) at first order in ${\it\epsilon}$ . This gives a new streamfunction field at $Z=0$ that is used to advect the potential temperature, and so on. The potential temperature in the whole domain can be reconstructed at each time step by differentiating the vertical profile of the streamfunction, but since the values of ${\it\theta}$ in the interior are not involved in the dynamics, this is done only in post-processing. All the fields are then transformed in physical space in post-processing.

The advection of potential temperature at $Z=0$ is performed with the semi-spectral method employed by Constantin et al. (Reference Constantin, Lai, Sharma, Tseng and Wu2012) to study the formation of singular solutions in the SQG equations. Taking advantage of the periodic boundary conditions in the horizontal directions, the potential temperature field at $Z=0$ and time $t$ is approximated as

where $\boldsymbol{X}=(X,Y)$ at $Z=0$ , $\boldsymbol{k}=(k_{X},k_{Y})$ and $N$ is the number of grid points in the horizontal directions. Time evolution is then performed by integrating with a fourth-order Runge–Kutta scheme the prognostic equation for the Fourier transform of the potential temperature $\hat{{\it\theta}}(\boldsymbol{k},t)$ at $Z=0$ . The Jacobian is computed in geostrophic space, after a fast Fourier transform, using the Arakawa discretization (Arakawa Reference Arakawa1966), which guarantees conservation of the energy and enstrophy invariants. In order to remove numerical instabilities, it is common practice to introduce a dissipation operator, typically in the form of a hyperdiffusion. Following Constantin et al. (Reference Constantin, Lai, Sharma, Tseng and Wu2012), we employ instead an exponential filter introduced by Hou & Li (Reference Hou and Li2007), multiplying each Fourier mode $\hat{{\it\theta}}(\boldsymbol{k},t)$ by ${\it\rho}(2k_{X}/N){\it\rho}(2k_{Y}/N)$ , where ${\it\rho}(x)=\exp [-{\it\alpha}x^{m}]$ . We take ${\it\alpha}=36$ and $m=19$ , as in Constantin et al. (Reference Constantin, Lai, Sharma, Tseng and Wu2012). This high-order exponential spectral filter has been developed by Hou & Li (Reference Hou and Li2007) specifically for the discretization of systems of equations developing singular or nearly singular solutions (as the semi-geostrophic equations). With this choice of the parameters, the application of the filter keeps the two-thirds lower-wavenumber modes unchanged and suppresses the one-third higher-wavenumber modes. Note that this filter is a less sharp version of the classical Orszag $2/3$ method, where the coefficients of the highest one third of the wavenumbers are simply put equal to zero. Hou & Li (Reference Hou and Li2007) showed that the use of this high-order exponential filter captures up to $15\,\%$ more effective Fourier modes than the $2/3$ Orszag method, producing more accurate approximations of the solutions. We have performed experiments with different choices of the parameters of the filter, in particular using ${\it\alpha}=512$ and $m=24$ , such that the filter was sharper and thus more effective in removing small-scale features, without observing differences in the results.

The streamfunction is computed by solving the nonlinear Monge–Ampère equation with the method described in the previous section. We limit our investigation to small Rossby numbers, so that we consider only the first term of the expansion resulting from the iterative procedure. The solution of the problem is computed at each time step by numerically computing (3.11). While in SQG the streamfunction can be computed from the potential temperature at the top boundary through a simple inversion in Fourier space, effectively reducing the problem to two dimensions, in SSG the presence of the nonlinear term in the Monge–Ampère equation requires taking explicitly into account the full 3D structure of the streamfunction and including a discretization of the vertical coordinate. However, like in SQG, also in SSG the vertical integration is made from a diagnostic equation and does not need time integration at each vertical layer. We have performed simulations on a $512\times 512$ horizontal square grid and 20 vertical levels. As one can see from (3.11), high-wavenumber modes decay faster with depth, so that high vertical resolution is needed only close to the surface. The vertical levels are exponentially spaced, with the layer depths varying from ${\rm\Delta}Z=0.004$ to ${\rm\Delta}Z=0.18$ from the top to the bottom.

The initial condition is defined as a random field of surface ${\it\theta}$ in geostrophic space, as common practice in the literature on the SG approximation. The random field of surface ${\it\theta}$ is defined as in Hakim et al. (Reference Hakim, Snyder and Muraki2002), such that the corresponding zeroth-order surface streamfunction follows

with $m=25$ and $k_{0}=14$ . A random phase is given to each mode, in order to have a random initial condition with a prescribed kinetic energy spectrum. After several tests aimed at avoiding cases with too many points at which the invertibility condition was violated, we have selected the cases with surface kinetic energy normalized at $KE_{sg}=5$ . Note that here as well as in the following we consider the kinetic energy associated with the geostrophic velocity field $(u_{g},v_{g})$ , as this is the quantity that enters in the conserved energy integral (2.6). The time step is taken as $\text{d}t=0.005$ , and the simulations are performed up to $T=100$ , corresponding to approximately 300 eddy turnover times, computed according to Kerr (Reference Kerr1990). We performed a sensitivity analysis on the Rossby number, performing simulations for ${\it\epsilon}=(0.02,0.05,0.1,0.15,0.2)$ . The range of values of ${\it\epsilon}$ is in agreement with the range chosen by Snyder et al. (Reference Snyder, Skamarock and Rotunno1991), who consider a maximum value ${\it\epsilon}=0.3$ . We compare the results also with runs of the standard finite-depth SQG.

Note that our experiments do not include forcing and large-scale dissipation. Among other works that have performed freely decaying simulations with SQG and SQG-like models, Hakim et al. (Reference Hakim, Snyder and Muraki2002) limited the analysis at long times, while Capet et al. (Reference Capet, Klein, Hua, Lapeyre and McWilliams2008a ) limited the analysis at early times during which nonlinear interactions are rather strong. Both note that for long times the system shows kinetic energy spectra much steeper than those predicted by the theory. For shorter times the spectra of Capet et al. (Reference Capet, Klein, Hua, Lapeyre and McWilliams2008a ) are closer to the theoretical prediction. However, Capet et al. (Reference Capet, Klein, Hua, Lapeyre and McWilliams2008a ) study the transfer of surface kinetic energy in SQG flows, and their analysis focuses on a stage of the evolution that is far from approximating a stationary state of the SQG equations (they have almost no vortices formed yet). We have therefore chosen the approach of Hakim et al. (Reference Hakim, Snyder and Muraki2002). We have restricted our analysis to the last 10 time units (corresponding to approximately 30 eddy turnover times), where we have tested that the statistics of the system does not change with time in a significant way. In this regime, energy is reduced to approximately 30 % of its value at time zero, but further decays by less then 5 % within the time window where the analysis is performed. Still, our freely decaying flow does not reach a real stationary state, and we expect to observe SQG kinetic energy spectra deviating from the theoretical prediction, as in Hakim et al. (Reference Hakim, Snyder and Muraki2002).

Transformation from geostrophic to Cartesian coordinates is performed in post-processing, checking if the integration produces singularities (Schär & Davies Reference Schär and Davies1990). We have checked that in the full turbulent state the invertibility condition is violated in the worst cases in a few per cent of the total points of the domain only, in the upper layers of the model and for the largest values of the Rossby number. The local formation of singularities is unavoidable in numerical simulations of the SG equations, and the problem occurs also when employing different methods of solution that do not involve the coordinate transformation (see Cullen Reference Cullen2006, and references therein). However, in the numerical model these singularities are systematically eliminated by the spectral filter, so that they are not problematic from the stability point of view. We have performed a large number of tests varying the total kinetic energy of the initial condition (see above), and selected those simulations for which the points at which the invertibility condition was not satisfied in the fully turbulent phase of the flow remained few and isolated. Note that the problem of violation of the invertibility condition at the initial condition could be avoided by formulating the model in isentropic coordinates. In this way, defining $Z=g{\it\theta}/{\it\theta}_{0}$ , SSG would be determined by imposing the condition that $Q_{sg}$ is constant on a region periodic in $(X,Y)$ with $Z_{0}(X,Y,t)<Z<Z_{1}$ , where $Z_{1}$ is constant. The evolution equation would then apply as an equation for $Z_{0}$ rather than for ${\it\theta}$ (although with different boundary conditions). In this case one could use any random field $Z_{0}$ as initial condition without violating the invertibility condition, thus avoiding the initial dissipation. However, one should note that local violations of the invertibility condition and dissipation would anyway emerge during the evolution of the flow due to numerical effects. Further, this would not allow for comparison of the results of this study with classical studies of SQG turbulence.

Figure 1. Snapshots at $T=100$ of ${\it\theta}$ at $z=0$ for SQG (a), SQG under coordinate transformation with ${\it\epsilon}=0.2$  (b), SSG for ${\it\epsilon}=0.2$ in geostrophic coordinates (c) and SSG for ${\it\epsilon}=0.2$ in physical coordinates (d).

4.2 Surface statistics

Figure 1 shows snapshots of surface ${\it\theta}$ at $T=100$ for (a) SQG, (b) SQG under coordinate transformation with ${\it\epsilon}=0.2$ (corresponding to an SSG simulation as traditionally proposed in the literature where only the coordinate transformation is considered), (c) SSG for ${\it\epsilon}=0.2$ in geostrophic coordinates (where only the effect of the inclusion of the nonlinearity of the inversion equation is visible) and (d) SSG for ${\it\epsilon}=0.2$ in physical coordinates (the full SSG case). While the range of variability is one order of magnitude larger, the scale of ${\it\theta}$ is limited to the range $[-0.2,0.2]$ , in order to highlight structures with lower intensity rather than the vortices. Note that ${\it\theta}$ is the active conserved scalar in SQG and SSG, so that it takes the role normally taken by vorticity. At $Z=0$ , SQG (a) produces mostly localized coherent structures with filaments formed by secondary instabilities in between. In comparison, SSG turbulence (d) is a mixture of local structures, represented by the coherent vortices, and nonlocal structures, represented by fronts and filaments dominating the dynamics. The filaments are not produced by secondary instabilities but rather by the organization in features with skewed potential temperature, induced by the joint action of the nonlinear term and the transformation of coordinates. Cyclones (region with positive ${\it\theta}$ anomaly) are smaller than anticyclones and isolated. The appearance of an asymmetry in the distribution and size of cyclonic and anticyclonic regions is an expected result of the coordinate transformation. As can be seen by comparing the SQG snapshot (a) and the SQG snapshot under coordinate transformation (b), the effect of the coordinate transformation is that positive relative vorticity (potential temperature) is increased in magnitude and the areas where it occurs are compressed, while negative relative vorticity (potential temperature) is decreased in magnitude and the areas where it occurs are expanded (Hoskins Reference Hoskins1975). On the other hand, the different vortex dynamics induced by the nonlinear term is clear on comparing the SQG snapshot (a) and the SSG snapshot in geostrophic coordinates (c), which show how the nonlinear terms act to enhance the negative values of ${\it\theta}$ , corresponding to the strain-dominated regions of the flow.

Figure 2. (a) PDFs of surface ${\it\theta}$ for SQG (blue line) and SSG for different values of Rossby number in physical coordinates (red lines). In the inset we show as a function of the Rossby number the mean and median of the distributions. The SQG case corresponds to the origin. (b) Surface kinetic energy spectra for SQG (blue line) and SSG for different values of Rossby number in physical coordinates (red lines).

Figure 2(a) shows the probability density functions (PDFs) of surface ${\it\theta}$ for SQG (blue line) and SSG for different values of Rossby number in physical coordinates (red lines). While SQG produces a non-Gaussian but zero-centred PDF, SSG produces instead skewed PDFs that are centred around non-zero negative values. Since the isolated coherent structures control the tails of the distribution, it is the filaments that are responsible for this asymmetry in the bulk of the distribution. The increase in magnitude of the mean and median of the PDFs with ${\it\epsilon}$ shows an almost linear trend (in the inset), as a signature of a continuous shift towards negative values of the potential temperature anomaly as ${\it\epsilon}$ increases. This effect has been observed also in Hakim et al. (Reference Hakim, Snyder and Muraki2002), and it is linked to a net cooling of the surface. The net cooling is due to the increase in surface kinetic energy that accompanies the forward cascade of buoyancy variance, which implies a decrease in available potential energy, which lowers the centre of gravity of the fluid. In an oceanic case, where the temperature anomaly has the opposite sign, this effect is associated with the restratification effect associated with filamento- and frontogenesis.

Figure 2(b) shows the surface horizontal kinetic energy spectra for SQG (blue line) and SSG for different values of Rossby number and physical coordinates (red lines). In general the spectra at level $Z$ have been computed taking the radial average of the kinetic energy spectral density field,

Andrews & Hoskins (Reference Andrews and Hoskins1978) predicted a semi-geostrophic $k^{-8/3}$ power law for both the kinetic and available potential energy of a one-dimensional (1D) front reaching an infinite value for vorticity. Following the behaviour of the front past the singularity, Boyd (Reference Boyd1992) corrected the value to $k^{-2}$ . Both these values refer, however, to the case of an isolated 1D front.

Finite-depth SQG spectra are supposed, in the chosen set-up, to follow a $k^{-5/3}$ law, as the transition between $k^{-3}$ and $k^{-5/3}$ occurs at $k\approx 6$ , so that the regime $k^{-3}$ is not represented. However, our SQG spectra are steeper than $k^{-5/3}$ , due to the fact that our free decaying flows do not reach a real steady state, as the energy dissipation at small scales by the spectral filter is not compensated by injection of energy by an external forcing. Without a forcing, the flow does not stabilize in a scale-invariant stationary state, but tends to develop coherent structures that prevent the spectra from attaining the theoretical slope. The same effect is observed and discussed in Hakim et al. (Reference Hakim, Snyder and Muraki2002) and Capet et al. (Reference Capet, Klein, Hua, Lapeyre and McWilliams2008a ), who have performed similar freely decaying simulations of the SQG and $\text{SQG}^{+1}$ models. SSG spectra are less steep and tend to $k^{-5/3}$ as ${\it\epsilon}$ increases, in agreement with the emerging role of filaments, which implies that more energy is stored at smaller scales for larger values of ${\it\epsilon}$ . SSG spectra show also the presence of ‘bumps’, which might be an indication of the enhanced energy present in coherent structures associated with the fact that in SSG regions with positive and negative potential temperatures are not symmetric.

4.3 Comparison with the SQG $^{+1}$ model

The SSG model presents several similarities with the $\text{SQG}^{+1}$ model of Hakim et al. (Reference Hakim, Snyder and Muraki2002), although the two models are rather different mathematically. It is therefore of interest to compare the results of the two models. The similarities come from the fact that both models can be seen as extensions of the SQG model, taking into account the effect of ageostrophic advection at first order in the Rossby number. The way in which the ageostrophic advection is introduced is, however, very different in the two models. In particular, the SSG model involves a coordinate transformation that introduces a deformation dependent on the vorticity field, which has the effect of favouring the formation of frontal structures. Further, while both SSG and $\text{SQG}^{+1}$ are second-order accurate approximations to the Euler equations, $\text{SQG}^{+1}$ requires ${\it\epsilon}\sim Fr$ , where $Fr$ is the Froude number, while SSG requires ${\it\epsilon}\sim Fr^{2}$ . As a consequence of this, $\text{SQG}^{+1}$ employs a linearization of the static stability profile about a reference value, while SSG does not. It would thus be expected for SSG to perform better than $\text{SQG}^{+1}$ at describing the effect of static stability variations, but to be deficient in describing vortex instabilities.

The PDFs of ${\it\theta}$ of SSG and $\text{SQG}^{+1}$ at the surface present the same kind of asymmetry. In both cases, the PDFs are peaked at moderately negative values, with both the median and the mean taking negative values. This net cooling is due to the restratification effect. In SSG it is possible to disentangle the effects of the nonlinear term and of the coordinate transformation, showing that the net cooling is induced solely by the coordinate transformation, and thus the advection by the ageostrophic flow, while the shift of the peak of the distribution is due to the nonlinear term. It is not therefore possible to connect $\text{SQG}^{+1}$ to only one of these two aspects of SSG. On the other hand, SSG and $\text{SQG}^{+1}$ strongly differ in the kinetic energy spectra. In SSG the coordinate transformation leads to more energy being stored at small scales, with a consequent flattening of the SSG spectra with respect to the SQG case. On the contrary, the $\text{SQG}^{+1}$ has exactly the same slope as the SQG spectra. Both models similarly affect the population number and morphology of the cyclones and/or anticyclones. Hakim et al. (Reference Hakim, Snyder and Muraki2002) performed a detailed analysis of the vortex statistics and structure, making use of algorithms for vortex census and estimation of the vortex radius. It would therefore be of great interest to perform a similar analysis also on the SSG dynamics in a future work.

4.4 Interior statistics

The analysis of the interior statistics is of particular interest, as one of the reasons for the appeal of SQG-like models is the possibility to infer properties of the ocean at depth from observations of surface fields (LaCasce & Mahadevan Reference LaCasce and Mahadevan2006; Lapeyre & Klein Reference Lapeyre and Klein2006; Wang et al. Reference Wang, Flierl, LaCasce, McClean and Mahadevan2013; Liu et al. Reference Liu, Peng, Wang and Huang2014). Although rigorously SSG cannot be applied to represent quantitatively processes occurring at scales smaller than the Rossby radius of deformation, the ability to reproduce the formation of frontal structures and filaments goes in the right direction to represent submesoscale dynamics that are not captured by the SQG approximation. Indeed, SSG in a 2D vertical plane has proved to be successful in correcting the underestimation of buoyancy anomalies at depth that characterize the SQG interior profiles (Badin Reference Badin2013), even if the effects observed by Badin (Reference Badin2013) are due to ${\it\epsilon}\rightarrow 1$ , while this study is restricted to small values of ${\it\epsilon}$ . Moreover, the SG approximation allows for $O(1)$ variations of static stability, which are physically important in determining the properties in the interior of the ocean from the surface mixed layer.

Figure 4. Vertical profiles of moments of ${\it\theta}$ distributions for SQG (blue lines) and SSG for different values of the Rossby number in geostrophic (black lines) and physical (red lines) coordinates: median (a), mean (b), standard deviation (c), skewness (d).

Figure 5. (a,c) PDFs of ${\it\theta}$ for SQG (blue lines) and SSG in physical coordinates (red lines) for different values of Rossby number at $Z=0.1$  (a) and $Z=0.8$  (c). (b,d) Energy spectra for SQG (blue lines) and SSG in physical coordinates (red lines) for different values of Rossby number at $Z=0.1$  (b) and $Z=0.8$  (d). The solid black line shows the expected behaviour of SQG spectra at the corresponding depth assuming a $-5/3$ slope of the spectrum at the surface.

Figure 6. (a) Vertical profile of fraction of horizontal kinetic energy due to the first-order correction to the SQG-like solution, in geostrophic (black lines) and physical (red lines) coordinates. The blue line corresponds to the SQG case. (b) Fraction of horizontal kinetic energy due to the first-order correction to the SQG-like solution as a function of ${\it\epsilon}$ for different values of $Z$ , in geostrophic (black lines) and physical (red lines) coordinates. Note that at $Z=0.8$ the black and red lines are identical, as the coordinate transformation has almost no effect.

Figure 7. Two-dimensional PDFs of ${\it\theta}$ and ${\it\phi}$ for SQG (a,c) and SSG with ${\it\epsilon}=0.2$ (b,d), at surface (a,b) and at $Z=0.8$ (c,d) in physical coordinates. The contours are in logarithmic scale.

Figure 8. (a,b) Two-dimensional PDFs of ${\it\theta}$ and ${\it\phi}_{0}$ for SSG with ${\it\epsilon}=0.2$ , at surface (a) and at $Z=0.8$  (b) in physical coordinates. (c,d) Two-dimensional PDFs of ${\it\theta}$ and ${\it\phi}_{1}$ for SSG with ${\it\epsilon}=0.2$ , at surface (c) and at $Z=0.8$  (d) in physical coordinates. (e,f) Two-dimensional PDFs of ${\it\phi}_{0}$ and ${\it\phi}_{1}$ for SSG with ${\it\epsilon}=0.2$ , at surface (e) and at $Z=0.8$  (f) in physical coordinates. The contours are in logarithmic scale.

Figure 4 shows vertical profiles of the median (a), mean (b), standard deviation (c) and skewness (d) of ${\it\theta}$ . The black lines refer to geostrophic coordinates, while the red lines refer to physical coordinates. The results show that the asymmetry of the PDFs strongly weakens at depth. The effect of the coordinate transformation, visible in the differences between the red and black lines, disappears at depth earlier than the effect of the nonlinearity of the Monge–Ampère equation, visible in the difference between the black and blue lines. From (3.11), small scales (i.e. large $k$ ) decay faster with depth at leading order. Therefore, fronts and filaments tend to disappear at depth, restoring the symmetry of the distributions. This is in agreement with the physical interpretation on the role of enhanced ageostrophic, small-scale frontal structures in determining the differences between SSG and SQG.

Figure 5 shows the PDFs of ${\it\theta}$ (a,c) and the kinetic energy spectra (b,d) at $Z=0.1$ (a,b) and $Z=0.8$ (c,d) for SQG (blue lines) and SSG in physical coordinates for different values of the Rossby number (red lines). In the interior of the domain, but still close to the surface ( $Z=0.1$ ), SSG produces PDFs of ${\it\theta}$ that are highly skewed and again peaked at negative values, while at depth ( $Z=0.8$ ) the PDFs appear with zero mean and are nearly symmetric even for high Rossby number.

From (3.11) the kinetic energy spectra in the interior $\mathscr{K}(k,Z)$ are linked at zeroth order to the kinetic energy spectrum at the surface $\mathscr{K}(k,0)$ by a multiplicative factor rapidly decaying for large $Z$ and large $k$ (Callies & Ferrari Reference Callies and Ferrari2013),

Superimposed on the spectra of figure 5 is the SQG-like zeroth-order expected behaviour following (4.5), and assuming a $-5/3$ spectral slope at the surface. The results for SSG at $Z=0.1$ show higher energies at smaller scales as the Rossby number increases, with a slope approaching the $-5/3$ spectral slope at the surface, corrected with depth. At $Z=0.8$ , the SQG and SSG spectra converge for all values of the Rossby number. Figure 5 thus further confirms that the deviation of SSG from SQG disappears with depth.

In classic SQG, kinetic energy and density variance spectra $\mathscr{K}(k,Z)$ and $\mathscr{T}(k,Z)$ are identical (Blumen Reference Blumen1978; Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). In finite-depth SQG (Tulloch & Smith Reference Tulloch and Smith2006), $\mathscr{K}(k,Z)$ and $\mathscr{T}(k,Z)$ are equal at small scales, with a transition occurring at scales larger than a critical wavenumber from a classical SQG-like to a QG $/$ 2D-like behaviour. In our case, as discussed above, only the $-5/3$ regime is properly represented. In finite-depth SSG, things differ due to the presence of the nonlinear term and of the coordinate transformation. From (3.11) we have

Additional contributions to the SQG-like energy spectrum $\mathscr{K}_{0}$ appear due to the presence of the order- ${\it\epsilon}$ part of the SSG solution in (3.11). We have verified that the spectra of the different components of (4.6) all follow the same slope in physical coordinates (not shown), so that the change of the slope with ${\it\epsilon}$ is not due to a flatter slope of the part of the spectrum due to the deviations from the SQG-like solution. Figure 6(a) shows the vertical profile of the fraction of kinetic energy connected to the additional terms ${\it\epsilon}\mathscr{K}_{01}+{\it\epsilon}^{2}\mathscr{K}_{1}$ for different values of the Rossby number in geostrophic (black lines) and physical (red lines) coordinates. Again, the relative importance of the SSG first-order correction term weakens at depth, and the contribution of the coordinate transformation is negligible for $Z>0.2$ . The maximum of the effect, however, is reached not at the surface but at an intermediate depth between $Z=0$ and $Z=0.1$ . Figure 6(b) shows how the fraction of kinetic energy connected to the additional terms ${\it\epsilon}\mathscr{K}_{01}+{\it\epsilon}^{2}\mathscr{K}_{1}$ increases with ${\it\epsilon}$ for $Z=0.1$ (solid lines) and $Z=0.8$ (dashed lines), in geostrophic (black) and physical (red) coordinates. At $Z=0.1$ the increase with ${\it\epsilon}$ is clearly nonlinear; however, the values always remain of the order of ${\it\epsilon}$ , thus confirming the robustness of our partition of the solution for small ${\it\epsilon}$ into a zeroth-order part and a first-order correction term. At $Z=0.8$ the increase is instead much slower, due to the minor importance of ageostrophic processes at depth, and the coordinate transformation does not introduce any contribution.

This has some interesting implications regarding the computation of the flow in the interior of the oceanic mixed layer from knowledge of the surface data, one of the applications of SQG. Let us suppose that we observe a current and buoyancy field at the surface. We can deduce an estimate of the currents and buoyancy in the interior of the domain by using the SQG inversion. This is known to lead to underestimating the reconstructed fields. If we use instead the SSG inversion ((3.11) at first order), the estimate of the streamfunction consists of the zeroth-order term, which is the same as if we had used the SQG inversion, plus an additional contribution. Figure 6 shows that this additional contribution leads to a positive correction to the kinetic energy of the estimated current field corresponding to the reconstructed streamfunction. The additional kinetic energy predicted by the SSG inversion is up to 20–30 % in the upper layers and decays with depth; therefore it is connected with filaments and other small-scale structures not captured by SQG. This could qualitatively explain why the SQG inversion underestimates the fields at depth when applied to observations at the surface.

4.5 The ${\it\phi}{-}{\it\theta}$ relationship

It is interesting to see the relation between the potential temperature and the different components of the streamfunction. Figure 7 shows the 2D PDFs of potential temperature and streamfunction for SQG (a,c) and SSG with ${\it\epsilon}=0.2$ in physical coordinates (b,d), at the surface (a,b) and at $Z=0.8$ (c,d). The contours are in logarithmic scale. The asymmetry at the surface and the different properties highlighted above are easily visible also in these plots. The PDFs of SQG are symmetric in both ${\it\theta}$ and ${\it\phi}$ , while in SSG the cyclonic tails, characterized by large positive values of ${\it\theta}$ associated with large negative values of ${\it\phi}$ , are longer than the anticyclonic tails, characterized by large negative values of ${\it\theta}$ associated with large positive values of ${\it\phi}$ . Still, negative values of ${\it\theta}$ dominate in the mean, due to the strong asymmetries of the distributions. An additional feature visible in figure 7 is the clear signature of coherent structures (vortices) in the PDFs at the surface, showing up as ‘fingers’ in the tails of the distributions, as coherent structures are characterized by a functional relation between streamfunction and advected scalar ${\it\theta}={\it\theta}({\it\phi})$ . While in SQG cyclonic and anticyclonic coherent structures have the same properties, in SSG a clear asymmetry emerges, due to the different deformations of strong cyclonic and anticyclonic areas induced by the coordinate transformation that changes the morphology of the corresponding coherent structures.

Figure 8 shows for the case ${\it\epsilon}=0.2$ the relation between ${\it\theta}$ and the two components of the SSG solution for the streamfunction ${\it\phi}_{0}$ (a,b) and ${\it\phi}_{1}$ (c,d). At the surface (a,c), while ${\it\phi}_{0}$ behaves basically like ${\it\phi}$ , ${\it\phi}_{1}$ has a very clear functional form. The values of ${\it\phi}_{1}$ are almost always negative, and larger in magnitude for positive ${\it\theta}$ . To better understand the role of ${\it\phi}_{1}$ in correcting the SQG-like solution, figure 8(e,f) shows the 2D PDFs of ${\it\phi}_{0}$ and ${\it\phi}_{1}$ at the surface (e) and at $Z=0.8$  (f). We can see that at the surface the effect of ${\it\phi}_{1}$ is to strongly enhance regions of negative ${\it\phi}_{0}$ and weakly damp regions of positive ${\it\phi}_{0}$ . It is then clear how this additional term systematically modifies the streamfunction, and therefore the velocity field, resulting in the observed asymmetric distributions of ${\it\theta}$ and ${\it\phi}$ . Consistently with what is seen previously, at $Z=0.8$ (b,d,f) the correction term ${\it\phi}_{1}$ is totally negligible with respect to ${\it\phi}_{0}$ , as it is associated with small-scale structures that do not penetrate at depth, and SQG-like symmetry is restored.

4.6 Vertical velocity, ageostrophic divergence and lateral strain rate

Following Hoskins & Draghici (Reference Hoskins and Draghici1977), the diagnostic equation for the vertical velocity in non-dimensional variables can be written as

with boundary conditions $w^{\ast }=0$ at $Z=0$ and $Z=1$ , where $w=Jw^{\ast }$ , and the forcing vector $\boldsymbol{Q}$ in non-dimensional form is

The SG equation for the vertical velocity in geostrophic coordinates is formally identical to the equation one derives in physical coordinates in the QG case, with two differences: the use of geostrophic coordinates; and the fact that the equation is formulated for a vertical velocity $w^{\ast }$ rescaled by $J$ . Note also that in the derivation of (4.7), $J$ is approximated neglecting the nonlinear terms (Hoskins & Draghici Reference Hoskins and Draghici1977; Hoskins, Draghici & Davies Reference Hoskins, Draghici and Davies1978). Equation (4.7) is not therefore entirely compatible with the full form of the SSG equations that we are here investigating. In order to be consistent with the derivation of equation (4.7), we have computed $w$ using the approximated version of $J$ , although there is no substantial qualitative difference between the results that we would have obtained using the full form of $J$ . For further studies of semi-geostrophic vertical velocities, the reader should see, for example, Hoskins et al. (Reference Hoskins, Draghici and Davies1978), Hoskins & West (Reference Hoskins and West1979), Pinot, Tintoré & Wang (Reference Pinot, Tintoré and Wang1996), Pedder & Thorpe (Reference Pedder and Thorpe1999), Thorpe & Pedder (Reference Thorpe and Pedder1999) and Viudez & Dritschel (Reference Viudez and Dritschel2004).

Figure 9 shows snapshots at time $T=100$ of the vertical velocity $w$ at $Z=0.015$ for SQG (a) and SSG with ${\it\epsilon}=0.2$ (b) in physical coordinates. Comparing figure 9 with figure 1, we can see that, while in SQG vertical velocity is larger inside the vortices with a typical quadrupole structure, in SSG the largest values are obtained most often on the edges of the vortices and of the frontal structures, as captured in primitive equation simulations and that can be significant for the transfer of nutrients to the surface euphotic layer (Levy, Klein & Treguier Reference Levy, Klein and Treguier2001; Mahadevan & Tandon Reference Mahadevan and Tandon2006). This is due to the fact that $w$ is given by the product of $w^{\ast }$ and $J$ , and the latter is largest exactly on the edges of vortices and of frontal structures.

Figure 9. Snapshot at time $T=100$ of vertical velocity at $Z=0.015$ for SQG (a) and SSG in physical coordinates (b), with ${\it\epsilon}=0.2$ .

Figure 10. (a) Vertical profile of average absolute value of vertical velocity $w$ for SSG for different values of ${\it\epsilon}$ in physical coordinates. (b) Vertical profile of average absolute value of horizontal divergence ${\it\delta}$ for SSG for different values of ${\it\epsilon}$ in physical coordinates. Note the different vertical scale of the two plots.

Vertical velocity decays in magnitude rather fast with depth. Figure 10(a) shows the vertical profile of the average of the absolute value of $w$ for SSG for different values of the Rossby number in physical coordinates. We can see that for increasing Rossby number the vertical velocity becomes larger close to the surface, and penetrates deeper at depth. The largest values of $w$ in the upper layers are dominated by large values of $J$ (not shown), so that stirring-dominated regions become of great relevance for what concerns the vertical velocity. In comparison, vertical velocities in SQG (not shown) show larger values in the interior, due to the larger features that are present in SQG and due to the fact that the vertical decay is proportional to the wavenumber.

The largest values of vertical velocities found in the interior as ${\it\epsilon}$ increases is in agreement with the analytical results found by Badin (Reference Badin2013). This effect would be even more enhanced if, following Badin (Reference Badin2013), the potential vorticity would be expressed as a function of the Rossby number. Further, as stated in § 4.3, the increase of the vertical velocities with ${\it\epsilon}$ in the interior might be limited by the use of height coordinates, as the solutions are forced to be confined close to the boundary. Further work is required to understand this behaviour of SSG.

Figure 11. (a,b) PDFs of ageostrophic horizontal divergence at $Z=0.015$  (a) and at $Z=0.8$  (b), for SQG (blue lines) and SSG for different values of ${\it\epsilon}$ in physical coordinates (red lines). (c,d) PDFs of lateral strain rate at $Z=0.015$  (c) and at $Z=0.8$  (d), for SQG (blue lines) and SSG for different values of ${\it\epsilon}$ in physical coordinates (red lines).

Once the vertical profile of $w$ has been computed, one can easily reconstruct the horizontal divergence. The flow divergence is an important signal for frontogenetic dynamics, as frontogenesis is associated with a divergent flow of warm and cold water/air rising or sinking at the front through the secondary ageostrophic circulation. As the geostrophic velocity field is divergence-free, the horizontal divergence is purely ageostrophic. From the continuity equation, one can compute the horizontal divergence ${\it\delta}$ in physical coordinates as

where $u_{ag}=u-u_{g}$ and $v_{ag}=v-v_{g}$ are the ageostrophic components of the velocity field. Figure 10(b) shows the vertical profile of the absolute value of ${\it\delta}$ for SSG for different values of the Rossby number in physical coordinates. Values of the horizontal divergence substantially different from zero are found only very close to the surface, so that the vertical scale is limited between $Z=0$ and $Z=0.2$ . The magnitude of ${\it\delta}$ increases for increasing Rossby number, similarly to $w$ . Figure 11(a,b) shows the PDFs of ${\it\delta}$ at $Z=0.015$  (a) and $Z=0.8$  (b) for SSG for different values of the Rossby number in physical coordinates. We can see that close to the surface ( $Z=0.015$ ) the PDFs become flatter for increasing Rossby numbers, with larger values in the far tails of the distributions. At depth ( $Z=0.8$ ) the same behaviour occurs, but, given the extremely small values of ${\it\delta}$ below $Z=0.2$ , it is of relatively low interest.

Submesoscale filaments are also characterized by large values of the lateral strain rate. For simplicity, we compute the lateral strain rate ${\it\alpha}$ on geostrophic velocity field only as (Shcherbina et al. Reference Shcherbina, D’Asaro, Lee, Klymak, Molemaker and McWilliams2013)

Figure 11(c,d) shows the PDFs of ${\it\alpha}$ at $Z=0.015$  (c) and $Z=0.8$  (d) for SQG (blue lines) and SSG for different values of the Rossby number in physical coordinates (red lines). As expected, close to the surface ( $Z=0.015$ ) the lateral strain rate in SSG takes larger values than in SQG, increasing as ${\it\epsilon}$ increases, with PDFs characterized by longer tails and less concentrated around small values of ${\it\alpha}$ . This is again a signature of the enhanced role of filaments in the dynamics. The lateral strain rate at depth ( $Z=0.8$ ) shows instead the opposite behaviour, with smaller values of ${\it\alpha}$ at $Z=0.8$ for SSG for increasing Rossby number. Values of ${\it\alpha}$ at $Z=0.8$ are, however, smaller by one order of magnitude than close to the surface, as small-scale structures characterized by large values of the lateral strain rate disappear at depth. In general, it is close to the surface, where small-scale ageostrophic structures are active, that we expect to see a strong impact on the mixing of tracers.