Numerical investigation of the Danilov inequality for two-layer quasi-geostrophic systems

We consider a Boussinesq fluid consisting of two immiscible layers of equal depth on an f-plane and confined between rigid horizontal plates. We assume that the fluid's motion is nearly geostrophic. This means that the fluid's governing equations are the following QG potential vorticity equations:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}{q}_{j} & = & -J({\psi }_{j},{q}_{j})+{\left(-1\right)}^{\gamma -1}{\nu }_{j}{{\rm{\nabla }}}^{2\gamma }{q}_{j}-\kappa {\delta }_{2j}{{\rm{\nabla }}}^{2}{\psi }_{j}\\ & & +\alpha {k}_{{\rm{d}}}^{2}\left[({\psi }_{j}-{\psi }_{3-j})-({\psi }_{j}^{\dagger }-{\psi }_{3-j}^{\dagger })\right],\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{q}_{j}={{\rm{\nabla }}}^{2}{\psi }_{j}+\displaystyle \frac{{k}_{{\rm{d}}}^{2}}{2}({\psi }_{3-j}-{\psi }_{j}),\,\,\,\,(j=1,2),\end{eqnarray} \tag{ 5 }$$

where, the subscripts 1 and 2 refer to the upper and lower layers, respectively, and ${q}_{j}({\boldsymbol{x}},t)$ and ${\psi }_{j}({\boldsymbol{x}},t)$ are the potential vorticity and stream function, respectively, of the jth layer (see e.g. Vallis 2006). In addition, ${\boldsymbol{x}}=(x,y),J$ is the horizontal Jacobian, and δ_ij is the Kronecker delta. The flow domain ${ \mathcal D }$ is assumed to be a rectangle of size L_x × L_y. In addition, we assume the presence of dissipation processes: the second term on the right-hand side of (4) represents the hyper-viscosity of order γ, and the third term is the Ekman friction term, which represents the damping of the lower layer's relative vorticity by the lower boundary. ν_j is the hyper-viscosity coefficient for the jth layer and κ is the Ekman friction coefficient. The fourth term represents the Newtonian temperature damping due to radiative cooling towards the prescribed basic states ${\psi }_{j}^{\dagger }$, and α represents the Newtonian cooling coefficient. In addition, the x and y directions are considered to be eastwards and northwards, respectively, for geophysical cases.

We assume that the basic states are horizontally uniform steady flows with a vertical shear:

$$\begin{eqnarray}&&{\psi }_{j}^{\dagger }={\left(-1\right)}^{j}{Uy},\end{eqnarray} \tag{ 6 }$$

where the positive constant U represents the magnitude of velocities for the upper and lower layers. From the thermal wind balance (Vallis 2006),  the presence of basic states with a vertical shear given by (6) implies that there is a temperature difference along the y direction.

We decompose the stream functions into basic state and perturbation components, as follows:

$$\begin{eqnarray}&&{\psi }_{j}({\boldsymbol{x}},t)={\left(-1\right)}^{j}{Uy}+{\psi }_{j}^{\prime} ({\boldsymbol{x}},t),\end{eqnarray} \tag{ 7 }$$

where primed quantities represent perturbation and unknown variables. In addition, we introduce the barotropic and baroclinic stream functions, ψ and τ, as follows:

$$\begin{eqnarray}&&\psi =\displaystyle \frac{1}{2}\left({\psi }_{1}^{{\prime} }+{\psi }_{2}^{{\prime} }\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\tau =\displaystyle \frac{1}{2}\left({\psi }_{1}^{{\prime} }-{\psi }_{2}^{{\prime} }\right).\end{eqnarray} \tag{ 9 }$$

Then, the governing equations (4) become

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}{q}_{\psi } & = & -J(\psi ,{q}_{\psi })-J(\tau ,{q}_{\tau })-U\displaystyle \frac{\partial {{\rm{\nabla }}}^{2}\tau }{\partial x}\\ & & +{\left(-1\right)}^{\gamma -1}\nu {{\rm{\nabla }}}^{2\gamma }{q}_{\psi }+{\left(-1\right)}^{\gamma -1}\displaystyle \frac{{\rm{\Delta }}\nu }{2}{{\rm{\nabla }}}^{2\gamma }({q}_{\psi }-{q}_{\tau })\\ & & -\displaystyle \frac{\kappa }{2}{{\rm{\nabla }}}^{2}(\psi -\tau ),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}{q}_{\tau } & = & -J(\psi ,{q}_{\tau })-J(\tau ,{q}_{\psi })-U\displaystyle \frac{\partial {{\rm{\nabla }}}^{2}\psi }{\partial x}-{{Uk}}_{{\rm{d}}}^{2}\displaystyle \frac{\partial \psi }{\partial x}\\ & & +{\left(-1\right)}^{\gamma -1}\nu {{\rm{\nabla }}}^{2\gamma }{q}_{\tau }+{\left(-1\right)}^{\gamma -1}\displaystyle \frac{{\rm{\Delta }}\nu }{2}{{\rm{\nabla }}}^{2\gamma }({q}_{\tau }-{q}_{\psi })\\ & & -\displaystyle \frac{\kappa }{2}{{\rm{\nabla }}}^{2}(\tau -\psi )+2\alpha {k}_{{\rm{d}}}^{2}\tau ,\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{q}_{\psi }={{\rm{\nabla }}}^{2}\psi ,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{q}_{\tau }={{\rm{\nabla }}}^{2}\tau -{k}_{{\rm{d}}}^{2}\tau ,\end{eqnarray} \tag{ 13 }$$

where q_ψ and q_τ are the barotropic and baroclinic potential vorticities, respectively. Here, we have introduced the notation ν ≡ ν₁ and Δν ≡ ν₂ − ν₁. To simplify both the numerical simulations and theoretical discussions, we adopt doubly periodic boundary conditions for both ψ and τ. These are the fundamental equations considered in this study.

Models similar to (10) and (11) have been used in studies of geostrophic turbulence by many researchers (see e.g. Salmon 1978, 1980, Haidvogel and Held 1980, Larichev and Held 1995, Tung and Orlando 2003, Arbic and Flierl 2004, Thompson and Young 2006, Gkioulekas 2012, 2014) and it is thus a very standard model for geostrophic turbulence studies. One of the notable characteristics of our model is the presence of an asymmetric hyper-viscosity term. The hyper-viscosity has been considered as a useful computational tool for numerical simulations to mimic sub-grid scale dissipations. There are no physical reasons for the asymmetry of the hyper-viscosity. Thus, previous studies have generally set the viscosity coefficients ν₁ and ν₂ to be of equal value.  However, Gkioulekas (2014)  discussed the effect of the asymmetric hyper-viscosity on the Danilov inequality and indicated that the Danilov inequality might be violated due to the presence of the asymmetric hyper-viscosity. To consider this possibility, we adopt the asymmetric hyper-viscosity in our model. Another characteristic of our model is that we neglect the β effect, which represents the dependence of the reference frame's rotation on the y coordinate or that on the Earth's curvature. As stated in section 1, the main purpose of the present study is to investigate the unsolved problem in Gkioulekas (2012, 2014). Thus, we follow the setting up of the problem of his studies, in which the Danilov inequality on an f-plane was investigated. The β effect does not directly affect the sign of (1) because it does not arise explicitly in the evolution equations for the total energy and total potential enstrophy spectra. In addition, the impact of the β effect on the NG spectrum's k^−5/3 region would be negligible, since the so-called Rhines scale (Rhines 1975)  is much larger than the wavenumber range of the k^−5/3 spectrum. Therefore, similar to Gkioulekas (2012, 2014), we neglect the β effect in the present study.

Since we have adopted the doubly periodic boundary conditions, we can expand the stream functions, ψ and τ, as Fourier series:

$$\begin{eqnarray}&&\psi ({\boldsymbol{x}},t)=\displaystyle \sum _{{\boldsymbol{k}}}\hat{\psi }({\boldsymbol{k}},t){{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\tau ({\boldsymbol{x}},t)=\displaystyle \sum _{{\boldsymbol{k}}}\hat{\tau }({\boldsymbol{k}},t){{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{x}}},\end{eqnarray} \tag{ 15 }$$

where ${\boldsymbol{k}}=({k}_{x},{k}_{y})=\left(\tfrac{2\pi }{{L}_{x}}{n}_{x},\tfrac{2\pi }{{L}_{y}}{n}_{y}\right)$ and n_x and n_y are integers. The summation ${\sum }_{{\boldsymbol{k}}}$ is taken over all of the integer vectors (n_x, n_y). The evolution equation for the modal total energy is then as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial { \mathcal E }({\boldsymbol{k}})}{\partial t}={T}_{{ \mathcal E }}({\boldsymbol{k}})+{F}_{{ \mathcal E }}({\boldsymbol{k}})+{D}_{{ \mathcal E }}({\boldsymbol{k}}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{ \mathcal E }({\boldsymbol{k}})=\displaystyle \frac{1}{2}\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}}){| }^{2}+({k}^{2}+{k}_{{\rm{d}}}^{2})| \hat{\tau }({\boldsymbol{k}}){| }^{2}\right],\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{T}_{{ \mathcal E }}({\boldsymbol{k}})={T}_{\psi }^{{\rm{I}}}({\boldsymbol{k}})+{T}_{\psi }^{\mathrm{II}}({\boldsymbol{k}})+{T}_{\tau }^{{\rm{I}}}({\boldsymbol{k}})+{T}_{\tau }^{\mathrm{II}}({\boldsymbol{k}})+{T}_{\tau }^{\mathrm{III}}({\boldsymbol{k}}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{F}_{{ \mathcal E }}({\boldsymbol{k}})=-{k}_{x}{{Uk}}_{{\rm{d}}}^{2}{\mathfrak{Im}}\left[\hat{\psi }({\boldsymbol{k}}){\hat{\tau }}^{* }({\boldsymbol{k}})\right],\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}\begin{array}{rcl}{D}_{{ \mathcal E }}({\boldsymbol{k}}) & = & -\displaystyle \frac{\kappa }{2}{k}^{2}\left(| \hat{\psi }({\boldsymbol{k}})-\hat{\tau }({\boldsymbol{k}}){| }^{2}\right)\\ & & -\nu {k}^{2\gamma }\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}}){| }^{2}+({k}^{2}+{k}_{{\rm{d}}}^{2})| \hat{\tau }({\boldsymbol{k}}){| }^{2}\right]\\ & & -\displaystyle \frac{{\rm{\Delta }}\nu }{2}{k}^{2\gamma }\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}})-\hat{\tau }({\boldsymbol{k}}){| }^{2}+{k}_{{\rm{d}}}^{2}\left(| \hat{\tau }({\boldsymbol{k}}){| }^{2}-{\mathfrak{Re}}[\hat{\psi }({\boldsymbol{k}})\hat{\tau }{\left({\boldsymbol{k}}\right)}^{* }]\right)\right]\\ & & -2\alpha {k}_{{\rm{d}}}^{2}| \hat{\tau }({\boldsymbol{k}}){| }^{2}.\end{array}\end{eqnarray} \tag{ 20 }$$

Here, ${T}_{{ \mathcal E }}({\boldsymbol{k}}),{F}_{{ \mathcal E }}({\boldsymbol{k}})$ and ${D}_{{ \mathcal E }}({\boldsymbol{k}})$ are the transfer function, forcing term and dissipation term for the modal total energy equation, respectively. The derivation of (16)–(20) is summarized in appendix A. The terms that make up (18) are given by (A.4)–(A.7) and (A.9)–(A.14). If we sum (16) over all wavenumber vectors, the term ${T}_{{ \mathcal E }}({\boldsymbol{k}})$, which arises from the governing equations' nonlinear terms, vanishes because the total energy is an inviscid invariant. The sign of ${F}_{{ \mathcal E }}({\boldsymbol{k}})$ is indefinite, but this is the system's only energy source, so we expect ${F}_{{ \mathcal E }}({\boldsymbol{k}})$ to be positive if we sum it over all wavenumber vectors and average the result over some interval. If ${\rm{\Delta }}\nu =0,{D}_{{ \mathcal E }}({\boldsymbol{k}})$ is negative definite. The sign of the asymmetric hyper-viscosity term in the modal total energy equation, which is proportional to Δν in (20), is indefinite. In contrast, the Ekman friction and Newtonian cooling terms (the terms proportional to κ and α, respectively, in (20)) are both negative definite.

In addition, the evolution equation for the modal total potential enstrophy can be obtained similarly:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{ \mathcal Z }({\boldsymbol{k}})={T}_{{ \mathcal Z }}({\boldsymbol{k}})+{F}_{{ \mathcal Z }}({\boldsymbol{k}})+{D}_{{ \mathcal Z }}({\boldsymbol{k}}),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{ \mathcal Z }({\boldsymbol{k}})=\displaystyle \frac{1}{2}\left[{k}^{4}| \hat{\psi }({\boldsymbol{k}}){| }^{2}+{\left({k}^{2}+{k}_{{\rm{d}}}^{2}\right)}^{2}| \hat{\tau }({\boldsymbol{k}}){| }^{2}\right]\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{T}_{{ \mathcal Z }}({\boldsymbol{k}})={k}^{2}\left({T}_{\psi }^{{\rm{I}}}({\boldsymbol{k}})+{T}_{\psi }^{\mathrm{II}}({\boldsymbol{k}})\right)+({k}^{2}+{k}_{{\rm{d}}}^{2})\left({T}_{\tau }^{{\rm{I}}}({\boldsymbol{k}})+{T}_{\tau }^{\mathrm{II}}({\boldsymbol{k}})+{T}_{\tau }^{\mathrm{III}}({\boldsymbol{k}})\right),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{F}_{{ \mathcal Z }}({\boldsymbol{k}})=-{k}_{x}{{Uk}}_{{\rm{d}}}^{4}{\mathfrak{Im}}\left[\hat{\psi }({\boldsymbol{k}}){\hat{\tau }}^{* }({\boldsymbol{k}})\right],\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{rcl}{D}_{{ \mathcal Z }}({\boldsymbol{k}}) & = & -\displaystyle \frac{\kappa }{2}{k}^{2}\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}})-\hat{\tau }({\boldsymbol{k}}){| }^{2}+{k}_{{\rm{d}}}^{2}\left(| \hat{\tau }({\boldsymbol{k}}){| }^{2}-{\mathfrak{Re}}[\hat{\psi }({\boldsymbol{k}}){\hat{\tau }}^{* }({\boldsymbol{k}})]\right)\right]\\ & & -\nu {k}^{2\gamma }\left[{k}^{4}| \hat{\psi }({\boldsymbol{k}}){| }^{2}+{\left({k}^{2}+{k}_{{\rm{d}}}^{2}\right)}^{2}| \hat{\tau }({\boldsymbol{k}}){| }^{2}\right]\\ & & -\displaystyle \frac{{\rm{\Delta }}\nu }{2}{k}^{2\gamma }\left[| {k}^{2}\hat{\psi }({\boldsymbol{k}})-({k}^{2}+{k}_{{\rm{d}}}^{2})\hat{\tau }({\boldsymbol{k}}){| }^{2}\right]\\ & & -2\alpha {k}_{{\rm{d}}}^{2}({k}^{2}+{k}_{{\rm{d}}}^{2})| \hat{\tau }({\boldsymbol{k}}){| }^{2}.\end{array}\end{eqnarray} \tag{ 25 }$$

Here, ${T}_{{ \mathcal Z }}({\boldsymbol{k}}),{F}_{{ \mathcal Z }}({\boldsymbol{k}})$ and ${D}_{{ \mathcal Z }}({\boldsymbol{k}})$ are the transfer function, forcing term and dissipation term for the modal total potential enstrophy equation, respectively. The ratio of the forcing terms for the modal total potential enstrophy and modal total energy, ${F}_{{ \mathcal Z }}({\boldsymbol{k}})/{F}_{{ \mathcal E }}({\boldsymbol{k}})$, is exactly ${k}_{{\rm{d}}}^{2}$. Unlike for the modal total energy equation, here the sign of the asymmetric hyper-viscosity term is definite and the sign of the Ekman friction term is indefinite. As we will show later, the fact that the signs of the asymmetric hyper-viscosity term for the modal total energy equation and Ekman friction term for the modal total potential enstrophy equation are not definite is one of the potential causes of the Danilov inequality violation in this system.

In this section, we conduct direct numerical simulations of turbulence governed by (10) and (11). First, we outline the numerical method, initial conditions and calculation resolution. The values used for the parameters in (10) and (11) were essentially the same as those used in Larichev and Held (1995) and Arbic and Flierl (2004) (hereinafter LH95 and AF04, respectively).

The pseudo-spectral method was used with double-precision arithmetic at a resolution of N² (the number of grid points in the square computational domain [0, L] × [0, L]), and the truncation wavenumber k_T was chosen according to the two-thirds dealiasing rule. The domain size was L = 2π. Time integration was performed using the Runge–Kutta–Gill scheme with a fixed time step Δt, ensuring that Δt values used satisfied the Courant–Friedrichs–Lewy condition.

The initial conditions were set by generating uniform random numbers in the interval [0, 2π] for the phases of $\hat{\psi }({\boldsymbol{k}})$ and $\hat{\tau }({\boldsymbol{k}})$. The initial amplitudes of $\hat{\psi }({\boldsymbol{k}})$ and $\hat{\tau }({\boldsymbol{k}})$ were normalized so that the barotropic and baroclinic energy spectra, E_ψ(k) and E_τ(k), were independent of k:

$$\begin{eqnarray}&&\,{E}_{\psi }\left(k\right)\equiv \left\langle \displaystyle \frac{1}{2}{k}^{2}{\left|\hat{\psi }\left({\boldsymbol{k}}\right)\right|}^{2}\right\rangle =3.0\times {10}^{-8},{E}_{\tau }\left(k\right)\equiv \left\langle \displaystyle \frac{1}{2}\left({k}^{2}+{k}_{{\rm{d}}}^{2}\right){\left|\hat{\tau }\left({\boldsymbol{k}}\right)\right|}^{2}\right\rangle =3.0\times {10}^{-8}.\end{eqnarray} \tag{ 40 }$$

First, the system was spun up at a resolution of N = 256 until t = t₁, and then, numerical integration was performed at a higher resolution of N = 512 until t = t₂. The results of the higher resolution calculations were analyzed and averaged over the interval between ${t}_{2}-{t}_{\mathrm{ave}}$ and t₂ to improve the statistical convergence of the results. Since, as shown in the previous section, the Newtonian damping term in (11) does not contribute to violating the Danilov inequality, it is omitted in the following discussion, i.e. α = 0.

Scaling the spatial coordinates by ${k}_{{\rm{d}}}^{-1}$, the time by ${(2{{Uk}}_{{\rm{d}}})}^{-1}$, the stream functions by $2{{Uk}}_{{\rm{d}}}^{-1}$ and the potential vorticities by $2{{Uk}}_{{\rm{d}}}$ leads to the governing equations being controlled by the following three dimensionless parameters:

$$\begin{eqnarray}&&\displaystyle \frac{2{{Uk}}_{{\rm{d}}}}{\kappa },\,\,\displaystyle \frac{2U}{\nu {k}_{{\rm{d}}}^{2\gamma -1}},\,\,\displaystyle \frac{2U}{{\rm{\Delta }}\nu {k}_{{\rm{d}}}^{2\gamma -1}}.\end{eqnarray} \tag{ 41 }$$

The parameter $(2{{Uk}}_{{\rm{d}}})/\kappa$, which Arbic and Flierl (2003)  call the throughput, is particularly significant because it controls the partition of the total energy into barotropic and baroclinic components or into kinetic energy and available potential energy. For example, when the throughput is small, the total energy is dominated by the available potential energy and the system behaves as if it is equivalent barotropic. In contrast, when the throughput is large, the total energy is dominated by the kinetic energy and the system behaves as if it is barotropic.

To resolve the high-wavenumber dynamics in the k > k_d region, our region of interest, and ensure that there is sufficient spectral room for the upward energy cascade in k < k_d, the region of primary interest for most geostrophic turbulence studies, we set k_d = 10 for all simulations, even though LH95 used k_d = 50 for their k_T = 128 resolution simulations. In addition, we fixed U = 2.5 × 10⁻², making our Uk_d value equivalent to that used in LH95. By changing the value of κ, we could vary the throughput from 1.25 × 10⁻² to 25, an equivalent range to that studied in AF04. Note in passing that we could not obtain from numerical simulations the statistical steady state of the system when we adopted $\tfrac{2{{Uk}}_{{\rm{d}}}}{\kappa }=50$. Since these experiments are those for changing the throughput value, we shall refer to these as the TH experiments. Note that the estimated throughput in the TO03 experiments was nearly 27.3 (see appendix B).

The hyper-viscosity parameters used in the present study were similar to those used in LH95. They used γ = 4 and ν = 5.12 × 10⁻¹⁶ for k_T = 128 resolution simulations. We decided to maintain the e-folding time for the hyper-viscosity at the truncation wavenumber at the value used in LH95, so we set ν = 5.12 × 10⁻¹⁶ × (128/k_T)^2γ. The numerical parameters used in the TH experiments are summarized in table 1.

We made sure to resolve the dissipation ranges for the total energy and total potential enstrophy spectra, because this affects the accuracy of the inertial range statistics (see e.g. Smith 2004, Watanabe and Gotoh 2007).  If the viscosity coefficients ${\nu }_{1}$ and ν₂ are not large enough to remove the total energy and total potential enstrophy transferred to higher wavenumbers, their spectra will exhibit transitions to shallower slopes near the truncation wavenumber. One of our criteria for resolving the dissipation ranges is that the dissipation wavenumbers have to be smaller than the truncation wavenumber k_T. By dimensional analysis, the dissipation wavenumber k_ε, which is based on the rate of total energy dissipation due to hyper-viscosity (ε_ν) and is known as the energy dissipation wavenumber, can be defined as

$$\begin{eqnarray}&&{k}_{\varepsilon }\equiv {\left(\displaystyle \frac{{\varepsilon }_{\nu }}{{\tilde{\nu }}^{3}}\right)}^{1/(6\gamma -2)},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\varepsilon }_{\nu }\equiv {\left(-1\right)}^{\gamma -1}\displaystyle \sum _{i=1}^{2}\overline{{\nu }_{i}{\psi }_{i}{{\rm{\nabla }}}^{2\gamma }{q}_{i}}.\end{eqnarray} \tag{ 43 }$$

Here, $\overline{A}$ indicates the spatial average of a given quantity A, and $\tilde{\nu }$ is the average viscosity coefficient, i.e. $\tilde{\nu }\equiv ({\nu }_{1}+{\nu }_{2})/2$. The enstrophy dissipation wavenumber k_η, which is similarly based on the rate of total potential enstrophy dissipation due to hyper-viscosity (η_ν), can be defined as

$$\begin{eqnarray}&&{k}_{\eta }\equiv {\left(\displaystyle \frac{{\eta }_{\nu }}{{\tilde{\nu }}^{3}}\right)}^{1/6\gamma },\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{\eta }_{\nu }\equiv {\left(-1\right)}^{\gamma }\displaystyle \sum _{i=1}^{2}\overline{{\nu }_{i}{q}_{i}{{\rm{\nabla }}}^{2\gamma }{q}_{i}}.\end{eqnarray} \tag{ 45 }$$

Here, table 2 gives the values of ε_ν, η_ν, the rates of total energy and total potential enstrophy dissipations due to Ekman friction, ε_κ and η_κ, where ε_κ and η_κ are defined as

$$\begin{eqnarray}&&{\varepsilon }_{\kappa }\equiv -\overline{\kappa {\psi }_{2}{{\rm{\nabla }}}^{2}{\psi }_{2}},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\eta }_{\kappa }\equiv \overline{\kappa {q}_{2}{{\rm{\nabla }}}^{2}{\psi }_{2}}.\end{eqnarray} \tag{ 47 }$$

These results show that the bulk of the total energy is dissipated due to Ekman friction, independent of the throughput. In contrast, hyper-viscosity is the main cause of potential enstrophy dissipation for high throughputs, but as the throughput decreases, the contributions of hyper-viscosity and Ekman friction become broadly equal. From the evolution equations of the potential enstrophy in the individual layers, the potential enstrophies in the individual layers are inviscid invariants. However, the present system is forced by the baroclinic instability of basic flow. Then, the same amount of potential enstrophy, ${\left(-1\right)}^{j+1}\tfrac{1}{2}{{Uk}}_{{\rm{d}}}^{4}{\int }_{{ \mathcal D }}{\psi }_{3-j}\tfrac{\partial {\psi }_{j}}{\partial x}\,{\rm{d}}x{\rm{d}}y$, is supplied to the individual layers. We can interpret these results in terms of the energy and potential enstrophy balances in the individual layers. Since the viscosity is the only dissipation mechanism of the potential enstrophy in the upper layer, the potential enstrophy in the upper layer is dissipated by the hyper-viscosity, independent of the throughput. In contrast, the Ekman friction and hyper-viscosity are the dissipation mechanisms in the lower layer. Then, the potential enstrophy in the lower layer is mainly dissipated by the hyper-viscosity for high throughputs. As the throughput decreases, the Ekman friction tends to dominate more and more. In the case of the lowest throughput realized in the present simulations, the potential enstrophy in the lower layer is mainly dissipated by the Ekman friction. In contrast, the upper layer energy is transferred to the lower layer and the upper and lower layers' energies are dissipated by the Ekman friction, independent of the throughput.

Figure 1 shows how the energy and enstrophy dissipation wavenumbers depend on the throughput, indicating that both dissipation wavenumbers are less than the truncation wavenumber k_T = 170.

Zoom In Zoom Out Reset image size

Another criterion for measuring the dissipation range is that of the spectrum of total energy dissipation due to hyper-viscosity, which is defined by

$$\begin{eqnarray}\begin{array}{rcl}{D}_{\nu }(k) & = & \left\langle \Space{0ex}{3.0ex}{0ex}\langle \right.\nu {k}^{2\gamma }\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}}){| }^{2}+({k}^{2}+{k}_{{\rm{d}}}^{2})| \hat{\tau }({\boldsymbol{k}}){| }^{2}\right]\\ & & \left.+\displaystyle \frac{{\rm{\Delta }}\nu }{2}{k}^{2\gamma }\left[{k}^{2}| \hat{\psi }({\boldsymbol{k}})-\hat{\tau }({\boldsymbol{k}}){| }^{2}+{k}_{{\rm{d}}}^{2}\left(| \hat{\tau }({\boldsymbol{k}}){| }^{2}-{\mathfrak{Re}}[\hat{\psi }({\boldsymbol{k}})\hat{\tau }{\left({\boldsymbol{k}}\right)}^{* }]\right)\right]\right\rangle ,\end{array}\end{eqnarray} \tag{ 48 }$$

which decays sufficiently rapidly in the wavenumber range k_ε ≲ k < k_T. Figure 2 shows the dissipation spectra for four selected experiments, TH2, TH5, TH7 and TH9. Here, the dissipation spectra take their maximum values around k_ε and then decay rapidly as the wavenumber increases further. Although we have only shown the dissipation spectra for some of the experiments, the dissipation spectra for the other TH experiments behave similarly. We can thus conclude that our experiments resolve the dissipation range well.

Zoom In Zoom Out Reset image size

This study focuses on the energy and enstrophy transfers in wavenumber space. To ensure the validity of our numerical simulations and analysis, we also verified that the transfer functions that make up (18) satisfy conservation laws,  and this is discussed in appendix C.

Figure 3 shows snapshots of the barotropic and baroclinic potential vorticities for the same four TH experiments as in figure 2. As the throughput decreases, the barotropic potential vorticity fields change drastically from circular to more ribbon-like shapes, while the baroclinic potential vorticity fields change from circular to more blob-like shapes. On the other hand, the kinetic and available potential energy fields have circular structures for large throughputs, but become ribbon-shaped (surrounding the available potential energy extrema) and blob-shaped, respectively, as the throughput decreases, as pointed out by AF04 (figures are not shown). In addition, the structures of the barotropic and baroclinic stream functions become similar as the throughput decreases (figures are not shown ). Since the lower layer's stream function is ψ₂ = ψ − τ, this indicates that its motion is reduced because of the bottom friction as the throughput decreases.

Zoom In Zoom Out Reset image size

Figure 4 shows the dependence of the ratios of kinetic energy K to total energy E and barotropic energy E_ψ to total energy E on the throughput. Here, K, E_ψ and E are defined by

$$\begin{eqnarray}&&K=-\displaystyle \frac{1}{2}\overline{\left(\psi {{\rm{\nabla }}}^{2}\psi +\tau {{\rm{\nabla }}}^{2}\tau \right)},\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{E}_{\psi }=-\displaystyle \frac{1}{2}\overline{\psi {{\rm{\nabla }}}^{2}\psi }=\displaystyle \sum _{k}{E}_{\psi }(k),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&E=-\displaystyle \frac{1}{2}\overline{\left[\psi {{\rm{\nabla }}}^{2}\psi +\tau ({{\rm{\nabla }}}^{2}-{k}_{{\rm{d}}}^{2})\tau \right]}=\displaystyle \sum _{k}E(k).\end{eqnarray} \tag{ 51 }$$

Note that the difference between E_ψ/E and K/E is equal to K_τ/E, where K_τ is the baroclinic kinetic energy. The potential energy dominates for small throughputs, but the kinetic energy increases with throughput. In addition, the kinetic energy is dominated by the barotropic component for large throughputs. Taken together, those results indicate that our numerical experiments cover a wide range of two-layer geostrophic turbulence dynamics.

Zoom In Zoom Out Reset image size

Although this paper's main aim is to examine whether or not the Danilov inequality holds for two-layer geostrophic turbulence, in this section we shall also investigate fundamental turbulence properties such as the total energy spectrum, total energy and total potential enstrophy fluxes. First, figure 5 gives the total energy spectra for all of the TH experiments, showing that they follow a power law of the form E(k) ∼ k^−δ with 3 ≲ δ ≲ 4 in the wavenumber range k_d ≲ k ≲ k_ε and smoothly decay as k increases further. Note that there is no transition from E(k) ∼ k^−δ to E(k) ∼ k^−5/3. The exponent δ is close to 3 for larger throughputs and increases as the throughput decreases. The exponents are very close to those of the enstrophy inertial range spectra for forced-dissipative 2D NS turbulence.

Zoom In Zoom Out Reset image size

Now, we develop this analysis of the total energy spectra further. TO03 defined the transition wavenumber k_trans as

$$\begin{eqnarray}&&{k}_{\mathrm{trans}}\equiv \sqrt{\displaystyle \frac{{\eta }_{\nu }}{{\varepsilon }_{\nu }}}.\end{eqnarray} \tag{ 52 }$$

They stated that the enstrophy inertial range occurs in the wavenumber range k < k_trans and the energy inertial range appears in the wavenumber range k > k_trans. Using (42) and (44), the transition wavenumber is related to the dissipation wavenumbers via

$$\begin{eqnarray}&&{k}_{\mathrm{trans}}={\left(\displaystyle \frac{{k}_{\eta }}{{k}_{\varepsilon }}\right)}^{3\gamma }{k}_{\varepsilon }.\end{eqnarray} \tag{ 53 }$$

For the energy inertial range to exist, the transition wavenumber k_trans must be less than the energy dissipation wavenumber k_ε and (53) indicates that the transition wavenumber is less than the energy dissipation wavenumber if k_ε > k_η. However, figure 1 shows that k_η was very similar to but slightly larger than k_ε in all of the TH experiments, which means that the energy inertial range was not present in these experiments. This is consistent with the total energy spectrum analysis. Note that it is impossible to control the values of k_ε and k_η a priori, since our model is internally forced by an unstable basic flow. Furthermore, as shown in the previous subsection, our numerical experiments cover a wide regime of two-layer geostrophic turbulence dynamics. Thus, the relation k_ε < k_η obtained by the numerical simulations should be interpreted as one of the intrinsic properties of the present system.

Next, we examine the total potential enstrophy and total energy fluxes, which are shown in figure 6, normalized by their dissipation rates due to hyper-viscosity. Again, results are only presented for four TH experiments, but the results of all of the TH experiments were analyzed. In addition, these figures only show the fluxes in the range k > k_d, since this is our wavenumber range of interest. The total potential enstrophy flux, which is equal to the total potential enstrophy dissipation rate η_ν, exists over a wide range of wavenumbers less than k_ε. Likewise, the total energy flux, which is equal to the total energy dissipation rate ε_ν, is present over a similar but narrower wavenumber range. In this sense, there is the so-called hidden energy cascade, as suggested by TO03. We note in passing that the fluxes in the range k ≲ k_d are consistent with classical theories of two-layer geostrophic turbulence, as shown in figure 2 of Thompson and Young (2006).

Zoom In Zoom Out Reset image size

Finally, we examine whether or not the Danilov inequality holds by computing the flux difference (1) and considering their sign. Figure 7 shows that the flux difference is negative over the whole wavenumber range 1 ≤ k ≤ k_T and this was confirmed by checking the numerical flux difference values. Thus, the Danilov inequality held in these experiments. As shown in (39), functional forms of the forcing spectra F_E(k) and F_Z(k) are key factors in determining the sign of the flux difference. The numerical results show that those spectra are significantly distributed in the range k ≤ k_d and are mainly positive there. (figures are not shown.) Moreover, the relation ${k}^{2}{F}_{E}(k)\lt {F}_{Z}(k)$ is held in the same wavenumber range, since ${F}_{Z}(k)={k}_{{\rm{d}}}^{2}{F}_{E}(k)$. As a result, the forcing spectra make the flux difference positive there, as shown in figure 8, which shows the contributions of the forcing spectra on the flux difference, ${\left[{k}^{2}{{\rm{\Pi }}}_{E}(k)-{{\rm{\Pi }}}_{Z}(k)\right]}_{{\rm{forc}}}\equiv -{\int }_{k}\infty \left[{k}^{2}{F}_{E}({k}^{{\rm{{\prime} }}})-{F}_{Z}({k}^{{\rm{{\prime} }}})\right]\,{\rm{d}}{k}^{{\rm{{\prime} }}}$. In contrast, the dissipation spectra due to the Ekman friction are significant in the same wavenumber range, compensate the contribution of the forcing spectra on the flux difference and make the flux difference negative.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have shown that the slopes of the simulated total energy spectra in the range k > k_d are very close to those of the enstrophy inertial range and show no evidence of a transition to those of the energy inertial range. Furthermore, the Danilov inequality held for the whole wavenumber range. However, it is natural to ask whether or not either of these results could have been affected by the simulation resolution. To verify that these results were independent of the simulation resolution, we doubled the resolution and found no significant changes.

As described in section 3, asymmetric hyper-viscosity may lead to the Danilov inequality being violated. To investigate this issue, we performed six additional numerical simulations under the same conditions as for experiment TH2 but with different values of Δν. We shall refer to these experiments as the HV experiments, and the values of Δν used are listed in table 3. Table 4 lists the total energy and total potential enstrophy dissipation rates, as well as the total energies obtained for the HV experiments.

Similar to the TH experiments, both the energy and enstrophy dissipation wavenumbers are lower than the truncation wavenumber. In addition, the enstrophy dissipation wavenumber is higher than the energy dissipation wavenumber, and the transition wavenumber lies in the dissipation range, as found by the TH experiments. This therefore implies that the energy inertial range suggested by TO03 did not appear for k_d ≲ k < k_T in these experiments.

Figure 9 shows the dissipation spectra for all of the HV experiments, together with the spectrum for experiment TH2 for comparison. This shows that the dissipation spectra for TH2, HV4, HV5 and HV6 take their maximum values around the dissipation wavenumber k_ε and then decay sufficiently rapidly as the wavenumber increases further. In contrast, although the dissipation spectrum for HV3 also has a local extreme around k_ε, it starts to increase rapidly near the truncation wavenumber and takes its global maximum value there. Furthermore, the dissipation spectra for HV1 and HV2 do not even have local maxima around k_ε. We can therefore conclude that two of the experiments at least (HV1 and HV2) could not resolve the dissipation range.

Zoom In Zoom Out Reset image size

Figure 10 shows the total energy spectra for all of the HV experiments, together with experiment TH2 for comparison. In all cases, these spectra follow power laws of the form $E(k)\sim {k}^{-\delta }$ with δ ≃ 3 in the wavenumber range k_d ≲ k ≲ k_ε. The spectra for HV4, HV5 and HV6 then decay smoothly as the wavenumber increases further, almost converging on the TH2 spectrum. In contrast, the slopes of the spectra for HV1, HV2 and HV3 become shallower for k ≳ k_ε as Δν decreases. In the cases of HV1 and HV2, we attribute this behavior to insufficient resolution in the dissipation range. As with the TH experiments, there is no evidence of the transition from E(k) ∼ k^−δ to E(k) ∼ k^−5/3 seen by TO03.

Zoom In Zoom Out Reset image size

Finally, we consider the sign of the flux difference (1), as shown in figure 11. Again, the sign is negative in the entire range 1 ≤ k ≤ k_T, i.e. the Danilov inequality held over the entire wavenumber range, even for the HV experiments.

Zoom In Zoom Out Reset image size