3.1. Linearization

Let us consider that instantaneous variables ( ${\it\zeta},{\it\eta}$ and ${\it\omega}$ ) consist of small fluctuations (denoted $^{\prime }$ ) around a steady-state solution (denoted $^{-}$ ):

(3.1a−c ) $$\begin{eqnarray}{\it\zeta}=\bar{{\it\zeta}}+{\it\zeta}^{\prime },\quad {\it\eta}=\bar{{\it\eta}}+{\it\eta}^{\prime },\quad {\it\omega}=\bar{{\it\omega}}+{\it\omega}^{\prime }.\end{eqnarray}$$

Equations (2.4), (2.5) and (2.8) are then linearized at the first order, leading to

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\it\zeta}^{\prime }}{\text{d}{\it\tau}}=-2{\it\alpha}{\it\lambda}{\it\zeta}^{\prime }+{\it\varepsilon}{\it\omega}^{\prime }, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\it\eta}^{\prime }}{\text{d}{\it\tau}}=-2{\it\alpha}{\it\lambda}\frac{1+\bar{{\it\eta}}}{1-\bar{{\it\zeta}}}(\bar{{\it\eta}}{\it\zeta}^{\prime }+\bar{{\it\zeta}}{\it\eta}^{\prime }), & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\omega}^{\prime }=-\frac{C_{d}^{2}}{{\it\Theta}\,\bar{{\it\omega}}}\left[\bar{{\it\eta}}{\it\zeta}^{\prime }-(1-\bar{{\it\zeta}}){\it\eta}^{\prime }\right]. & \displaystyle\end{eqnarray}$$

By combining them, we obtain the following second-order differential equation for the fluctuating variable ${\it\zeta}^{\prime }$ :

(3.5) $$\begin{eqnarray}\frac{\text{d}^{2}{\it\zeta}^{\prime }}{\text{d}{\it\tau}^{2}}+A_{1}\frac{\text{d}{\it\zeta}^{\prime }}{\text{d}{\it\tau}}+A_{2}{\it\zeta}^{\prime }=0.\end{eqnarray}$$

It is important to note that the same differential equation applies to the variables ${\it\eta}^{\prime }$ and ${\it\omega}^{\prime }$ , so that, in the case of oscillating regimes, the fluctuating variables ( ${\it\zeta}^{\prime }$ , ${\it\eta}^{\prime }$ and ${\it\omega}^{\prime }$ ) will have the same period and the same rate of damping. However, they could be phase-shifted. Using the relationships (2.11) between variables in steady state, coefficients $A_{1}$ and $A_{2}$ in (3.5) can be expressed as functions of $\bar{{\it\zeta}}$ and ${\it\Theta}$ alone:

(3.6a,b ) $$\begin{eqnarray}A_{1}={\it\alpha}{\it\lambda}\left(2+3\,\frac{\bar{{\it\zeta}}+{\it\Theta}}{1-\bar{{\it\zeta}}}\right)\quad \text{and}\quad A_{2}=2{\it\alpha}^{2}{\it\lambda}^{2}\,\frac{\bar{{\it\zeta}}+{\it\Theta}}{1-\bar{{\it\zeta}}}\left[3+\frac{\bar{{\it\zeta}}+{\it\Theta}}{1-\bar{{\it\zeta}}}+\frac{{\it\Theta}}{\bar{{\it\zeta}}}\right].\end{eqnarray}$$

3.2. Oscillating regimes

Equation (3.5) shows that underdamped oscillations are possible provided that ${\it\Delta}={A_{1}}^{2}-4\,A_{2}<0$ . By replacing $A_{1}$ and $A_{2}$ by their expressions (3.6), this condition can be written as

(3.7) $$\begin{eqnarray}{\it\Delta}={\it\alpha}^{2}{\it\lambda}^{2}\frac{17\bar{{\it\zeta}}^{3}+(22{\it\Theta}-20)\bar{{\it\zeta}}^{2}+(9{\it\Theta}^{2}-20{\it\Theta}+4)\bar{{\it\zeta}}-8{\it\Theta}^{2}}{\bar{{\it\zeta}}(1-\bar{{\it\zeta}})^{2}}<0,\end{eqnarray}$$

or alternatively

(3.8) $$\begin{eqnarray}{\it\Theta}^{2}(9\bar{{\it\zeta}}-8)+{\it\Theta}\bar{{\it\zeta}}(22\bar{{\it\zeta}}-20)+\bar{{\it\zeta}}(17\bar{{\it\zeta}}^{2}-20\bar{{\it\zeta}}+4)<0.\end{eqnarray}$$

This second formulation immediately indicates that, for high values of ${\it\Theta}$ , oscillations will not be possible if $\bar{{\it\zeta}}>8/9\approx 0.889$ .

When ${\it\Delta}<0$ , we can define the two following indicators: the dimensionless half-period ${\it\tau}_{1/2}=2{\rm\pi}/\sqrt{-{\it\Delta}}$ , and the associated damping ratio ${\it\beta}=|{\it\zeta}^{\prime }({\it\tau}+{\it\tau}_{1/2})/{\it\zeta}^{\prime }({\it\tau})|=\exp \left(-{\rm\pi}A_{1}/\sqrt{-{\it\Delta}}\right)$ , which is defined by the ratio of any two successive amplitudes on a half-period (i.e. between two successive overshoots). These indicators are expressed as

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{1/2}=\frac{2{\rm\pi}}{{\it\alpha}{\it\lambda}}\frac{1-\bar{{\it\zeta}}}{\sqrt{-17\bar{{\it\zeta}}^{2}-(22{\it\Theta}-20)\bar{{\it\zeta}}-(9{\it\Theta}^{2}-20{\it\Theta}+4)+8\displaystyle \frac{{\it\Theta}^{2}}{\bar{{\it\zeta}}}}}, & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\beta}=\exp \left[-{\rm\pi}\frac{\bar{{\it\zeta}}+3{\it\Theta}+2}{\sqrt{-17\bar{{\it\zeta}}^{2}-(22{\it\Theta}-20)\bar{{\it\zeta}}-(9{\it\Theta}^{2}-20{\it\Theta}+4)+8\displaystyle \frac{{\it\Theta}^{2}}{\bar{{\it\zeta}}}}}\right]. & \displaystyle\end{eqnarray}$$

As mentioned earlier, the time evolution of the fluctuating quantities ${\it\zeta}^{\prime }$ , ${\it\eta}^{\prime }$ and ${\it\omega}^{\prime }$ may be out of phase. Denoting ${\it\Phi}_{1}$ as the phase shift between ${\it\zeta}^{\prime }$ and ${\it\eta}^{\prime }$ , and ${\it\Phi}_{2}$ as the phase shift between ${\it\zeta}^{\prime }$ and ${\it\omega}^{\prime }$ , we obtain

(3.11a,b ) $$\begin{eqnarray}{\it\Phi}_{1}={\rm\pi}-\arctan \left[\frac{A_{2}}{{\it\alpha}{\it\lambda}\left(1+\displaystyle \frac{1}{2}\displaystyle \frac{\bar{{\it\zeta}}+{\it\Theta}}{1-\bar{{\it\zeta}}}\right)\sqrt{-{\it\Delta}}}\right]\quad \text{and}\quad {\it\Phi}_{2}={\rm\pi}-\arctan \left[\frac{A_{2}}{{\it\alpha}{\it\lambda}\sqrt{-{\it\Delta}}}\right],\end{eqnarray}$$

using the linearized equations (3.2)–(3.4). These phase shifts and the damping ratio ${\it\beta}$ (the other intrinsic indicator of the oscillatory system) are independent of ${\it\alpha}{\it\lambda}$ , i.e. of the box width $L$ (since $A_{2}\propto {\it\alpha}^{2}{\it\lambda}^{2}$ and $\sqrt{-{\it\Delta}}\propto {\it\alpha}{\it\lambda}$ ). By contrast, the quantities related to the physical time (period and time shifts) obviously do depend upon the value of $L$ . When the oscillating regime appears ( ${\it\Delta}\lesssim 0$ ), relations (3.11) indicate that both signals for layer density and outflow velocity are $90^{\circ }$ phase-shifted with respect to the interface height signal. For lower (negative) values of ${\it\Delta}$ , these shifts are expected to increase. In all cases, the phase shift between ${\it\zeta}^{\prime }$ and ${\it\omega}^{\prime }$ will be smaller than the one between ${\it\zeta}^{\prime }$ and ${\it\eta}^{\prime }$ , indicating that the outflow velocity should react before the layer density to an interface height fluctuation. This observation is consistent with the order in which the first overshoots appear for ${\it\zeta}$ , ${\it\omega}$ and ${\it\eta}$ in figure 3(a).

3.3. Case of ‘small plumes’

Before addressing the general non-Boussinesq case, it is interesting to first focus on ‘small plumes’ for which ${\it\Theta}\longrightarrow 0$ . Indeed, for ${\it\Theta}=0$ , the results so far are those that would have been obtained if we had tackled the problem under the Boussinesq approximation, following Kaye & Hunt (Reference Kaye and Hunt2004). In this case, from inequality (3.7), the existence of oscillations simply leads to

(3.12) $$\begin{eqnarray}17\bar{{\it\zeta}}^{2}-20\bar{{\it\zeta}}+4<0\quad \Longrightarrow \quad 0.255<\bar{{\it\zeta}}<0.921.\end{eqnarray}$$

The upper value $0.921$ was already found by Kaye & Hunt (Reference Kaye and Hunt2004) from a numerical search aimed at quantifying the critical condition in which a (downward) ‘overshoot’ is possible. These authors used a geometrical parameter ${\it\mu}=2{\it\alpha}H/e$ (whose critical value ${\it\mu}_{c}$ was $0.32$ ) related to our parameter ${\it\Lambda}$ ( ${\it\Lambda}=1/{\it\mu}^{2}$ ). They also reinterpreted this parameter as being the ratio of two competing time scales (a filling time scale and a draining time scale), which are given for a line plume by $T_{f}\sim L/B^{1/3}$ and $T_{d}\sim LH/eB^{1/3}$ , respectively. In this light, one can deduce that:

1. (i) for $\bar{{\it\zeta}}>0.921$ , draining dominates the process ( $T_{d}\ll T_{f}$ ), and the perturbation ${\it\zeta}^{\prime }$ varies as $+\!\exp (-a{\it\tau})$ (overdamped regime for ${\it\Delta}>0$ ) in the vicinity of the horizontal asymptote corresponding to the steady state;

2. (ii) for $\bar{{\it\zeta}}<0.255$ , filling dominates the process ( $T_{f}\ll T_{d}$ ), which creates an initial strong impulse followed by an overdamped regime. In this case, the perturbation ${\it\zeta}^{\prime }$ tends to the (steady-state) horizontal asymptote as $-\!\exp (-a{\it\tau})$ ;

3. (iii) for $0.255<\bar{{\it\zeta}}<0.921$ draining and filling compete, allowing oscillatory behaviour to appear. According to equation (2.9a,b ) (for ${\it\Theta}=0$ , ${\it\alpha}=0.1$ and $C_{d}=1/\sqrt{2}$ ), this corresponds to a top opening width $e$ ranging from 3 % to 63 % of the box height $H$ .

Additionally, from relations (3.9) and (3.10) with ${\it\Theta}$ set to $0$ , it can be shown that the oscillation frequency $1/{\it\tau}_{1/2}$ and the damping ratio ${\it\beta}$ will be maximal for $\bar{{\it\zeta}}=6/7$ and $\bar{{\it\zeta}}=6/11$ , respectively. The maximal value that ${\it\beta}$ can reach (at $\bar{{\it\zeta}}=6/11$ ) is $0.28\,\%$ . This leads to the conclusion that, for ‘small plumes’, although oscillations can exist, their damping will remain very high. The main results discussed here for ‘small plumes’ are summarized in figure 4.

Figure 4.

Main regimes for ‘small plumes’ as a function of the stabilized interface location (the corresponding ratio $e/H$ is calculated for ${\it\alpha}=0.1$ and for $C_{d}\approx 0.7$ ).

From relation (3.11), in which ${\it\Theta}$ is set to 0, figure 5 presents the evolution of the phase shifts ${\it\Phi}_{1}$ and ${\it\Phi}_{2}$ as a function of the stabilized interface location $\bar{{\it\zeta}}$ in the range where oscillations are possible ( $0.255<\bar{{\it\zeta}}<0.921$ ). As expected, ${\it\Phi}_{1}$ and ${\it\Phi}_{2}$ are minimal and close to $90^{\circ }$ at the limits of the oscillating regime. Whatever the value of $\bar{{\it\zeta}}$ , the phase shift between ${\it\zeta}^{\prime }$ and ${\it\eta}^{\prime }$ will always be greater than the one between ${\it\zeta}^{\prime }$ and ${\it\omega}^{\prime }$ . These conclusions on phase shifts are still valid in the general case (non-‘small plumes’) presented in the next subsection.

Figure 5.

Phase shifts ${\it\Phi}_{1}$ (between ${\it\zeta}^{\prime }$ and ${\it\eta}^{\prime }$ ) and ${\it\Phi}_{2}$ (between ${\it\zeta}^{\prime }$ and ${\it\omega}^{\prime }$ ) as a function of the stabilized interface location $\bar{{\it\zeta}}$ for ‘small plumes’.

Figure 6.

Main regimes in the general case (non-‘small plumes’) as a function of the dimensionless parameter ${\it\Theta}$ and the geometrical parameter ${\it\Lambda}$ .

3.4. General case – non-‘small plumes’

From relation (2.9a,b ) and inequality (3.8), conditions for which ${\it\Delta}<0$ are calculated in the general case ( ${\it\Theta}\neq 0$ ). In the ranges $10^{-3}<{\it\Theta}<1$ and $10^{-3}<{\it\Lambda}<100$ , the three regimes already mentioned for ‘small plumes’ are shown in figure 6. It can be seen that for ${\it\Theta}$ approximately greater than $0.1$ , the ‘overshoot’ overdamped regime may disappear, indicating that oscillations are possible provided that the upper layer thickness at steady state $h_{ss}$ is greater than $H/9$ (i.e. $\bar{{\it\zeta}}<8/9$ ). In order to quantify these oscillations, figure 7 compares the evolution of the damping ratio ${\it\beta}$ (in %) as a function of $\bar{{\it\zeta}}$ for several values of ${\it\Theta}$ including the envelope cases ${\it\Theta}=0$ and ${\it\Theta}\longrightarrow \infty$ (though an infinite value of ${\it\Theta}$ is unrealistic for the applications targeted). Note that these asymptotic cases allow the damping ratio to be expressed analytically as

(3.13a,b ) $$\begin{eqnarray}\displaystyle {\it\beta}_{({\it\Theta}=0)}=\exp \left[-{\rm\pi}\frac{\bar{{\it\zeta}}+2}{\sqrt{-17\bar{{\it\zeta}}^{2}+20\bar{{\it\zeta}}-4}}\right]\quad \text{and}\quad {\it\beta}_{({\it\Theta}\longrightarrow \infty )}=\exp \left[-{\rm\pi}\sqrt{\bar{{\it\zeta}}/(8/9-\bar{{\it\zeta}})}\right]\!. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Figure 7 shows that the damping ratio ${\it\beta}$ increases with ${\it\Theta}$ for a given value of $\bar{{\it\zeta}}$ . For high enough values of ${\it\Theta}$ (typically ${\it\Theta}\gtrsim 0.2$ ), ${\it\beta}$ continuously increases when $\bar{{\it\zeta}}$ decreases from an upper limit value (that is less than $0.921$ but greater than $0.889$ ) to $0$ . From this upper limit value to $\bar{{\it\zeta}}\approx 0.7$ , the damping ratio remains very weak and similar to the one obtained for ‘small plumes’. Beyond this value, it can ‘take off’ to reach more significant values as $\bar{{\it\zeta}}\longrightarrow 0$ . In this situation, even if the filling process normally dominates the draining process, the interface descends sufficiently to positions where the vertical density gradient in the (non-Boussinesq) plume is strong. The (draining) potential energy of the layer then undergoes marked variations that contribute to reinforcing the oscillations. This conclusion must be qualified, however, by the fact that values of $\bar{{\it\zeta}}$ close to $0$ are no longer realistic when we consider that, in the present model, the buoyant fluid emerges from a source of finite area. This issue is discussed in the next subsection.

Figure 7.

Damping ratio ${\it\beta}$ (in %) as a function of $\bar{{\it\zeta}}$ for several values of ${\it\Theta}$ (axes are reversed for clarity).

3.5. Sources of finite area

For a buoyant fluid emitted (upwards) from a source of finite area (round or 2D), the vertical evolution of the plume variables (velocity $u$ , radius or thickness $b$ and density ${\it\rho}$ ) can be exactly and explicitly calculated from a so-called ‘plume function’ ${\it\Gamma}(z)$ , which quantifies locally the ratio between buoyancy and momentum (see Hunt & Kaye (Reference Hunt and Kaye2005) and Michaux & Vauquelin (Reference Michaux and Vauquelin2008) for extension to non-Boussinesq plumes). Unfortunately, ${\it\Gamma}(z)$ is expressed in terms of an incomplete beta function, which limits its usefulness in solving analytically the emptying–filling box problem. Alternatively, and more simply, it is possible to modify self-similar solutions (2.1) by introducing a virtual origin correction ( $z$ is replaced by $z+z_{v}$ where $z_{v}$ is the virtual source location). Convenient analytical expressions for the virtual origin location have been offered by Hunt & Kaye (Reference Hunt and Kaye2001) and by Candelier & Vauquelin (Reference Candelier and Vauquelin2012) amongst others. From a mathematical point of view, such a shift is equivalent to considering a line (or point source) plume in a box of height $H+z_{v}$ , but with an upper layer thickness that cannot exceed $H$ . As a consequence, the domain of variation of $\bar{{\it\zeta}}$ is strictly limited between $1$ and $z_{v}/(H+z_{v})$ , instead of being between $1$ and $0$ as shown on the ordinate axis in figure 7. For initially buoyancy-dominated plumes (lazy plumes), the virtual origin is expected not to be very far from the actual source, whereas for initially momentum-dominated plumes (forced plumes), $z_{v}$ would represent a non-negligible fraction of the box height. This may limit the range of application of the results presented in this article.