---
abstract: 'The *equations of Boussinesq approximation* (EBA) for an incompressible and inhomogeneous in density fluid are analyzed from a viewpoint of the asymptotic theory. A systematic scaling shows that there is an infinite number of related asymptotic models. We have divided them into three classes: ‘poor’, ‘reasonable’ and ‘good’ Boussinesq approximations. Each model can be characterized by two parameters $q$ and $k$, where $q =1, 2, 3, \dots$ and $k=0, \pm 1, \pm 2,\dots$. Parameter $q$ is related to the ‘quality’ of approximation, while $k$ gives us an infinite set of possible scales of velocity, time, viscosity, *etc.* Increasing $q$ improves the quality of a model, but narrows the limits of its applicability. Parameter $k$ allows us to vary the scales of time, velocity and viscosity and gives us the possibility to consider any initial and boundary conditions. In general, we discover and classify a rich variety of possibilities and restrictions, which are hidden behind the routine use of the Boussinesq approximation. The paper is devoted to the multiplicity of scalings and related restrictions. We do not study any particular solutions and particular failures of EBA.'
author:
- 'V.A.VladimirovandN.Al-Salti'
date: November 12th 2016
title: Many Faces of Boussinesq Approximations
---

Introduction \[sect01\]
=======================

In the *equations of Boussinesq approximation* (EBA) the density variations are neglected in the inertial terms of the equations of motion. EBA are actively employed to describe the flows of a stratified fluid, as well as convective flows. The use of EBA is so common that it is often accepted as a ‘starting’ or ‘original’ system of governing equations in many studies. A number of papers and books gives purely physical justification of the EBA (or does not give any justifications at all), see *e.g.* [@GershuniC; @Craik; @Moffatt; @Lighthill; @YihF; @Grimshaw; @DrazinStab]. In this note we present a general analysis of EBA from a viewpoint of the asymptotic theory. The asymptotic approach to EBA has already been exploited in few papers. In [@Ver; @Zey] the EBA for a compressible fluid are analyzed. The aim of the authors was to ‘filter out’ acoustic waves with the Mach number used as a small parameter. [@Long] analysed the EBA for an incompressible fluid, however, he used it alone with the long-wave approximation. The generalization to magneto-hydrodynamic flows has been presented by [@Dave]. In our paper we avoid all those very interesting generalizations and concentrate our attention at the simplest case of an incompressible fluid.

We derive EBA with the use of the same ‘asymptotic spirit’ as, say, in the derivations of KDV equation, Boussinesq equation for surface waves, Schrodinger equation in various areas of fluid dynamics. Our analysis shows the multiplicity of the scalings and, simultaneously, their internal restrictions. For example, one such restriction shows that if the velocity is small (of order $\varepsilon$) then the time-scale must be of order $1/\varepsilon$.

The aim of this paper is to establish and analyse the multiplicity of scalings leading to EBA, the quality of EBA models, and their restrictions. The sensitive point to consider is: the authors, who consider EBA as a ‘starting’ system of governing equations, introduce a scaling required for their particular problems. At the same time EBA require their own scaling, which is not unique but rather restrictive. It is apparent that in a systematic approach these two independently introduced scalings should be incorporated into one scaling, which simultaneously produces EBA and other desired properties of the equations (long waves, multi-scale, *etc.*). The possibility of such an incorporation should be analysed separately, especially for multi-scale theories. Our analysis is concentrated at the general conditions of applicability and quality of EBA with the message that EBA should not be used ‘blindly’, as an original governing equation. In our view, the discovered multiplicity and conditions of applicability rise interesting questions to some existing models and results on nonlinear waves, nonlinear stability, and solitons in a continuously stratified flows, see [@Craik; @Grimshaw1]. Another area, where our results could be implemented, is vibrational convection, see [@GershuniV]. Contrary to [@Long], we do not target the particular solutions of EBA and the failures of EBA.

General scaling and basic equations
===================================

The dimensional equations describing three-dimensional flows of an incompressible stratified viscous fluid are $$\begin{aligned}
&&\rho^* {\boldsymbol{a}}^*=\rho^*[{\boldsymbol{u}}^*_{t^*}+({\boldsymbol{u}}^*\cdot\nabla^*){\boldsymbol{u}}^*]=-\nabla^* p^*-\rho^*\nabla^*\Phi^*+\mu^*\Delta^*{\boldsymbol{u}}^*\label{maineq-dim}\\
&& \rho^*_{t^*}+ ({\boldsymbol{u}}^*\cdot\nabla^*)\rho^*=0,\quad \nabla^*\cdot{\boldsymbol{u}}^*=0\label{maineq-dim1}\\
&& \mu^*={\mathrm{const}},\quad\nabla^*\equiv (\partial/\partial x_1^*, \partial/\partial x_2^*, \partial/\partial x_3^*),\quad \Delta^*\equiv\nabla^{*2}\nonumber\end{aligned}$$ where $t^*$ and ${\boldsymbol{x}}^*=({\boldsymbol{x}}_1^*,{\boldsymbol{x}}_2^*, {\boldsymbol{x}}_3^*)$ are time and cartesian coordinates; ${\boldsymbol{a}}^*({\boldsymbol{x}}^*,t^*)$, ${\boldsymbol{u}}^*({\boldsymbol{x}}^*,t^*)$, $\rho^*({\boldsymbol{x}}^*,t^*)$, $p^*({\boldsymbol{x}}^*,t^*)$, $\mu^*$ are the fields of acceleration, velocity, density, pressure, and constant viscosity correspondingly; $\Phi^*({\boldsymbol{x}}^*)$ is a given potential of external force, such that the acceleration of gravity is ${\boldsymbol{g}}^*=-\nabla^*\Phi^*$; the subscripts of independent variables denote partial derivatives and asterisks stand for dimensional variables. The exact solution of these equations is an equilibrium of a homogeneous fluid: $$\begin{aligned}
 && \rho^*({\boldsymbol{x}}^*,t^*)=R={\mathrm{const}},\quad {\boldsymbol{u}}^*({\boldsymbol{x}}^*,t^*)=0,\quad p^*({\boldsymbol{x}}^*,t^*)=-R\,\Phi^*({\boldsymbol{x}}^*)+{\mathrm{const}}\label{mainsoln}
 \end{aligned}$$ It is also called *a reference solution*. Below we consider the class of solutions $$\begin{aligned}
\rho^*=\rho^*({\boldsymbol{x}}^*,t^*),\quad {\boldsymbol{u}}^*={\boldsymbol{u}}^*({\boldsymbol{x}}^*,t^*), \quad p^*=p^*({\boldsymbol{x}}^*,t^*)\label{genflow}\end{aligned}$$ which represent small deviations from , characterised by two small parameters $$\begin{aligned}
 && \delta_1\equiv\frac{||\rho^*-R||}{R}\ll 1,\quad \delta_2\equiv\frac{||{\boldsymbol{a}}^*||}{||{\boldsymbol{g}}^*||}\sim \frac{||{\boldsymbol{u}}^*_{t^*}||}{||\nabla^*\Phi^*||}\sim \frac{||({\boldsymbol{u}}^*\cdot\nabla^*){\boldsymbol{u}}^*||}{||\nabla^*\Phi^*||}\ll 1\label{physics}\end{aligned}$$ where the particular choice of mathematical norms is not important; for example they can be $C$-norms. Two key EBA assumptions are, as a rule, accepted by all authors, often implicitly. The list of characteristic parameters for solutions is taken as $$\begin{aligned}
{L},\quad R,\quad{G}, \quad \delta_1, \quad \delta_2\label{params0}\end{aligned}$$ where ${L}$ is a characteristic length, ${G}$ is a characteristic value for $|{\boldsymbol{g}}^*|$. We also take $\mu^*=R\nu^*=const$, where $\nu^*$ is a constant value of dimensional kinematic viscosity. Parameter $\nu^*$ is not included into the list , since our intention is to consider an inviscid fluid $\nu^*=0$ as a primary case and a viscous fluid $\nu^*\neq 0$ only as an indication of related possibilities. For the latter case $\nu^*$ must be included into the list , while the dependent viscous length scale is $L_\nu=(\nu^{*2}/G)^{1/3}$. Then, we consider $L=L_\nu$, in order to avoid a consideration of boundary layers, *etc.*

The imposed conditions appear as two independent small parameters $\delta_1$ and $\delta_2$. They must be expressed as a path $(\delta_1(\varepsilon), \delta_2(\varepsilon))$ in the plane $(\delta_1,\delta_2)$ with the limit $(\delta_1,\delta_2)\to(0,0)$ in order to form a distinguished path (or distinguished limit). This path must be parameterized with a single small parameter $\varepsilon\to 0$. Hence, the list must be replaced with $$\begin{aligned}
{L},\quad R,\quad{G}, \quad \varepsilon\label{params}\end{aligned}$$ A very general *ansatz* of asymptotic analysis is: all independent parameters and functions are allowed to be scaled with the use of a basic dimensionless parameter, which, in our case, is $\varepsilon$. It means, that we can introduce, say, dimensionless velocity ${\boldsymbol{u}}$ by an expression ${\boldsymbol{u}}^*=\varepsilon^k{L}^{1/2}{G}^{1/2}{\boldsymbol{u}}$ with an arbitrary real number $k$ (not just ${\boldsymbol{u}}^*={L}^{1/2}{G}^{1/2}{\boldsymbol{u}}$), *etc*. This kind of scaling is common in the asymptotic analysisand in perturbation methods, and, in particular, in fluid dynamics, where, say, the length-scale for viscous flows can be chosen as $L\sqrt{Re}$ ($Re$ is Reynolds number). Hence, as a starting point of our asymptotic consideration, we introduce a ‘flexible’ scaling, depending on indefinite real parameters $q,k,m,n,r$ $$\begin{aligned}
&&{\boldsymbol{x}}^*=L{\boldsymbol{x}}\label{scale-x}\\
&&\rho^*=R(1+\varepsilon^q\sigma),\quad q>0\label{scale-rho}\\
&& {\boldsymbol{u}}^*=\varepsilon^k{L}^{1/2}{G}^{1/2}{\boldsymbol{u}}\label{scale-u}\\
&&  t^*= {\varepsilon}^{-k}{L}^{1/2}{G}^{-1/2}\ {t}\label{scale-time}\\
&& \Phi^*=\varepsilon^m{L}\, {G}\, \Phi\label{scale-Phi}\\
&& p^*=R{{L}{G}}(-\varepsilon^m\Phi+\varepsilon^n \pi)\label{scale-p}\\
&& \mu^*=\varepsilon^rR {L}^{3/2}{G}^{1/2}\nu\label{scale-mu}\end{aligned}$$ where $\nu={\mathrm{const}}$. The following *imposed requirements* (IR1-4) are crucial for the understanding of -:

(IR1) The key restriction $q>0$ provides the smallness of density deviations (from constant density $R$), required by $\delta_1\ll 1$ ; one can take $\delta_1=\varepsilon^q$.

(IR2) The introduction of dimensionless time $t$ with the same parameter $k$ as in ${\boldsymbol{u}}^*$ is required for the balance $$\begin{aligned}
{||{\boldsymbol{u}}^*_{t^*}||}\sim {||({\boldsymbol{u}}^*\cdot\nabla^*){\boldsymbol{u}}^*||} \label{nonlin}\end{aligned}$$ which represents a part of . The restriction can be revoked in the papers which concentrate solely on linearized models, see [@Wood]; we do not study this simplified direction here. Hence, the field of acceleration is $${\boldsymbol{a}}^*=\varepsilon^{2k} G{\boldsymbol{a}}({\boldsymbol{x}},t);\quad {\boldsymbol{a}}\equiv{\boldsymbol{u}}_t+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}$$

(IR3) The appearance of the same power $m$ in and is necessary for the mutual cancelation of a gravity term and a pressure term, describing a hydrostatic equilibrium in and .

(IR4) We set all dimensionless functions and parameters $\sigma, {\boldsymbol{u}}, \Phi, \pi, \nu$ in - of order one $$\begin{aligned}
\label{O1}
&&\sigma({\boldsymbol{x}},t)\sim O(1),\ {\boldsymbol{u}}({\boldsymbol{x}},t)\sim {\boldsymbol{O}}(1),\ \pi({\boldsymbol{x}},t)\sim O(1),\ \Phi({\boldsymbol{x}})\sim O(1),\ \nu\sim O(1)\end{aligned}$$ which means that dimensionless magnitudes of the related fields and parameters are described only by related powers of $\varepsilon$. We also accept that all required for our study derivatives of $\sigma, {\boldsymbol{u}}, \Phi, \pi$ are of order one.

The substitution of - into , yields $$\begin{aligned}
&&(1+\varepsilon^q \sigma)\varepsilon^{2k}[{\boldsymbol{u}}_t+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}]=-\varepsilon^n\nabla \pi-\varepsilon^{q+m}\sigma\nabla\Phi+\varepsilon^{k+r}\nu\Delta {\boldsymbol{u}}\label{pre-Buss-0}\\
&& \sigma_t+({\boldsymbol{u}}\cdot\nabla)\sigma=0,\quad \nabla\cdot{\boldsymbol{u}}=0\nonumber\end{aligned}$$ where the required (for the obtaining of EBA) setting is to choose all involved terms of the same order in $\varepsilon$: $$\begin{aligned}
2k=n=q+m=k+r\label{relation}\end{aligned}$$ It leads to the equations $$\begin{aligned}
&&(1+\varepsilon^q \sigma)[{\boldsymbol{u}}_t+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}]=-\nabla \pi-\sigma\nabla\Phi+\nu\Delta {\boldsymbol{u}}\label{pre-Buss-1}\\
&& \sigma_t+({\boldsymbol{u}}\cdot\nabla)\sigma=0,\quad \nabla\cdot{\boldsymbol{u}}=0\nonumber\end{aligned}$$ which below are referred to as the *main equations*, which explicitly contains a small parameter $\delta\equiv \varepsilon^q$. Since the relations give us three connections between five parameters, we rewrite - with the use of two parameters $q$ and $k$ only: $$\begin{aligned}
&&{\boldsymbol{x}}^*=L{\boldsymbol{x}}\label{scale-x-1}\\
&&\rho^*=R(1+\varepsilon^q\sigma),\quad q>0\label{scale-rho-1}\\
&& {\boldsymbol{u}}^*=\varepsilon^k{L}^{1/2}{G}^{1/2}{\boldsymbol{u}}\label{scale-u-1}\\
&&  t^*= {\varepsilon}^{-k}{L}^{1/2}{G}^{-1/2}\ {t}\label{scale-time-1}\\
&& \Phi^*=\varepsilon^{2k-q}{L}\, {G}\, \Phi\label{scale-Phi-1}\\
&& p^*=R{{L}{G}}(-\varepsilon^{2k-q}\Phi+\varepsilon^{2k} \pi)\label{scale-p-1}\\
&& \mu^*=\varepsilon^kR {L}^{3/2}{G}^{1/2}\nu\label{scale-mu-1}\end{aligned}$$ The *observations* O1-3 are:

(O1) Equations , taken for an arbitrary value of $\varepsilon$, still give us the exact governing equations describing incompressible stratified flows. However, for $\varepsilon\to 0$ these equations describe only some special classes of stratified flows. Taken together, equations - give us an infinite number of asymptotic models of incompressible stratified flows, one model for each pair $(q,k)$.

(O2) We notice, that the dimensionless equations of motion contain only parameter $q$, while the underlining scaling - possesses an additional freedom given by two parameters $q$ and $k$. As soon as $q$ is chosen, the restriction on the choice of $k$ is in related physics: we have to decide what is the asymptotic behaviour of $\Phi^*\sim\varepsilon^{2k-q}$ at $\varepsilon\to 0$. It is ‘singular’ for ${2k-q}<0$ however it does not prevent us against of using it (see Discussion).

(O3) For any $k$, the amplitudes of velocity and viscosity are of the order $\varepsilon^k$. It brings us to an expected conclusion: the Reynolds numbers ($Re$) are of order one in all viscous flows. Indeed, in we have accepted that a viscous term is of the same order as the other terms in the Navier-Stokes equations, which corresponds to $Re\sim 1$. The consideration of high or low $Re$ will take us to different asymptotic models, which we avoid in this paper.

Three levels of Boussinesq approximations
=========================================

Due to $q>0$, the main equation explicitly contains a small parameter; hence its solutions for $\varepsilon\to 0$ can be considered as series in $\varepsilon$: $$\begin{aligned}
\label{series1}
(\sigma,{\boldsymbol{u}},\pi)=\sum_{\alpha=0}^\infty \varepsilon^\alpha (\sigma_\alpha,{\boldsymbol{u}}_\alpha,\pi_\alpha),\quad \alpha=0,1,2,\dots\end{aligned}$$ The substitution of into produces the equations of successive approximations in $\varepsilon$. The use of integers $\alpha$ means that we should operate only with the integers $q=1,2,\dots$ and $k=0, \pm1,\pm2,\dots$. The leading approximation ($\alpha=0$) is $$\begin{aligned}
&&{\boldsymbol{u}}_{0t}+({\boldsymbol{u}}_0\cdot\nabla){\boldsymbol{u}}_0=-\nabla \pi_0-\rho_0\nabla\Phi+\nu\Delta {\boldsymbol{u}}_0\label{pre-Buss-1-0}\\
&& \rho_{0t}+({\boldsymbol{u}}_0\cdot\nabla)\rho_0=0,\quad \nabla\cdot{\boldsymbol{u}}_0=0\nonumber\end{aligned}$$ which can be immediately recognized as *the equations of Boussinesq approximation* (EBA) for unknown fields $\rho_0,{\boldsymbol{u}}_0,\pi_0$. To emphasise the use of EBA as the ‘governing’ equations of stratified flows, we change the notations in as $\sigma_0\mapsto \rho,\ {\boldsymbol{u}}_0\mapsto{\boldsymbol{v}},\ \pi_0\mapsto p$. That produces a conventional dimensionless form of EBA: $$\begin{aligned}
&&{\boldsymbol{v}}_{t}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}=-\nabla p-\rho\nabla\Phi+\nu\Delta{\boldsymbol{v}}\label{EBA}\\
&& \rho_{t}+({\boldsymbol{v}}\cdot\nabla)\rho=0,\quad \nabla\cdot{\boldsymbol{v}}=0\nonumber\end{aligned}$$ Their solutions also can be presented as the amplitude series $$\begin{aligned}
\label{seriesB}
(\rho,{\boldsymbol{v}},p)=\sum_{\beta=0}^\infty \varepsilon^\beta (\rho_\beta,{\boldsymbol{v}}_\beta,p_\beta),\quad \beta=0,1,2,\dots\end{aligned}$$ where we have accepted, that the amplitudes are described by the same small parameter $\varepsilon$. Here we should explain, why the amplitude parameter for the main equations and EBA is the same $\varepsilon$. Indeed, it is possible to introduce different amplitude parameters, however such a complication can be seen as an introduction of two different small parameters for the same physical entity (for the amplitude of a solution).

The linearization of (which corresponds to the term $\beta=1$ in ) is $$\begin{aligned}
&&{\boldsymbol{v}}_{1t}+({\boldsymbol{v}}_0\cdot\nabla){\boldsymbol{v}}_1+({\boldsymbol{v}}_1\cdot\nabla){\boldsymbol{v}}_0=-\nabla p_1-\rho_1\nabla\Phi+\nu\Delta{\boldsymbol{v}}_1\label{EBA-lin}\\
&& \rho_{1t}+({\boldsymbol{v}}_0\cdot\nabla)\rho_1+({\boldsymbol{v}}_1\cdot\nabla)\rho_0=0,\quad \nabla\cdot{\boldsymbol{v}}_1=0\nonumber\end{aligned}$$ At the same time, the first approximation (which corresponds to the term $\alpha=1$ in ) of the main equations (written for the lowest value $q=1$), is: $$\begin{aligned}
&&{\boldsymbol{u}}_{1t}+({\boldsymbol{u}}_0\cdot\nabla){\boldsymbol{u}}_1+({\boldsymbol{u}}_1\cdot\nabla){\boldsymbol{u}}_0 +{\boldsymbol{A}}_0=-\nabla \pi_1-\sigma_1\nabla\Phi+{\boldsymbol{u}}\Delta{\boldsymbol{u}}_1\label{pre-Buss-1-1}\\
&& \sigma_{1t}+({\boldsymbol{u}}_0\cdot\nabla)\sigma_1+({\boldsymbol{u}}_1\cdot\nabla)\sigma_0=0,\quad \nabla\cdot{\boldsymbol{u}}_1=0\nonumber\end{aligned}$$ Comparison of and shows that $$\begin{aligned}
{\boldsymbol{A}}_0\equiv \rho_0 {\boldsymbol{a}}_0=\rho_0 [{\boldsymbol{u}}_{0t}+({\boldsymbol{u}}_{0}\cdot\nabla){\boldsymbol{u}}_{0}]\label{spoil}\end{aligned}$$ represents a ‘spoiling term’, making the linearized version of the main equations different from the linearized version of EBA . If we take $q= 2$ or above, then this ‘spoiling term’ does not appear in the equations corresponding to $\alpha=1$; for the first time it appears in the approximation of the order $\alpha=q$. Hence, we introduce *three different levels for the quality of approximations of the main equations given by EBA:* $$\begin{aligned}
&&\text{\underline{Level 1}}:\quad  q=1\quad \text{\texttt{poor approximation}}\label{levels}\\
&&\text{\underline{Level 2}}:\quad q=2\quad \nonumber \text{\texttt{reasonable approximation}}\\
&&\text{\underline{Level 3}}:\quad q\geq 3\quad \text{\texttt{good approximation}} \nonumber\end{aligned}$$ For a poor approximation (Level 1) only EBA itself coincides with the leading term $\alpha=0$ of the main equation, while the linearized versions $\alpha=1$ and $\beta=1$ are different from each other. For a reasonable approximation (Level 2) we have the coinciding of EBA and its linearization $\beta=1$ with the approximations $\alpha=0$ and $\alpha=1$ of the main equations. Physically, it guarantees the same results for linear waves and linear stability for the EBA and for the main equations. For a good approximation (Level 3) the series for EBA coincides with that for the main equations up to $\alpha=2$. Physically, this means that both linear perturbations and weakly nonlinear perturbations for EBA and for the main equations are the same.

Classification of the EBA models.
=================================

Asymptotic models with $k=0$, $k>0$, and $k<0$ are qualitatively different, hence we consider them separately. For equations - yield: $$\begin{aligned}
&&{\boldsymbol{x}}^*=L{\boldsymbol{x}}\label{scale-x-10}\\
&&\rho^*=R(1+\varepsilon^q\sigma),\quad q>0\label{scale-rho-10}\\
&& {\boldsymbol{u}}^*={L}^{1/2}{G}^{1/2}{\boldsymbol{u}}\label{scale-u-10}\\
&&  t^*= {L}^{1/2}{G}^{-1/2}\ {t}\label{scale-time-10}\\
&& \Phi^*=\varepsilon^{-q}{L}\, {G}\, \Phi\label{scale-Phi-10}\\
&& p^*=R{{L}{G}}(-\varepsilon^{-q}\Phi+ \pi)\label{scale-p-10}\\
&& \mu^*=R {L}^{3/2}{G}^{1/2}\nu\label{scale-mu-10}\end{aligned}$$ One can see, that an advantage of this model is the finite limits (at $\varepsilon\to 0$) for ${\boldsymbol{u}}^*$, $t^*$, $\mu^*$, and for the deviation $RLG\pi$ of pressure from a hydrostatic reference state. At the same time, the limit for $\Phi^*$ becomes more ‘singular’ with the increase of $q$ from Level 1 to Level 3. In other words, these models correspond to ‘order one’ velocity and viscosity and to ‘usual physical time’ in the case of strong gravity. It is important to notice, that a ‘singular’ limit for gravity is perfectly permitted and actively used, see Discussion.

The are characterised by ‘slow time’ $t^*$, small velocity ${\boldsymbol{u}}^*$ and small viscosity $\nu^*$, which appear in - for $k>0$. The gravity is weak for $2k-q>0$ and ‘strong’ for $2k-q<0$.

, giving a reasonable approximation , corresponds to $k=1$ and $q=2$ $$\begin{aligned}
&&{\boldsymbol{x}}^*=L{\boldsymbol{x}}\label{scale-x-main}\\
&&\rho^*=R(1+\varepsilon^2\sigma)\label{scale-rho-main}\\
&& {\boldsymbol{u}}^*=\varepsilon\sqrt{{L}{G}}{\boldsymbol{u}}\label{uscale}\\
&&  t^*= \frac{1}{\varepsilon}\sqrt{{L}/{G}}\ {t}\label{scale-time-main}\\
&& \Phi^*={L}\, {G}\, \Phi\label{scale-Phi-main}\\
&& p^*=R{{L}{G}}(-\Phi+\varepsilon^{2} \pi)\label{scale-p-main}\\
&& \mu^*=\varepsilon R {L}^{3/2}{G}^{1/2}\mu\label{scale-mu-main}\end{aligned}$$ when the main equations are: $$\begin{aligned}
&&(1+\varepsilon^2 \sigma)[{\boldsymbol{u}}_t+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}]=-\nabla \pi-\sigma\nabla\Phi+\mu\Delta {\boldsymbol{u}}\label{pre-Buss-main}\\
&& \sigma_t+({\boldsymbol{u}}\cdot\nabla)\sigma=0,\quad \nabla\cdot{\boldsymbol{u}}=0\nonumber\end{aligned}$$ This case is characterised by: (i) a small amplitude of velocity $|{\boldsymbol{u}}^*|\sim\varepsilon$; (ii) small (of order $\varepsilon^2$) deviations of density and pressure from hydrostatic (reference) density and pressure; (iii) a finite limit for gravity $\Phi^*\sim\varepsilon^0$; (iv) a long time-scale $t^*\sim 1/\varepsilon$ (or a low characteristic frequency $\sim \varepsilon$); (v) small viscosity $\mu^*\sim \varepsilon$.

One may say that properties (i) and (iii) can be seen as physical advantages. Indeed, the smallness of velocity ${\boldsymbol{u}}^*\sim\varepsilon$ is natural for the waves of small amplitude or for generic small perturbations; such a smallness is accepted in the majority of related publications. The gravity $\Phi^*\sim\varepsilon^0$ also may be considered as ‘physically natural’. Some authors consider only this scaling, which is always written differently due to the presence of additional small parameters, see *e.g.* [@Long]). These additional parameters are related to, say, the long-wave approximation or small compressibility.

Finally, are characterised by the ‘fast time’ variable and the large amplitude velocity and viscosity, see - for $k<0$. The gravity is always strong. This class of models has been never studied or used, however the possible application could be to rapidly oscillating flows or to turbulence.

How to choose parameters $q$ and $k$?
=====================================

As one have seen the main equations contain only parameter $q=1,2,3,\dots$, while the underlining scaling - is defined by two parameters $q$ and $k=\pm 1,\pm 2,\pm 3,\dots$. In particular problems, these parameters can only appear from the prescribed boundary conditions and initial conditions.

since in we do not re-scale the spatial variable ${\boldsymbol{x}}$, the dimensionless initial conditions can be always written as $$\begin{aligned}
\rho({\boldsymbol{x}},0)=\rho^{\dag}({\boldsymbol{x}}),\quad {\boldsymbol{u}}({\boldsymbol{x}},0)={\boldsymbol{u}}^{\dag}({\boldsymbol{x}})\label{IC}\end{aligned}$$ with given dimensionless functions $\rho^{\dag}({\boldsymbol{x}})$ and ${\boldsymbol{u}}^{\dag}({\boldsymbol{x}})$. Hence, here we must accept that the magnitudes of the solution , coincide with that of the initial conditions .

is slightly different for time-independent and time-dependent cases. The former case requires the same orders of functions , , while the latter case additionally requires the same orders in the time-scale .

Summarising, we write that the introduced multiplicity of scalings allows us to consider any amplitudes of the functions and variables, as they appear in the initial and boundary conditions.

An interesting question arises from the particular appearance of $\varepsilon$ in the main equations , which contain only a small parameter $\delta\equiv\varepsilon^q$. It suggests that solutions can be represented as power series of either $\delta$ or $\varepsilon$. In the former case and should be replaced with $$\begin{aligned}
&&(1+\delta\, \rho)[{\boldsymbol{u}}_t+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}]=-\nabla p-\rho\nabla\Phi\label{pre-BussD}\\
&& \rho_t+({\boldsymbol{u}}\cdot\nabla)\rho=0,\quad \nabla\cdot{\boldsymbol{u}}=0\nonumber\end{aligned}$$-5mm $$\begin{aligned}
\label{seriesD}
(\rho,{\boldsymbol{u}},p)=\sum_{\gamma=0}^\infty \delta^\gamma (\rho_\gamma,{\boldsymbol{u}}_\gamma,p_\gamma),\quad \gamma=0,1,2,\dots\end{aligned}$$ which corresponds to the poor approximation , with $q=1$. However, one can also use parameter $\varepsilon$, take $q=2$ or $q\geq 3$ and obtain the reasonable or good approximations. The cost of this improvement will be the increasing of a small parameter. Indeed, if in (for a poor approximation) we choose, say $\delta=10^{-3}$, then for a good approximation with $\delta=\varepsilon^3$ we obtain $\varepsilon=10^{-1}$. Of course, the increasing of a small parameter makes the converging of series worse. Hence, *we have a counterbalance*: the improving of the model quality makes the convergence worse and for a lower model quality we improve the convergence of series. The optimization of this counterbalance represents an interesting practical question. At the same time we notice here that the conversion from $\varepsilon$ to $\delta$ is impossible if $\varepsilon$ appears explicitly in boundary conditions or in initial conditions, through -.

Another restriction on the values of $q$ and $k$ may appear from physics: for a particular application one may ‘prefer’ strong gravity or weak gravity. For example, if the value of $k$ (enforced, say, by boundary conditions) is $k=0$, then immediately gives us $\Phi^*\sim \varepsilon^{-q}$, hence, with the increase of the quality of the EBA model we have to consider the increasingly strong gravity field. On the other hand, if we wish to consider $\Phi^*\sim \varepsilon^{0}$ then we have to take $2k=q$, and the option to choose is $k=1, q=2$ - with slow time . Other options for $t^*$ will require the change of the magnitude of a gravity field.

The general conclusion of this section is: a flexible choice of $q$ and $k$ allows to consider a broad variety of problems with different orders of unknown functions and scales involved; simultaneously one can choose the asymptotic models of different (desired) quality.

Discussion
==========

(D1) It is surprising, that we were unable to find the scaling - and even the main equations in the literature on stratified flows.

(D2)A question could be raised: what is the relation between the internal time scales (say, the inverse Brunt-Vaisala frequency) and that of ? The answer is: after the internal time-scale is properly scaled, it coincides with .

(D3) Another question could be asked: why the main equations are ‘worse’ or ‘more difficult to solve’ than the EBA ? We do believe that the direct solving of the main equation is preferable for many studies.

(D4) The asymptotic for gravity $\Phi^*\sim\varepsilon^{2k-q}$ at $\varepsilon\to 0$ in is ‘singular’ for ${2k-q}<0$ however it does not prevent one against of using it. For example, the ‘singular’ asymptotic $\Phi^*\sim 1/\varepsilon$ has been actively used in the vibrational convection, see [@Zenkovskaya; @Simonenko; @GershuniV; @Levenshtam]. It is apparent that, without the use of this ‘singular’ asymptotic, the theory of vibrational convection can not be constructed. The physical meaning of the ‘singular’ asymptotic is: we deal with a strong (not infinite) gravity, since in applications $\varepsilon$ is always finite.

(D5) The separation of the EBA models into poor, reasonable, and good approximations requires further clarification. The particular impacts of ‘spoiling terms’ should be investigated separately. The target here could be the obtaining of an estimation which could be written as an upper bound for a norm $$\begin{aligned}
\label{close}
||{\boldsymbol{u}}-{\boldsymbol{v}},\rho-\sigma||\leq C(\varepsilon)||\sigma{\boldsymbol{a}}, {\boldsymbol{u}}-{\boldsymbol{v}},\rho-\sigma||_{t=0}\end{aligned}$$ where constant $C$ depends on $\varepsilon$ such that $C\to 0$ when $\varepsilon\to 0$. The presence of $\rho{\boldsymbol{a}}$ in represents an assumption that the difference between solutions is caused by the ‘spoiling term’, as well as by the differences in initial conditions. It is especially interesting to compile a list of examples, where particular solutions of EBA are qualitatively different from the related solutions of the main equations (taken for small $\varepsilon$). One such example is presented in [@Long].

(D6) Particular solutions, illustrating the results of this paper, must be essentially unsteady, since the steady versions of EBA and the exact equations for steady flows are mathematically equivalent to each other, see [@YihS]. The most useful for our studies would be an exact unsteady solution of EBA and the main equations, satisfying the same initial and boundary conditions. However, the unsteady exact solutions for stratified flows are almost unknown, a noticeable exception is [@CraikExact].

(D7) The presented derivations raise some interesting questions about the scalings used in the studies of nonlinear waves (including solitons) and nonlinear stability of shear flows in the presence of density stratification, such as considered in [@Craik; @Grimshaw1]. It would be interesting to incorporate the presented scaling required for EBA to the elaborated scalings of the related studies. The combining of two scalings, which have been introduced independently from each other, may bring some surprises.

(D8) Interesting developments can appear in various cases where the Boussinesq approximation is used simultaneously with additional independent small parameters (for example, long waves or an inverse frequency of externally imposed oscillations). It is likely that some developments (and even some interesting contradictions between different requirements to scalings used for different small parameters) may appear in the problem of vibrational convection and in the problem of stability and waves in a vibrating stratified fluid [@GershuniV; @Yudovich; @Zenkovskaya], where two different time scales are employed from very beginning.

(D9) The ‘filtering out’ of acoustic oscillations was performed in [@Ver; @Zey]. An alternative to this studies approach could be the deriving of averaged equations with the use of two-timing averaging over the high-frequency oscillations, as it has been done in [@VladIlin; @VladimirovL; @VladProc].

The authors would like to express special gratitude to Profs. A.D.D. Craik, R.H.J. Grimshaw, and H.K. Moffatt for reading this manuscript and making useful critical comments. Many thanks to Profs. I.A. Eltayeb, D.W. Hughes, A.R. Kacimov, M.R.E. Proctor, and M.M. Rahman for helpful discussions. This research is supported by the grant IG/SCI/ DOMS/16/13 from Sultan Qaboos University, Oman.

<span style="font-variant:small-caps;">Bowker, J.A., Hughes, D.W. and Kersale, E.</span> 2014 Incorporating velocity shear into the magneto-Boussinesq approximation. *Geophysical and Astrophysical Fluid Dynamics,* **108**, 553-567.

<span style="font-variant:small-caps;">Drazin, P.G. & Reid, W.H.</span> 1981 *Hydrodynamic Stability.* Cambridge, CUP.

<span style="font-variant:small-caps;">Craik, A.D.D.</span> 1985 *Wave Interactions and Fluid Flows.* Cambridge, CUP.

<span style="font-variant:small-caps;">Craik, A.D.D.</span> 1989 The stability of unbounded two- and three-dimensional flows subject to body forces: some exact solutions. *J. Fluid Mech.*, **198**, 275-292.

<span style="font-variant:small-caps;">Gershuni, G.Z. & Zhukhovitsky, E.M.</span> 1972 *Convective Stability of Incompressible Fluid.* Moscow: Nauka (in Russian). Translated from Russian by D. Louvish. Keter Publications, Jerusalem/Wiley, 1976.

<span style="font-variant:small-caps;">Gershuni, G.Z. & Lyubimov, D.V.</span> 1997 *Thermal Vibrational Convection.* NY: John Wiley & Sons.

<span style="font-variant:small-caps;">Grimshaw, R.H.J.</span> 1984 Wave action and wave-mean flow interactions, with applications to stratified shear flow. *Ann. Rev. Fluid Mech.*, **16**, 11-44.

<span style="font-variant:small-caps;">Grimshaw, R.H.J.</span> 2001 Internal solitary waves. *Proceedings of the International Conference dedicated to the 100th anniversary of A.A. Andronov; “Progress in Nonlinear Science”,* Volume II, Frontiers of Nonlinear Physics, ed. A.G. Litvak, Inst. Applied Physics, Nizhny Novgorod, 139-153.

<span style="font-variant:small-caps;">Levenshtam, V.B.</span> 1996 Justification of averaging method for thermal vibrational convection. *Dokl. Akad. Nauk.*, **349**, 5, 621-623 (in Russian).

<span style="font-variant:small-caps;">Lighthill, M. J.</span> 1978 *Waves in Fluids.* Cambridge: CUP.

<span style="font-variant:small-caps;">Long, R.R.</span> 1965 On the Boussinesq approximation and its role in the theory of internal waves. *Tellus*, **XVII**, 1, 46-52.

<span style="font-variant:small-caps;">Moffatt, H.K.</span> 1978 *Magnetic Field Generation in Electrically Conducting Fluid.* Cambridge: CUP.

<span style="font-variant:small-caps;">Simonenko, I.B.</span> 1972 Justification of averaging method for convection problem in the field of rapidly oscillating forces and for other parabolic equations. *Math. Sbornik*, **87(129)**, 2, 236-253 (in Russian).

<span style="font-variant:small-caps;">Spiegel, E.A. & Veronis, G.</span> 1960 On the Boussinesq approximation for a compressible fluid. *Astrophysical Journal,* **131**, 5, 442-447.

<span style="font-variant:small-caps;">Vladimirov, V.A.</span> 2012 Magnetohydrodynamic drift equations: from Langmuir circulations to magnetohydrodynamic dynamo? [*J. Fluid Mech.*]{}, **698**, 51-61.

<span style="font-variant:small-caps;">Vladimirov, V.A. & Ilin, K.I.</span> 2013, An asymptotic model in acoustics: Acoustic-drift equation. *J.Acoust.Soc.Am.*, **134**, 5, 3419-3424.

<span style="font-variant:small-caps;">Vladimirov, V.A., Proctor, M.R.E. and Hughes, D.W</span> 2015, Vortex dynamics of oscillating flows. *Arnold Math J.*, **239**, 2, 113-126.

<span style="font-variant:small-caps;">Wood, T.S. & Bushby, P.J.</span> 2016 Oscillatory convection and limitations of the Boussinesq approximation. *arXiv:* 1607.06005v1, \[physics.flu-dyn\].

<span style="font-variant:small-caps;">Yih, C.S.</span> 1965 *Dynamics of nonhomogeneous fluid.* NY: Macmillian.

<span style="font-variant:small-caps;">Yih, C.S.</span> 1977 *Fluid Mechanics.* Ann Arbor: West River Press.

<span style="font-variant:small-caps;">Yudovich, V.I.</span> 2006 Vibrodynamics and vibrogeometry of mechanical systems with constrains. *Uspehi Mekhaniki,* **4**, 3, 26-158 (in Russian).

<span style="font-variant:small-caps;">Zenkovskaya, S.M. & Simonenko, I.B.</span> 1966 On the influence of high frequency vibrations on the convection initiation. *Fluid Dynamics,* **1**, 5, 35-37.

<span style="font-variant:small-caps;">Zeytonian, R.K.</span> 2003 Joseph Boussinesq and his approximation: a contemporary view. *C.R. Mecanique,* **331**, 575-586.
