will refer to this system as the Boussinesq model, or BQ model. For $\mathbf{F}=0$ and $\nu=0$, the dispersion relation obtained by linearizing these equations is $$\label{rel_disp_bous}
    \omega = \frac{c_0 k }{\sqrt{1 + \frac{h^2_0 k^2}{3}}} ,$$ which, up to the fourth order, coincides with Eq. .

Note that there are other choices for the extra term in Eq.  that result in many formulations of the Boussinesq model, all compatible up to fourth order in a Taylor expansion in terms of $h_0 k$ [@whitham_linear_1974]. The formulation we use here was employed in previous studies of wave turbulence [@onorato_four-wave_2008], and is also easy to solve numerically using pseudospectral methods by writing Eq.  as $$\frac{\partial \mathbf{u'}}{\partial t} = -(\mathbf{u} \cdot {\boldsymbol{\nabla}}) \mathbf{u}  -g
    {\boldsymbol{\nabla}}h + \frac{\nu}{h} {\boldsymbol{\nabla}}\cdot ( h {\boldsymbol{\nabla}}\mathbf{u}) +
    \mathbf{F},
\label{eq:Helmholtz}$$ where $\mathbf{u'} = {\cal H}\mathbf{u}$, and where ${\cal H} = (1-h_0^2 {\nabla^2}/3)$ is the Helmholtz operator. This operator can be easily inverted in Fourier space [@mininni_numerical_2005; @mininni_numerical_2005-1], and the resulting equations can be efficiently solved by means of pseudospectral codes. It is interesting that the same operator appears in Lagrangian-averaged models [@foias_navier-stokes-alpha_2001]. In these models, and in regularized versions of the shallow water equations [@camassa_integrable_1993], it introduces dispersion that results in 