Numerical methods

Table 1: Numerical parameters.

	modes	$f$ ( $\times 10^{-4}$rad/sec)	$L_{\mathrm{v}}$ ($\times 10$kg/m$^3$)	forcing	initial condition
I	$512^2 \times 256$	$\sqrt{2}/3$	2.7	none	GM
II	$1024^2 \times 512$	$\sqrt{2}/3$	2.7	$\omega \sim 3f$	white noise
III	$256^3$	$0.25$	5	$\vert\bldk\vert^2 \! + m^2 \! \leq 6^2$	white noise
IV	$256^3$	$1$	5	$\vert\bldk\vert^2 \! + m^2 \! \leq 6^2$	white noise
V	$512^2 \times 256$	$\sqrt{2}/3$	2.7	none	GM w/o long waves

In direct numerical simulations, we add external forcing and hyper-viscosity to the canonical equation (5), as forcing and damping are needed to achieve statistically steady states. Thus the resulting dynamic equation is given by

$$\frac{\upartial a(\bldp)}{\upartial t} = - \mathrm{i}$$ (17)

$$end{tex2html_deferred}, \omega(\bldp) a(\bldp) + {\cal N}(a(\bldp)) + F(\bldp) - D(\bldp) a(\bldp)$$ (18)

$$end{tex2html_deferred}, .$$ (19)

In the simulations, the linear terms, which are a linear dispersion term and a dissipation term, $-D(\bldp) a(\bldp)$, are explicitly calculated. The nonlinear terms, ${\cal N}(a(\bldp))$, derived from the nonlinear parts of the canonical equation (5) with Hamiltonian (6) are obtained by a pseudo-spectral method with the phase shift under the periodic boundary conditions for all three directions. The external forcing, $F(\bldp)$, is implemented by fixing the amplitudes of several small wavenumbers to be constant in time. The dissipation is also added as $D(\bldp) = D_{\mathrm{h}} \vert\bldk\vert^8 + D_{\mathrm{v}} \vert m\vert^4$. Here, $D_{\mathrm{h}}$ and $D_{\mathrm{v}}$ are chosen so that the dissipation is effective for wavenumbers larger than the half of maximum wavenumbers. The wavenumbers are discretized as $\bldp = (2\upi/ L_{\mathrm{h}} \bldk, \: 2\upi/ L_{\mathrm{v}} m)$, where $L_{\mathrm{h}}$ and $L_{\mathrm{v}}$ are horizontal periodic length and vertical period in the isopycnal coordinates, and $\bldk$ and $m$ are integer-valued wavenumbers afterward. Time-stepping is implemented with the fourth-order Runge-Kutta method. In all the simulations, the buoyancy frequency and horizontal period are fixed at $N_0 = 10^{-2}$rad/sec and $L_{\mathrm{h}}=10^5$m, respectively.

We perform a series of five numerical experiments that are listed in Table 1. The total energy per unit periodic box of all the numerical experiments except Run II is around $3 \times 10^3$J/(kg $\cdot$ m$^2$) which is characteristic of real oceans. The total energy in the Run II is around $1.2 \times 10^3$J/(kg $\cdot$ m$^2$).

Two-dimensional energy spectra are measured in the shell of radius $k$ as

$$E(k, \vert m\vert) = \int\limits_{k- 1/2 \leq \vert\bldk^{\prime}... ...hrm{d}\bldk^{\prime} \sum_{s=\pm 1} \omega \vert a(\bldk^{\prime}, s m)\vert^2$$ (20)

$$end{tex2html_deferred}, .$$ (21)

Integrated energy spectra are defined as $\overline{E}_{\mathrm{int}}(k) = \int \mathrm{d}m E(k,\vert m\vert)$ and $\overline{E}_{\mathrm{int}}(\vert m\vert) = \int \mathrm{d}k E(k,\vert m\vert)$. Cross-sectional spectra, $E_m(k)$ and $E_k(\vert m\vert)$, are obtained from the two-dimensional energy spectra as a function of isopycnal wavenumbers $k$ along a certain density wavenumber $m$ and as a function of density wavenumbers $m$ along a certain isopycnal wavenumber $k$, respectively.

In general in real oceans only integrated spectra, $\overline{E}_{\mathrm{time}}(\omega)$ obtained from time series and $\overline{E}_{\mathrm{vertical}}(m)$ from vertical series, can be measured. Then the two-dimensional spectrum, $E(\omega,m)$, is composed from integrated spectra. To do it, the assumption of separability of two-dimensional spectra is made. Namely, the functional form, $E(\omega,m) \propto \overline{E}_{\mathrm{time}}(\omega) \overline{E}_{\mathrm{vertical}}(m)$, is assumed. As we will see below in our numerical simulations, our resulting spectra are not always separable. Such non-separability is also observed in the real oceans (Polzin, 2007). Note that the power-law exponents of one-dimensional spectra and those of two-dimensional spectra do not always coincide in our direct numerical simulation, as explained below. We do not have any means to see whether this observation will also hold in the real ocean.

Figure 1: Run I. Cross-sectional spectra $E_k(\vert m\vert)$ of GM spectrum as the initial condition (left) and of energy spectrum after about 35 hours (right). Significant differences in the region $16 < \vert m\vert < 64$ indicates that GM spectrum is not statistically steady.