Flow regimes in emptying–filling boxes with two buoyancy sources of differing strengths and elevations

Abstract

Emptying–filling boxes have been studied in a wide range of configurations for decades, but the flow created in the box by two plumes rising from sources of arbitrary strength and elevation was previously unsolved. Guided by experiments and simplified analytical modelling, we reveal a rich array of two- and three-layer stratifications across seven possible flow regimes. The governing equations for these regimes show how the prevailing regime and stratification properties vary with three key parameters: the relative strength of the plumes, the height difference between their sources and a parameter characterising the resistance of the box to emptying. We observe and explain new behaviours not described in previous studies that are crucial to understanding emptying–filling boxes with multiple plumes. In particular, we demonstrate that the oft-assumed premise that $n$ plumes leads to a stratification with $n+1$ layers is not necessarily true, even in the absence of mixing. Two emptying–filling box models are developed: an analytical model addressing all combinations of the governing parameters and an extended model for three-layer stratifications that incorporates two mixing processes observed in the experiments. The predictions of these two models are in generally excellent agreement with measurements from the experimental campaign covering 69 combinations of the governing parameters. This study improves our understanding of emptying–filling boxes and could facilitate improvements to natural ventilation building design, as demonstrated by an example scenario in which occupants feel cooler upon the addition of a second source of heat.

1. Introduction

The stratification established by turbulent miscible plumes within sealed containers of fluid (the so-called ‘filling box’ after Baines & Turner 1969) or containers in which the interior connects to an exterior environment via openings (the ‘emptying–filling box’ after Linden, Lane-Serff & Smeed 1990) has garnered considerable interest given the range of buoyancy-driven stratified flows in the manmade and natural environments that fall into these categories. Examples of filling boxes in the natural environment include underground chambers supplied by magma plumes (Campbell & Turner 1989) and ocean basins stratified by aqueous plumes (Wong, Griffiths & Hughes 2001). Our focus herein concerns the aptly named emptying–filling box, for which buoyant fluid is supplied to and simultaneously empties from the box – the latter via openings in the base and top that link the interior and exterior environments.

The wide range of possible variations in box geometry and in the manner in which buoyancy is supplied results in a myriad of different emptying–filling box configurations. We consider the flow patterns and stratifications established by two buoyancy sources of differing strengths where one source is elevated relative to the other. This configuration covers scenarios of practical interest including indoor spaces with multiple occupied levels, such as atria, mezzanines and tiered auditoria. In contrast to the large area distributed buoyancy sources considered by Livermore & Woods (2007), we address the flow patterns resulting from localised buoyancy sources that give rise to turbulent plumes. The predictions of our theoretical model (§ 2), backed up by complementary flow visualisations and density measurements (§ 5), reveal this to be a rich and complex problem in which a multiplicity of distinct flow regimes are possible.

Mathematical descriptions of emptying–filling box flows have been developed to predict the steady stratification established by turbulent plumes rising from single or multiple localised buoyancy sources at the same elevation with identical strengths (Linden et al. 1990) and for sources with different strengths (Cooper & Linden 1996; Linden & Cooper 1996, hereafter CL96 and LC96). The transient flows that lead to a steady state are of considerable interest due to the complex behaviours they exhibit, including interface and thermal overshooting phenomena and sidewall overturning motions (Baines & Turner 1969; Kaye & Hunt 2004,  2007; Bower et al. 2008; Shrinivas & Hunt 2014a). Herein the transients have been tracked during experiments solely to confirm that a unique steady flow is reached, irrespective of the order in which the plume sources are activated – a key result on which our complementary steady-state theoretical analysis relies. The transient processes by which plumes create layers that evolve towards steady thicknesses and homogeneous densities is detailed elsewhere (e.g. Baines & Turner 1969; Linden et al. 1990; Kaye & Hunt 2004). Our focus hereinafter is exclusively on steady flows, on sources that produce turbulent plumes and on boxes with openings to the exterior in the base and top.

In the absence of plume–plume or plume–sidewall interactions, CL96 postulate that two buoyancy sources of unequal strengths that are located at an identical elevation produce a three-layer stratification (comprising homogeneous upper, intermediate and lower layers each separated by horizontal interfaces). For the case of positive buoyancy sources producing rising plumes, the stronger plume supplies the upper layer, the weaker plume the intermediate layer, and the lower layer is comprised of fluid drawn in from the exterior through the base opening (figure 1a). In addition to this basic arrangement, they modelled two secondary processes: a downward vertical transport of fluid caused by the impingement of the weaker plume on the upper interface, and the change in plume entrainment due to the reduction of buoyancy flux in the stronger plume as it enters the intermediate layer. These processes raised the height of the upper interface relative to the predictions made if neglected (figure 1b). Flow visualisations, such as the shadowgraph (figure 1c) from our experimental campaign (§ 5), confirm the general appropriateness of their theoretical conceptualisation.

Figure 1. (a) With both buoyancy sources on the box base, the plume with greater buoyancy flux forms the upper layer and the weaker plume forms the intermediate layer. (b) The coupling between plume and stratification results in two secondary effects: vertical transport of fluid from the upper to the intermediate layer of volume flux $Q_{I}^*$, and a modification of the entrainment into the stronger plume within the intermediate layer. These secondary effects increase the height of the upper interface. (c) Flow visualisation showing two buoyant layers and the impingement of the weaker plume (left) on the upper interface (experiment 43, § 5).

Building on the analysis of CL96, LC96 considered the more complex case of $n \geq 2$ equal-elevation plume sources. It was demonstrated by LC96 that accounting for the aforementioned reduction of the plume buoyancy flux in the stratification is not necessarily important. Neglecting both of the secondary processes, their simplified model provided a robust prediction of the stratification created by $n$ sources and offered design guidance for naturally ventilated buildings with multiple localised floor-level heat sources. In § 2 we initially adopt a similar approach by neglecting these secondary processes so as to develop an analytical model for the more complicated general case where two sources are at distinct elevations.

The remainder of this paper is structured as follows. The analytical model is presented in three parts, commencing in § 2 with the theoretical framework that includes a description of the flow structure and the associated governing equations for six distinct flow regimes. In § 3 we consider the limiting behaviour in each regime to identify the regime boundaries and reveal the dependence on the relative elevations of the two sources, the ratio of the source strengths and the resistance of the box to emptying. Regime boundaries are not straightforward to identify a priori and the regime maps we develop aid in determining which regime prevails for a given set of parameters and, hence, which governing equations are appropriate. There are stark differences between regimes in stratification and flow behaviour, and some example predictions are shown in § 4. In § 5 we describe the experimental campaign conducted to measure the stratification and provide confirmation of the predicted regimes. The campaign broadly affirms the predictions of the analytical model, and the flow visualisations reveal the significance of plume-induced mixing. This mixing results in a seventh regime as well as two mechanisms for vertical transport of fluid across a density interface, one of which concerns previously unreported mixing near the top of the box. Incorporating phenomenological models for the observed mixing into a more sophisticated model (§ 6) yields qualitative and quantitative improvements in stratification predictions. A brief discussion section (§ 7) considers some practical implications of our findings, including the potential merits of introducing an elevated heat source in a room in order to provide a cooler environment for occupants below. Finally, our conclusions are drawn in § 8.

2. Theoretical framework

We consider the range of possible steady stratification regimes that are established in an emptying–filling box by the turbulent plumes that originate from two localised buoyancy sources of differing strength and elevation. The ‘base’ and ‘elevated’ sources are located inside the box at the vertical coordinates $z=0$ and $z=k(\geq 0)$ with positive buoyancy fluxes $B_1$ and $B_2$, respectively. Both sources have zero volume and momentum flux.

The horizontal cross-sectional area of the box is independent of height and sufficiently large that the plumes are far away from each other and the walls, and thus, their vertical development is unaffected by the lateral confinement or plume–plume interaction. Similarly, the flow through the base openings (of total area $a_b$) and the top openings (of total area $a_t$) is assumed to have a negligible effect on plume development. The vertical distance between the base source and the top is equal to the box height $H$. We assume the plumes establish a stratification in the box consisting of one or two horizontal buoyant layers. Given the plumes cannot stratify the region below the lower source, a simple vertical offset in source elevation transforms cases where both sources are elevated to a configuration where the lower source can be treated as if on the base of the box.

Buoyant fluid rises through openings at the top at a flow rate $Q_n$ and, simultaneously, an inflow of fluid from a quiescent exterior of constant and uniform density is drawn through openings in the base. We restrict our analysis to cases for which the inflowing fluid does not mix with the buoyant layer(s) and unidirectional outflow is maintained at the top opening(s). These restrictions are routinely assumed for emptying–filling box models and the resulting flow is referred to as a displacement flow – the incoming fluid displacing, rather than mixing with, the buoyant fluid in the interior. Hunt & Coffey (2010) developed the constraints on the opening areas, layer depths and layer buoyancies that ensure displacement flow and, broadly speaking, this flow is expected for $a_t/a_b\lesssim 3$ for circular openings.

Three possible displacement flow configurations with a three-layer stratification are sketched in figure 2. The layers are delineated by the first and second interfaces at heights $h_1$ and $h_2$, respectively. The buoyancy of the upper, intermediate and lower layers are denoted $g'_2$, $g'_1$ and $g'_e$, respectively. Given we define buoyancy relative to the exterior density, $g'_e=0$. The plume volume flux $Q_{ij}$ at a given interface is identified by two integer subscripts, the first subscript identifies the plume source and the second identifies the interface – so, for example, $Q_{21}$ is the volume flux of the elevated plume at the level of the first interface. Each distinct flow configuration is referred to as a regime, and the three-layer regimes A3, B3 and C3 are distinguished by whether the elevated plume source supplies the intermediate layer (regime A3, figure 2a) or the upper layer (regimes B3 and C3), and whether the elevated source is below the first interface (regime B3, figure 2b) or above it (regime C3, figure 2c). We also identify three two-layer stratifications (referred to as regimes A2, B2 and C2, figure 3).

Figure 2. The idealised three-layer flows for (a) regime A3, (b) regime B3 and (c) regime C3. The horizontal dashed lines represent steady interfaces delineating layers of uniform buoyancy.

Figure 3. The idealised two-layer flows for (a) regime A2, (b) regime B2 and (c) regime C2.

Whilst, in general, a three-layer stratification might be anticipated with two plume sources, there are three physical explanations for the formation of only two layers. First, the source conditions can be such that the average buoyancy in each plume is identical at the height of the lower interface, a scenario discussed further in § 3.1.1 and that may be regarded as a limiting case of regimes A3 and B3 rather than novel behaviour. Second, the plumes can mix the upper portion of the box with sufficient vigour such that only a single well-mixed buoyant layer forms; this mixing is incorporated into the extended model (§ 6), but not the analytical model discussed at present. Third, the opening areas can be so large and, thus, provide so little resistance to flow through them that only a single buoyant layer forms, as idealised in the sketches of the two-layer regimes in figure 3. Regimes A2, B2 and C2 are similar to their three-layer counterparts, with the exception that the plume that was supplying the upper layer with buoyancy $g'_2$ now rises to the top of the box, and fluid from this plume exits without forming this layer. In the two-layer regimes the total flow rate through the upper openings ($Q_n$) is the sum of the volume flux of the exiting plume and the flow rate $Q_\omega$.

This third mechanism is not necessarily obvious and, indeed, CL96 were seemingly only partially aware of it in their analysis of emptying–filling boxes with two plume sources. When accounting for mixing in the box, CL96 explained that the second interface was pushed upwards and, depending on the opening areas and source conditions, predicted that it would reach the top of the box. At this point, they claimed that one of the plumes would reach the top of the box and exit immediately, without forming the third layer. However, they did not comment that their model also predicted that the second interface would reach the top of the box if the openings were large enough without the need for any mixing: governing systems of equations in both CL96 and LC96 have non-physical solutions where the second interface is predicted to be above the height of the box. Two of the key theoretical developments of this paper are the wider recognition of this mechanism, which applies in any emptying–filling box with multiple buoyancy sources with sufficiently large openings, and a more detailed description and model for how one of the plumes can exit the box without forming a layer.

Having described the idealised stratifications, the analysis proceeds by considering conservation of volume and buoyancy. At steady state, $h_1$ and $h_2$ are constant, and the net volume flux into any layer must be zero. Applying continuity for the three-layer regimes:

(2.1a{–}c)\begin{equation} \underbrace{Q_n=Q_{11}+Q_{21}=Q_{12}}_{\textit{Regime A3}}, \quad \underbrace{Q_n=Q_{11}+Q_{21}=Q_{22}}_{\textit{Regime B3}}, \quad \underbrace{Q_n=Q_{11}=Q_{22}}_{\textit{Regime C3}}. \end{equation}

For the two-layer regimes,

(2.2a,b)$$\begin{gather} \underbrace{Q_n = Q_{11} + Q_{21} = Q_{12} + Q_\omega}_{\textit{Regime A2}}, \quad \underbrace{Q_n = Q_{11} + Q_{21} = Q_{22} + Q_\omega}_{\textit{Regime B2}}, \end{gather}$$

(2.2c)$$\begin{gather} \underbrace{Q_n = Q_{11} = Q_{22} + Q_\omega}_{\textit{Regime C2}}. \end{gather}$$

Conservation of buoyancy requires that the layers have uniform buoyancy and

(2.3a{–}c)\begin{equation} \underbrace{g'_2=\frac{B_1+B_2}{Q_n}}_{\textit{Three-layer regimes}}, \quad \underbrace{g'_1=\frac{B_2}{Q_{21}}}_{\textit{A regimes}}, \quad \underbrace{g'_1=\frac{B_1}{Q_{11}}}_{\textit{B & C regimes}}. \end{equation}

The buoyant layers drive a flow rate $Q_L$ that is determined by balancing the pressure drop associated with the flow entering and exiting the box with the difference in hydrostatic pressure between the openings (cf. Linden et al. 1990). On determining the pressure drop through the openings using Bernoulli's equation, it is readily shown that

(2.4)\begin{equation} Q_L=A^*\left(\int_0^H g'(z)\,\mathrm{d} z\right)^{1/2} \ \textrm{with `effective area'}\ A^* = \left( \frac{1}{2c_t^2 a_t^2} + \frac{1}{2c_b^2 a_b^2} \right)^{{-}1/2}, \end{equation}

where $c_t$ and $c_b$ are coefficients associated with the top and base opening areas, respectively, incorporating the effects of the vena contracta and any frictional losses.

In the three-layer regimes the flow through the openings results solely from the hydrostatic pressure force from the buoyant layers, and evaluating the integral in (2.4) gives

(2.5)\begin{equation} Q_n=Q_L=A^*\sqrt{g'_2(H-h_2)+g'_1(h_2-h_1)}. \end{equation}

In regimes A2, B2 and C2, one plume exits the box directly and, thus, the layer-driven component accounts for only the portion $Q_\omega$ of the flow through the top opening(s) (figure 3). It is not obvious how the combination of these two flows should be best modelled, so we propose that the layer-driven flow occupies a proportion $\varOmega$ $(0<\varOmega <1)$ of the effective area. Accordingly, in the two-layer regimes,

(2.6)\begin{equation} Q_\omega=Q_L=\varOmega A^*\sqrt{g'_1(H-h_1)}. \end{equation}

Furthermore, we propose that $\varOmega =Q_L/Q_n$ such that the layer-driven flow uses a fraction of $A^*$ in proportion to its contribution to the total flow rate $Q_n$ (this proposal is consistent with the layer-driven flow using the areas $a_b$ and $a_t$ in this same proportion and the coefficients $c_t$ and $c_b$ remaining unmodified). This proposal implies that

(2.7)\begin{equation} Q_n=A^*\sqrt{g'_1(H-h_1)}, \end{equation}

a relationship which is identical to the form for the displacement flow driven by the two-layer stratification resulting from a single plume (cf. Linden et al. 1990), and may have been expected from the outset. Indeed, CL96 use this form without discussion of how the exiting plume affects the layer-driven flow. It is also a particularly appealing form for modelling purposes as the resulting solutions (§ 4) for quantities such as $h_1$ and $Q_n$ are continuous at the boundaries between the two- and three-layer regimes. The associated assumption that $c_t$ and $c_b$ are unchanged ensures that $A^*$ is only a function of the box geometry rather than the flow conditions. We note that a similar assumption is made in many emptying–filling box models, regardless of the presence of an exiting plume, where the known variation of $c_t$ with layer thickness and buoyancy (Hunt & Holford 2000) is neglected for simplicity.

In reality, we expect that $c_t$ and $\varOmega$ would vary depending on the flow conditions at the opening(s) at the top of the box, and, in particular, the position of the exiting plume in relation to the opening(s). One could imagine differences in flow behaviours for two limiting geometries: an exiting plume directly beneath a large opening such that the plume simply flows through it or, conversely, a plume impinging on the top of the box and forming a buoyant outflow current that flows to the opening(s) before exiting. However, we do not speculate herein on the significance of possible differences in these flow behaviours and, to avoid including the horizontal positions of the sources and openings in the model, treat these flows identically. Our complementary experiments (§ 5) were conducted in a box where the plume sources and openings were not vertically aligned and these confirmed that (i) a two-layer stratification can form even if the exiting plume impinges on the top of the box and (ii) that predictions for $h_1$ in experiments with a two-layer stratification based on (2.7) are very good. Accordingly, it seems that the simple proposed model for the layer-driven flow in the presence of the exiting plume for a two-layer stratification is not unreasonable.

In principle, it would be possible to have an opening so large that both plumes exit without forming a layer, yielding the trivial solution of one unstratified layer. This scenario is not explicitly modelled herein – indeed, it does not include the filling aspect of the emptying–filling box – although we note that two-layer solutions can approximate this scenario with an arbitrarily thin buoyant layer above a region at ambient density comprising nearly all of the box height.

For all regimes in the simplified model, we follow LC96 in neglecting the effects of the stratification on the plume volume flux and, thus,

(2.8a{–}c)\begin{equation} Q_{11}=C_QB_1^{1/3}h_1^{5/3}, \quad Q_{12}=C_QB_1^{1/3}h_2^{5/3}, \quad Q_{22}=C_QB_2^{1/3}(h_2-k)^{5/3}, \end{equation}

where $C_Q=(6\alpha _{p}/5)(9\alpha _{p}/10)^{1/3}{\rm \pi} ^{2/3}$ is the coefficient for the normalised plume volume flux and $\alpha _{p}$ is the top-hat entrainment coefficient of a plume (cf. Morton, Taylor & Turner 1956). Based on the van Reeuwijk & Craske (2015) summary of experimental and direct numerical simulation data, we take $\alpha _{p}=0.13$. When modelling regimes A2, B2 and C2, $h_2=H$. If the elevated plume source is above the first interface, $Q_{21}=0$, and

(2.9)\begin{equation} Q_{21}=C_QB_2^{1/3}(h_1-k)^{5/3} \quad \textrm{for } h_1\geq k. \end{equation}

Defining the non-dimensional quantities

(2.10a{–}e)\begin{align} \zeta_1=h_1/H, \quad \zeta_2=h_2/H, \quad \phi=k/H, \quad \psi=B_2/B_1 \quad \text{and} \quad R=C_Q^{3/2}H^2/A^*, \end{align}

the following relationship for the behaviour in regime A2 is obtained on combining (2.3), (2.1a{–}c), (2.7), (2.8) and (2.9):

(2.11)\begin{equation} \textrm{Regime A2:} \quad R=\sqrt{\frac{\psi(1-\zeta_1)}{\psi^{1/3}(\zeta_1-\phi)^{5/3}(\zeta_1^{5/3}+\psi^{1/3}(\zeta_1-\phi)^{5/3})^2}}. \end{equation}

The quantity $R$ represents the relative resistance of the box to the plume-driven flow through the openings: increasing the plume volume flux (by increasing $C_Q$ or $H$) or reducing $A^*$ results in higher velocities and, hence, pressure losses, at the openings. Alternatively, $R$ can be regarded as the ratio of a characteristic plume volume flux $Q_p=C_QB_1^{1/3}H^{5/3}$ to a characteristic layer-driven volume flux through the openings $Q_O$ where the typical layer buoyancy is taken as $B_1/Q_p$:

(2.12)\begin{equation} R=\frac{Q_p}{Q_O}=\frac{Q_p}{A^*\sqrt{H B_1/Q_p}}=\frac{C_Q^{3/2}H^2}{A^*}. \end{equation}

As continuity requires that the actual plume and layer-driven volume fluxes be equal in a given stratification, then larger values of $R$ imply that deeper layers are required to sufficiently reduce the sum of the plume volume fluxes at the interfaces below the characteristic value. Conversely, smaller values of $R$ imply shallower buoyant layers (to reduce the layer-driven volume flux below its characteristic value), and, for sufficiently small values of $R$, only one buoyant layer is required. Combining (2.1), (2.3), (2.5), (2.8) and (2.9) leads to the following pair of governing equations for regime A3:

(2.13)\begin{align} \textrm{Regime A3:} \quad & \zeta_2^{5/3}=\zeta_1^{5/3}+\psi^{1/3}(\zeta_1-\phi)^{5/3}, \end{align}

(2.14)\begin{align} & R=\sqrt{\frac{(1+\psi)(1-\zeta_2)+\psi^{2/3}\left(\dfrac{\zeta_2}{\zeta_1-\phi}\right)^{5/3} (\zeta_2-\zeta_1)}{\zeta_2^5}}. \end{align}

The governing equations for regime B, obtained by manipulation of (2.1)-(2.9), are

(2.15)\begin{align} \textrm{Regime B2:} \quad & R=\sqrt{\frac{1-\zeta_1}{\zeta_1^{5/3}(\zeta_1^{5/3}+\psi^{1/3}(\zeta_1-\phi)^{5/3})^2}}, \end{align}

(2.16)\begin{align} \textrm{Regime B3:} \quad & \psi^{1/3}(\zeta_2-\phi)^{5/3}=\zeta_1^{5/3}+\psi^{1/3}(\zeta_1-\phi)^{5/3}, \end{align}

(2.17)\begin{align} & R=\sqrt{\frac{(1+\psi)(1-\zeta_2)+\psi^{1/3}\left(\dfrac{\zeta_2-\phi}{\zeta_1}\right)^{5/3} (\zeta_2-\zeta_1)}{\psi(\zeta_2-\phi)^5}}. \end{align}

Similarly, for regime C,

(2.18)\begin{align} \textrm{Regime C2:} \quad & R=\sqrt{\frac{1-\zeta_1}{\zeta_1^5}}, \end{align}

(2.19)\begin{align} \textrm{Regime C3:} \quad &\zeta_1^{5/3}=\psi^{1/3}(\zeta_2-\phi)^{5/3}, \end{align}

(2.20)\begin{align} & R=\sqrt{\frac{(1+\psi)(1-\zeta_2)+(\zeta_2-\zeta_1)}{\zeta_1^5}}. \end{align}

We note that the governing equation for regime C2 is identical to the equation for an emptying–filling box with a single plume source on the base; thus, changing the properties of the elevated source has no effect on the stratification (provided the changes do not result in a switch of regimes).

3. Regime boundaries

By considering the limiting behaviour for each regime, we identify which regime occurs for a given combination of $\psi$, $\phi$ and $R$. The appropriate governing equation(s) from (2.11)–(2.20) can then be solved for $\zeta _1$ and, if applicable, $\zeta _2$. The plume volume fluxes at each interface can then be determined via (2.8)–(2.9) and these values used to calculate $g'_1$ and $g'_2$ via (2.3), and $Q_n$ via (2.1) or (2.2), as appropriate.

3.1. Boundaries of the A regimes

3.1.1. Boundary between the A and B regimes

The crossover between the A and B regimes occurs when the buoyancy of the base plume is equal to that of the elevated plume at the first interface, i.e. when $B_1/Q_{11}=B_2/Q_{21}$. Substitution of (2.8)–(2.9) into this condition shows that it is met if

(3.1)\begin{equation} \zeta_1=\frac{\phi}{1-\psi^{2/5}} \quad \textrm{for} \ 0<\phi<1 \quad \textrm{or} \quad \psi=1 \ \textrm{for} \ \phi=0. \end{equation}

If both buoyancy sources are on the base, then the relative strength of the plume sources alone determines which regime prevails, with $\psi =B_2/B_1 > 1$ resulting in a B regime. For ${\phi \neq 0}$, substitution of (3.1) into either the governing equations for regime A2 (2.11) or regime A3 (2.13)–(2.14) shows that the boundary between the A and B regimes is defined by

(3.2)\begin{equation} R_{AB}=\frac{(1-\psi^{2/5})^{2}}{1+\psi}\sqrt{\frac{1-\psi^{2/5}-\phi}{\phi^5}} \quad \textrm{for} \ 0<\phi<1, \end{equation}

where $R_{AB}(\psi,\phi )$ is the critical value that delineates the A and B regimes. If ${\psi ^{2/5}+\phi > 1}$ then $R_{AB}$ has an imaginary component and an A regime cannot occur; the buoyancy in the elevated plume will exceed the buoyancy in the base plume for any interface height $\phi <\zeta _1 < 1$. If $R < R_{AB}$ then an A regime will exist while if $R > R_{AB}$ then a B or C regime will occur. If $R=R_{AB}$ (or $\phi =0$ and $\psi =1$) then both plumes will reach $\zeta _1$ with the same average buoyancy. Consequently, no further stratification above $\zeta _1$ is possible, and thus, the boundary between the A and B regimes is characterised by a two-layer stratification.

3.1.2. Boundary between regimes A2 and A3

The boundary between regimes A2 and A3 can be found by substituting the condition $\zeta _2=1$ into (2.13)–(2.14). After simplification, this shows that

(3.3)\begin{equation} R_{A23}=\sqrt{\frac{\psi(1-\zeta_1)}{1-\zeta_1^{5/3}}} \quad \textrm{for}\ 1=\zeta_1^{5/3}+\psi^{1/3}(\zeta_1-\phi)^{5/3}, \end{equation}

where $R_{A23}(\psi,\phi )$ is the critical value delineating regime A2 ($R < R_{A23}$) from regime A3 ($R > R_{A23}$). Should a given pair of $\psi$ and $\phi$ require $\zeta _1\leq \phi$ to satisfy (3.3b) then neither regime A2 nor A3 can occur irrespective of the value of $R$ (as can also be shown by considering the boundary between the A and B regimes (§ 3.1.1).

Figure 4 shows regime maps for three example values of $\psi$. Focusing on the boundaries of the A regimes (which can only occur if $\psi < 1$), the map in figure 4(a) shows that the boundary between regimes A2 and A3 is nearly vertical (almost independent of $\phi$) and that, for most values of $\phi$, increasing $R$ results in regime A2 transitioning to regime A3, which transitions to regime B2 or B3 with further increase. However, the map also shows that there are some values of $\phi$ where $R_{A23} > R_{AB}$ and the flow will skip regime A3 and transition directly from regime A2 to B2 with increasing $R$ (e.g. $0.58 \lesssim \phi \lesssim 0.60$ for $\psi =1/10$). These qualitative boundary features are typical for regimes A2 and A3 across a range of $\psi$.

Figure 4. Regime maps in $R-\phi$ space for (a) $\psi =1/10$, (b) $\psi =1$ and (c) $\psi =10$. Regimes A (purple), B (green) and C (blue) are delineated by dashed lines. Dotted lines delineate two-layer regimes (lighter shading) from three-layer regimes (darker shading).

3.2. Boundaries between and within the B and C regimes

The boundary between the B and C regimes is set by substituting $\zeta _1=\phi$ into the governing equations and the boundaries within are set by substituting $\zeta _2=1$ into the three-layer governing equations. Applying both of these conditions simultaneously shows that all four boundaries intersect at a single critical point (figure 4). For a given value of $\psi$,

(3.4a,b)\begin{equation} R_{\textit{crit}}=\sqrt{\frac{(1+\psi^{1/5})^4}{\psi}} \quad \textrm{and} \quad \phi_{\textit{crit}}=\frac{\psi}{\psi+\psi^{4/5}}. \end{equation}

Comparison to these critical values greatly assists in determining which regime will occur for a given combination of $\psi$, $\phi$ and $R$. As is clearly shown by considering the critical points in figure 4, regime B2 cannot occur if $R > R_{\textit {crit}}$, regime C3 cannot occur if $R < R_{\textit {crit}}$, regime B3 cannot occur if $\phi > \phi _{\textit {crit}}$ and regime C2 cannot occur if $\phi < \phi _{\textit {crit}}$. Applying these conditions eliminates two possible regimes leaving two remaining possibilities (assuming that the A regimes have already been ruled out).

3.2.1. Boundary between regimes B2 and C2

This boundary should be considered if $R < R_{\textit {crit}}$ and $\phi > \phi _{\textit {crit}}$. For both regimes B2 and C2, upon substituting $\zeta _1=\phi$, the governing equations (2.15) and (2.18) simplify to

(3.5)\begin{equation} R_{BC2}=\sqrt{\frac{1-\phi}{\phi^5}}, \end{equation}

where $R_{BC2}(\phi )$ is the critical value delineating regime B2 ($R < R_{BC2}$) from regime C2.

3.2.2. Boundary between regimes B2 and B3

This boundary is important if $R < R_{\textit {crit}}$ and $\phi < \phi _{\textit {crit}}$. The boundary between regimes B2 and B3 is found by setting $\zeta _2=1$ in (2.16)–(2.17). This yields

(3.6)\begin{equation} R_{B23}=\sqrt{\frac{1-\zeta_1}{\zeta_1^{5/3}\psi^{2/3}(1-\phi)^{10/3}}}, \quad \textrm{for}\ (1-\phi)^{5/3}=\frac{\zeta_1^{5/3}}{\psi^{1/3}}+(\zeta_1-\phi)^{5/3}, \end{equation}

where $R_{B23}(\psi,\phi )$ is the critical value delineating regime B2 ($R < R_{B23}$) from regime B3.

3.3. Boundary between regimes B3 and C3

This boundary is important if $R > R_{\textit {crit}}$ and $\phi < \phi _{\textit {crit}}$. For both regimes B3 and C3, upon substituting $\zeta _1=\phi$, (2.16)–(2.17) and (2.19)–(2.20) simplify to

(3.7)\begin{equation} R_{BC3}=\sqrt{\frac{(1+\psi)(1-\phi)-\phi\psi^{4/5}}{\phi^5}}, \end{equation}

where $R_{BC3}(\psi,\phi )$ is the critical value delineating regime B3 ($R < R_{BC3}$) from regime C3.

3.4. Boundary between regimes C2 and C3

This boundary is important if $R > R_{\textit {crit}}$ and $\phi > \phi _{\textit {crit}}$. The boundary between regimes C2 and C3 is found by setting $\zeta _2=1$ in (2.19)–(2.20), yielding

(3.8)\begin{equation} R_{C23}=\sqrt{\frac{1-\psi^{1/5}(1-\phi)}{\psi(1-\phi)^5}}, \end{equation}

where $R_{C23}(\psi,\phi )$ is the critical value delineating regime C2 ($R < R_{C23}$) from regime C3.

3.5. Regime maps discussion

The links between the flow behaviour and a given set of values for $\psi$, $\phi$ and $R$ are encapsulated in the relations (3.3)–(3.8), but are better visually communicated via the use of regime maps. The maps in figure 4 show how flow behaviour varies with $\phi$ and $R$ at different constant values of $\psi$. While the behaviour when $\psi \geq 1$ is relatively straightforward (figure 4b,c), when $\psi < 1$, the regime map effectively conveys the eight different possible transitions between regimes that can occur upon changing $\phi$ or $R$.

As $R$ is a function of the box geometry, $R$ may take a fixed value in practical situations. In such cases, it is instructive to consider the regime boundaries as functions of $\phi$ and $\psi$, subject to constant $R$. Figure 5 shows regime maps for some example $R$ values: $R=0.1$; $R=1$, which is of the approximate order of magnitude for a room designed to be naturally ventilated; and $R=10$. These maps show that the existence of any of the regimes is plausible in a practical scenario, particularly when considering that, in the context of a room, opening or closing a window or door could change $R$ by an order of magnitude.

Figure 5. Regime maps in $\psi -\phi$ space for (a) $R=0.1$, (b) $R=1$ and (c) $R=10$. Regime A (purple), B (green) and C (blue) are delineated by dashed lines. Dotted lines delineate two-layer regimes (lighter shading) from three-layer regimes (darker shading).

4. Stratification predictions

Having expressed the governing equations (§ 2) for each regime and shown how to identify which regime prevails for a given set of values for $\psi$, $\phi$ and $R$ (§ 3), we can now calculate the stratification in a general emptying–filling box with two buoyancy sources. Figure 6 shows contours of interface heights, $\zeta _1$ and $\zeta _2$, and the normalised layer buoyancies,

(4.1a,b)\begin{equation} \widehat{g'_1}=\frac{g'_1H}{(B_1+B_2)^{2/3}} \quad \textrm{and} \quad \widehat{g'_2}=\frac{g'_2H}{(B_1+B_2)^{2/3}}, \end{equation}

as a function of $\psi$ and $\phi$ for $R=1$ and $R=10$. Figure 7 shows contours of the quantities

(4.2a{–}c)\begin{equation} \frac{g'_2}{g'_1}, \quad \hat{Q}_{n}=\frac{Q_n}{\left(B_1+B_2\right)^{1/3}H^{5/3}} \quad \textrm{and} \quad \frac{g'_2(1-\zeta_2)}{g'_1(\zeta_2-\zeta_1)}, \end{equation}

which are, respectively, the layer buoyancy ratio, the normalised flow rate through the box and the ratio of the relative importance of each buoyant layer in determining the flow rate. These contour plots reveal complex behaviours and highlight the importance of identifying the appropriate regime: note the cusps in the contours of $\zeta _1$ at the A3–B2 and B2–B3 boundaries in figures 6(a) and 7(a), the discontinuity in $\zeta _2$ across the A–B boundary in figure 6(c,d) and the stark differences in contour spacing and shape depending on the regime.

Figure 6. Contours of constant (a,b) $\zeta _1$, (c,d) $\zeta _2$, (e,f) $\widehat {g'_1}$ and (g,h) $\widehat {g'_2}$ as a function of $\psi$ and $\phi$ for (a,c,e,g) $R=1$ and (b,d,f,h) $R=10$. The regime map, with the same colour scheme as in figure 4, is shown in the background for reference. For results in two-layer regimes, $\zeta _2=1$ and $\widehat {g'_2}$ is undefined.

Figure 7. Contours of constant (a,b) $\hat {Q}_{n}$, (c,d) $g'_2/g'_1$ and (e,f) $g'_2(1-\zeta _2)/(g'_1(\zeta _2-\zeta _1))$ as a function of $\psi$ and $\phi$ for (a,c,e) $R=1$ and (b,d,f) $R=10$. The regime map, with the same colour scheme as in figure 4, is shown in the background for reference.

Despite the complexity shown in figures 6 and 7, we can still make general comments on the flow behaviour. Perhaps most significantly, the value of $R$ considerably affects both the shapes and magnitudes of the contours such that, in the context of a potential application, guidelines or advice for poorly ventilated buildings could be considerably different than for well-ventilated buildings. Also, most of the plotted quantities change with $\phi$ or $\psi$ more rapidly when in regimes B or C compared with regime A, thereby giving different routes for stratification control depending on regime.

5. Experimental campaign

Sixty-nine individual experiments with two saline sources were conducted to establish how the flow regime, the fractional interface heights ($\zeta _1$, $\zeta _2$) and buoyancy ratio ($g'_2/g'_1$) vary with $\psi$, $\phi$ and $R$. The campaign considered three values of $\psi$, and six values of $\phi$ and $R$, chosen to span regions of the parameter space where regime transitions were expected and subsequently observed to occur.

5.1. Experimental set-up

A clear acrylic box, of horizontal dimensions $50\ {\rm cm} \times 50\ {\rm cm}$ and vertical height $H_{p}=30$ cm, was submerged in a freshwater-filled, glass-sided visualisation tank (1.2 m height and 2.5 m by 1.2 m in plan), as sketched in figure 8. The tank was sufficiently large such that the ambient conditions surrounding the box could be considered constant for the entire duration of each experiment. The lower face of the box had one opening of diameter $\varnothing =3$ cm and eight openings of $\varnothing =5$ cm, and these openings were selectively blocked to vary the area available for outflow $a_t$. Experiments were conducted with either the single $\varnothing =3$ cm opening or 1, 2, 3, 4 or 8 openings of $\varnothing =5$ cm. The area of the inflow openings $a_b=573$ cm$^2$ was unchanged for all experiments and distributed across $22\times \varnothing =5$ cm and $20\times \varnothing =3$ cm openings on the upper face of the box. The ratio of the inflow and outflow opening areas exceeded $3.6$ for all experiments and, as confirmed by our observations, thereby ensured unidirectional flow at the openings and avoided mixing by the inflow (cf. Hunt & Coffey 2010).

Figure 8. Schematic of the experimental set-up showing the side view of the box inside the visualisation tank alongside the camera and lightbox (not to scale) and plan views of the upper and lower box faces (to scale). The view of the upper face shows the position of the two sources and the 42 inflow openings. The view of the lower face shows, in blue shading, the approximate size and position of the plume impingement regions (should the plume reach the lower face) and the locations of the outflow openings. These were selectively blocked, and the labels indicate which of the $\varnothing =5$ cm openings were used; for example, experiments with $3\times \varnothing =5$ cm openings had the openings labelled 1–3 unblocked with the remaining blocked.

To avoid concerns regarding the appropriate value of the coefficients $c_b$ and $c_t$ or how the effective area $A^*$ (2.4) is influenced by distributing the area over multiple openings, values of $R$ for each of the six opening configurations were determined in experiments with two sources of equal strength at the same elevation ($\phi =0$). For this benchmark condition, which intentionally excludes any effects due to differences in source parameters, the value of $R$ was determined from the steady height of the single interface via $R=(1-\zeta _1)^{1/2}(4\zeta _1^5)^{-1/2}$, which can be derived from theoretical consideration of this special case (cf. Linden et al. 1990) or by substituting $\phi =0$ and $\psi =1$ into the governing equations for regime A3 (2.13)–(2.14). There is some evidence that the value of $c_b$ or $c_t$ at an opening varies with the local flow conditions, particularly the density difference between the interior and exterior, such that $A^*$ and thereby $R$ might vary between experiments even if the opening geometry is unchanged (Hunt & Holford 2000; Radomski 2009; Vauquelin et al. 2017). These effects are not considered further herein, and the experiments are analysed taking $R$ as one of the six benchmark values, ranging from $R=7.6$ ($1\times \varnothing =3$ cm opening for outflow) to $R=0.40$ ($8\times \varnothing =5$ cm openings for outflow), as indicated in tables 1 and 2.

Table 1. Parameters (source strength ratio $\psi$, source height difference $\phi$, box resistance to emptying $R$ and source Richardson numbers $\varGamma _{0,1}$ and $\varGamma _{0,2}$ for the base and elevated plumes, respectively) and measurements (interface heights $\zeta _1$ and $\zeta _2$ and the layer buoyancy ratio $g'_2/g'_1$) for the 36 experiments with sources of nominally equal buoyancy flux. Entries for $\zeta _2$ and $g'_2/g'_1$ are given if two peaks in the gradient of the buoyancy profile could be identified. The entries in the ‘Expected’ and ‘Observed’ columns are the regimes predicted by the analytical model and observed in the experiments, respectively. If observations did not clearly show the existence of a single regime, then multiple regimes are listed in the order of best match (see Appendix B for discussion of ambiguous classifications).

Table 2. Parameters and measurements for the 33 experiments with unequal sources ($B_1\neq B_2$). Entries follow the same convention as table 1. By symmetry, the measurements at $\phi =0$ and $\psi =2.7$ (experiments 37, 43 and 49) also apply to $\phi =0$ and $\psi =1/2.7$; from this second perspective, the ‘Expected’ and ‘Observed’ regimes are A3 (instead of B3).

Each plume source had radius $b_0=2.5$ mm and was supplied from a constant head tank filled with saline solution of buoyancy $g'_0$, where $g'_0=g(\rho _0-\rho _e)/\rho _e$ and $\rho _0$ and $\rho _e$ denote the densities of the source saline and the freshwater environment, respectively. If the elevated source is within a buoyant layer, as for regimes C2 and C3, then the source buoyancy should instead be determined relative to the density of that layer, thereby reducing the source buoyancy flux compared with that calculated based on the definition of $g'_0$. However, this correction was not important herein given ${\rho _0-\rho _e\gg \rho _1-\rho _e}$, where $\rho _1$ denotes the density of the lower buoyant layer. Densities were measured using an Anton Paar DMA 4500 densitometer such that the source buoyancy was determined to within $1\,\%$. The source volume flux $Q_0$ was measured with a needle flow meter to within $4\,\%$. The source Reynolds numbers, based on the volume flux and diameter, were within the range 350–500, and the plumes were observed to be turbulent within two source diameters for all experiments, as expected from a nozzle based on the Cooper design (Hunt & Linden 2001). The base source was always on the upper face of the box, while the height of the elevated source was varied in $50$ mm increments, starting at the upper face until $50$ mm from the lower face. Virtual origin corrections were calculated for each source following the approach developed by Hunt & Kaye (2001) and the virtual sources (corresponding to zero volume and momentum flux) were 1.6–1.9 cm above the physical sources, depending on the source fluxes. The values for $H$ and $\phi$ were determined from the positions of the base and elevated virtual sources; for experiments where both sources were on the upper face, we neglect the small difference in the calculated virtual origins and report $\phi =0$. The source volume fluxes and buoyancies were chosen to set the buoyancy flux ratio as either $\psi =2.7, 1 \textrm { or } 1/2.7$ and such that each plume was near pure at source. For a pure plume, the scaled Richardson number

(5.1)\begin{equation} \varGamma_0=\frac{5}{8\sqrt{{\rm \pi} }\alpha_{p}}\frac{Q_0^2B_0}{M_0^{5/2}}=\frac{5{\rm \pi} ^2}{8\alpha_{p}} \frac{g'_0b_0^5}{Q_0^2}=1, \end{equation}

where top-hat profiles have been assumed (cf. Morton & Middleton 1973). The entries in tables 1 and 2 confirm that the plumes were nominally pure at source (within $15\,\%$ of $\varGamma _0=1$) and that the variation in $\phi$ and $\psi$ about each nominal value is small (within $3\,\%$).

Experiments were initiated by simultaneously supplying saline solution to the two plume sources. Measurements of the stratification were taken after interface height(s) and the average buoyancy in the layer(s) remained invariant for a period of 300 s and, as such, were regarded as representative of the steady state. The buoyancy distribution was determined using the dye attenuation technique to measure the dilution of methylene blue dye added to the source saline (Cenedese & Dalziel 1998; Allgayer & Hunt 2012). For these measurements, the visualisation tank was backlit uniformly by a lightbox comprising an array of closely spaced fluorescent tubes behind a translucent, diffusive acrylic sheet and the flows recorded with a CCD camera (JAI CV-M4+CL) positioned 3 m in front of the tank, aligned with the lower face of the box (figure 8). After making the buoyancy measurements, the lightbox was removed and the visualisation tank was lit by collimated light from a 35 mm slide projector for supplemental shadowgraph visualisations (with the same camera) on translucent paper placed on the front of the tank.

5.2. Diagnostics and interpretation

While the negatively buoyant saline plumes in the experiments descended vertically, further discussion and results are presented as if the plumes ascended for consistency and ease of comparison with the theoretical modelling.

Each row in figure 9 shows a shadowgraph visualisation, a time-averaged image of the buoyancy field and the gradient of the vertical buoyancy profile from a typical experiment. The vertical buoyancy profile was calculated by horizontally averaging the time-averaged buoyancy field across the width of the box, excluding the vertical strips occupied by the plumes or the physical source. Quantitative results for $\zeta _1$ and, for three-layer stratifications, $\zeta _2$ and $g'_2/g'_1$ are summarised in tables 1 and 2, and are based on the vertical buoyancy profile and its gradient. Note that the plotted profiles (e.g. figure 9c) correspond to the camera view and, thus, include the effects of parallax; the coordinate $z$ in these plots corresponds to the front of the box with the box base at $z=0$ and the top of the box at $z=H_{p}$ (the same height as the camera). Profiles are not shown for $z\lesssim 0.1 H_{p}$ as, in this region, light reflections from edges of openings in the box base contaminated the dye attenuation measurements; this did not impact any layer measurements. Interface positions were determined from the observed profiles following a parallax correction (Appendix A).

Figure 9. (a,d,g,j,m) Shadowgraphs, (b,e,h,k,n) time-averaged buoyancy fields and (c,f,i,l,o) vertical profiles of the buoyancy gradient. The value of $g'_*$ is such that the peak in the profile of the buoyancy gradient corresponding to $\zeta _1$ has a magnitude of unity, and the ‘No peak’ annotation on the profiles indicates the gradient continues to increase as $z/H_{p}\rightarrow 1$. Each row is an example of a different regime: (a–c) regime B2, experiment 29; (d–f) regime B3, experiment 45; (g–i) regime C2, experiment 69; (j–l) regime C3, experiment 57; (m–o) regime AB2, experiment 55.

Identifying the regime was straightforward for the majority of experiments as (i) the flow visualisations (e.g. figure 9a,b) and the gradient of buoyancy profile (e.g. figure 9c) clearly indicated whether there were one or two distinct buoyant layers, and (ii) the flow visualisations clearly showed which plume supplied a given layer. The first four rows of figure 9 show clear examples of regimes B2, B3, C2 and C3. The regime code AB2 is introduced to indicate a single buoyant layer supplied by both plume sources, as in figure 9(m–o). This two-layer regime exists as a result of plume-induced mixing that precludes a third layer forming when both plumes have similar buoyancies at $\zeta _1$. This could be considered the practical realisation of the theoretical special case when both plumes reach $\zeta _1$ with exactly the same buoyancy (§ 3.1.1) and, as such, this regime is expected for all cases where $\psi =1$ and $\phi =0$ (table 1). The regime observations for these five experiments are indicated by a single code in tables 1 and 2, as are other experiments with similarly straightforward classifications.

Other experiments were more difficult to classify because the shadowgraphs and buoyancy measurements showed characteristics of multiple regimes. While we did investigate metrics to help characterise ambiguous cases as one regime instead of another, it seemed more appropriate to indicate the uncertainty rather than arbitrarily introduce a distinction; the observation column for these experiments lists multiple regimes in the order of best match. The codes B2/B3, B3/B2, C2/C3 and C3/C2 reflect the similar appearances of a thin upper layer and of a horizontally propagating outflow current resulting from plume impingement with the top of the box. The codes AB2/B2, B2/AB2, AB2/B3 and B3/AB2 indicate experiments where mixing in the buoyant region disrupted flow features of regimes B2 or B3, but was not sufficient to result in the uniform layer required for an unambiguous regime AB2 classification. Appendix B has further discussion and examples of these difficult to classify cases, including experiment 19 with its unique AB3 classification.

5.3. Brief analysis

In general, the experimental observations match the analytical model. Every regime except A2 was observed, consistent with expectations based on the studied ranges of $\psi$, $\phi$ and $R$. The observations of regime AB2 for conditions where the two plumes did not reach the lower interface with the same buoyancy were not surprising; the model predicts $g'_2/g'_1\approx 1$ for these conditions and it was expected that plume-induced mixing could prevent such a weak stratification from forming or persisting. A more detailed quantitative comparison of the measurements and predictions follows in § 6.2, but we note that the agreement is generally good as the largest differences in the predicted and measured interface heights are only $0.05H$ for $\zeta _1$ and $0.16H$ for $\zeta _2$. There are some discrepancies between the observed and predicted regimes for some combinations of $\psi$, $\phi$ and $R$, but this is not unexpected given the deliberate choice to study parameter combinations that are predicted to be near regime boundaries.

The visualisations did show behaviours not included in the analytical model that may be responsible for some of the discrepancies. The shadowgraph in figure 10(a) shows evidence for two volume flux transfers across the upper interface resulting from plume-induced mixing. On the left, the ‘impinging’ plume – so named as it impinges on the upper layer – forms a turbulent mixing region on the upper interface. On the right, the ‘supplying’ plume – so named as it supplies the upper layer with volume flux $Q_{S}$ – ascends through the layer, impinges on the top of the box and forms a region of turbulent mixing that disturbs the upper interface. These two processes result in net transfers of fluid across the upper interface at rates $Q_{I}^*$ downwards and $Q_{S}^*$ upwards, as indicated in the idealised schematic in figure 10(b). Unsurprisingly, the significance of the mixing induced by the supplying plume appears to depend on the height of the upper interface; the shadowgraphs in figure 10(c,d) show less mixing on the interface for lower $\zeta _2$. While the former impingement behaviour was described and modelled by CL96 for the limiting case $k=0$, the additional mixing due to the impingement of the supplying plume was not, possibly because the mixing was not apparent in their experiments where both plume sources were at the base of the box and $\zeta _2$ was thus lower than observed in the present study.

Figure 10. (a) Shadowgraph visualisation (experiment 44) showing the base plume forming a region of turbulent mixing on the upper interface and showing the elevated plume impinging on the top of the box and disturbances on the upper interface. (b) Idealised model of fluid transfer across the upper interface with volume fluxes $Q_{I}^*$ and $Q_{S}+Q_{S}^*$. (c,d) Shadowgraph visualisations showing that the upper interface is more disturbed when it is closer to the impingement region at the top of the box (experiments 38 and 45, respectively).

The disruption of the upper interface by the observed mixing is also the primary source of uncertainty in the reported quantitative measurements (tables 1 and 2), as it causes the interface to span a region of finite thickness (in a time- and depth-averaged sense). While acknowledging that the significance of the mixing varies in each experiment, we estimate that the typical uncertainty in values of $\zeta _2$ is $4\,\%$. Disturbances also affect measurements of $\zeta _1$, and we estimate an uncertainty of $2\,\%$ for the majority of the experiments, but for those at lower values of $R$ where the interface is closer to the impingement regions, we estimate $3\,\%$ ($R=0.92$, experiments 19–24) and $4\,\%$ ($R=0.68$ and $0.40$, experiments 25–36). The value of the layer buoyancy ratio is sensitive to the heights over which the average layer buoyancies are calculated, and so it is also affected by the uncertainties in the interface positions, particularly when one of the layers is thin and/or an interface is not horizontal due to plume impingement; we estimate that the quantity $(g'_2-g'_1)/g'_1$ is within $5\,\%$ of the (implied) reported values for experiments with unambiguous regime classifications and within $10\,\%$ for those with ambiguous classifications.

6. Extended model and comparison with measurements

The analytical model for the three-layer regimes (§ 2) is now extended in order to improve both quantitative agreement and qualitative understanding of the flow behaviour. After presenting this extended model (§ 6.1), we show detailed comparisons with the experimental measurements (§ 6.2) and then briefly discuss how the model for the two-layer regimes could be similarly extended (§ 6.3).

6.1. Extension for the three-layer regimes

The extended model incorporates the vertical transport of fluid between layers caused by the impinging and supplying plumes, as sketched in figure 10(b), and is based on the same essential principles as the analytical model for the three-layer variants, namely, horizontal layers of uniform buoyancy driving a flow rate through the openings given by (2.5) subject to constraints imposed by the conservation of volume and buoyancy. Accounting for the mixing and the associated vertical transport does require modification of the conservation equations (2.1)–(2.3). For all regimes, continuity now requires

(6.1)\begin{equation} Q_n=Q_{11}+Q_{21}=Q_{S}+Q_{S}^*-Q_{I}^*. \end{equation}

Note that $Q_{21}=0$ in regime C3 and that whether the base or elevated plume is the supplying plume depends on the regime ($Q_{S}=Q_{12}$ in regime A3 and $Q_{S}=Q_{22}$ regimes B3 and C3). Conservation of buoyancy now requires

(6.2a{–}c)\begin{equation} \underbrace{g'_2=\frac{B_1+B_2}{Q_n}}_{\textit{All regimes}}, \quad \underbrace{g'_1=\frac{B_2+Q_{I}^* g'_2}{Q_{21}+Q_{I}^*}}_{\textit{A regimes}} \quad \textrm{and} \quad \underbrace{g'_1=\frac{B_1+Q_{I}^* g'_2}{Q_{11}+Q_{I}^*}}_{\textit{B & C regimes}}. \end{equation}

Given the transfer of buoyancy from the upper layer to the intermediate layer, $g'_1$ is greater than the average buoyancy of the impinging plume when it enters the layer (unlike in the analytical model where these buoyancies are equal).

Multiple theoretical and empirical studies have parameterised the vertical transport across a density interface due to turbulent mixing with an interfacial Froude number $\textit {Fr}$ relating the density difference to the scale and velocity of the impinging flow (e.g. Baines 1975; Kumagai 1984; Baines, Corriveau & Reedman 1993). Shrinivas & Hunt (2014b) developed a phenomenological model for the rate of transport across an interface based on a description of an impingement dome; given their predictions closely matched existing measurements and that this dome appears in our shadowgraph visualisations (figure 10), we applied their model to specify how $Q_{I}^*/Q_{I}$ varies with $\textit {Fr}$, where $Q_{I}$ is the volume flux of the impinging plume at the height of the upper interface ($Q_{I}=Q_{22}$ for regime A3 and $Q_{I}=Q_{12}$ for regimes B3 and C3). Here $\textit {Fr}$ and $Q_{I}$ are determined via numerical integration of the Morton et al. (1956) plume equations to account for the jet-like behaviour of the impinging plume as it crosses the intermediate layer. Further detail on the calculation of $Q_{I}^*$ is provided in Appendix C.

Unable to find a suitable existing phenomenological model to calculate the volume flux $Q_{S}^*$, we considered how the mixing resulting from the impingement of the supplying plume could be simply modelled. Looking at this flow in a general sense, namely, a stable stratification disturbed by turbulence, we expect some similarities with the Shrinivas & Hunt (2014b) model such that $Q_{S}^*$ scales with some representative volume flux and has some dependence on an interfacial Froude number. Specifying this Froude number dependence is not straightforward, particularly as the mixing at the interface is substantially different when the interface is very close to the plume impingement zone – where the precise structure of the impingement zone could be important – or when the interface is relatively far away. The latter case shares some features with grid turbulence experiments (e.g. Hopfinger & Toly 1976) and wind shear over a stratified water column (e.g. Mellor & Durbin 1975) where the turbulence generated away from a density interface causes vertical transport across it; Fernando (1991), Sullivan & McWilliams (2010) and Caulfield (2021) discuss the phenomena further. In these flows, the rate of vertical transport decays with increasing distance between the interface and the turbulence source, as indicated in our experiments (figure 10). We now propose a simple model to determine $Q_{S}^*$ based on these expected qualitative features.

Following existing models of plume impingement with a horizontal solid boundary (e.g. Kaye & Hunt 2007), we assume that the supplying plume creates a cylindrical impingement zone of height $\varDelta$. The radius $b_Z$ of the impingement zone and the volume flux $Q_Z$ into it are taken as the radius and volume flux of the supplying plume upon reaching the impingement zone, i.e. evaluated at $z=H-\varDelta$. Consideration of volume and kinetic energy fluxes into and out of the impingement zone shows that $\varDelta =b_Z/(2\gamma )$, where $\gamma ^2$ is the ratio of the kinetic energy flux leaving the impingement zone compared with the flux supplied to the zone by the plume; based on measurements of Ezzamel (2011), $\gamma \approx 0.8$. Taking the volume flux through the impingement zone as characteristic for the interfacial mixing, we scale $Q_{S}^*$ on $Q_Z$. We consider two types of flow behaviour depending on whether the interface is within the impingement zone or below. The distance below the impingement region is scaled on its thickness, yielding a relative distance $\beta =(H-\varDelta -h_2)/\varDelta$. To proceed, consider the following model, where, in the spirit of a simplified approach, we did not fine tune the form of the decay or introduce a prefactor (which is taken as 1 by default):

(6.3)\begin{equation} Q_{S}^*=\begin{cases} Q_Z, & h_2\geq H-\varDelta, \\ Q_Z\exp{(-\beta)}, & h_2< H-\varDelta. \end{cases} \end{equation}

Whilst the model is crude given the constant term within the impingement zone and the absence of any dependence on the stabilising density difference across the interface, we will show that incorporating this mixing model into the extended model yields predictions that closely match the experimental measurements. So, while we do not claim that the interfacial mixing itself has been well modelled, the predicted rates are sufficiently accurate in the context of the emptying–filling box that the inclusion of this simple model is worthwhile. This emptying–filling box context is significant, and we note that $\varDelta$ is a relatively small thickness compared with the box height and that $Q_{S}^*$ only appears in the conservation equations as part of the sum $Q_{S}+Q_{S}^*$, both factors that help to reduce the sensitivity of predictions for the interface heights on the details of this mixing model.

While we considered following the approach of CL96, namely, incorporating the effect of the stratification when calculating $Q_{S}$ and $Q_Z$, we show in Appendix C that this would change the calculation of $Q_{S}$ by less than $10\,\%$ for typical conditions. Given the uncertainty in the model for $Q_{S}^*$ (6.3) and, as noted, only the sum $Q_{S}+Q_{S}^*$ appears in (6.1), refining the calculation of the plume volume fluxes did not seem worth the added complexity, and they are determined using the pure plume equations (2.8), as in the analytical model (using $z=H-\varDelta$ rather than $z=h_2$ to determine $Q_Z$). Consistent with the calculations of $Q_{S}$ and $Q_Z$, the radius $b_Z$ was also determined using the pure plume equations, and upon substitution into $\varDelta =b_Z/(2\gamma )$, yields

(6.4a,b)\begin{equation} \underbrace{\varDelta=\frac{C_bH}{C_b+2\gamma}}_{\textit{Regime A3}} \quad \textrm{and} \quad \underbrace{\varDelta=\frac{C_b(H-k)}{C_b+2\gamma}}_{\textit{Regimes B3 & C3}}, \end{equation}

where $C_b=6\alpha _{p}/5$ is the coefficient for the normalised radius of a plume.

The stratification solution for a given set of $\psi$, $\phi$ and $R$ can be determined using the extended model using a numerical approach. Further details, including the numerical values of the various coefficients and the tolerance used for the plots herein, are described in Appendix C.

Note that a solution for a three-layer stratification is not necessarily expected at every combination of $\psi$, $\phi$ and $R$, and, as shown in figure 11, the shape of the regime maps in the extended model differs from the analytical model. The extended model does not dramatically change the regime map for the three-layer regimes, but it affirms quantitatively that stratifications at regime boundaries are sensitive to modelling of the mixing processes. In these examples, the solution space for regime C3 is larger because increasing $Q_{S}^*$ delays the approach of $\zeta _2\rightarrow 1$ while the spaces for A3 and B3 are smaller as increasing $Q_{I}^*$ causes a more rapid approach of $\zeta _2\rightarrow 1$ or $g'_2/g'_1\rightarrow 1$. While the extended model only provides a solution for a three-layer stratification, it does provide some insight on two-layer behaviour near to the boundary. For example, if $\zeta _2\rightarrow 1$ on approach to the boundary of a three-layer regime, then the two-layer counterpart would be expected on the other side, and regime AB2 would be expected following an approach where $g'_2/g'_1\rightarrow 1$. For boundary sections where $\zeta _2\rightarrow 1$ while $g'_2/g'_1\approx 1$ (key code 3, figure 11), there is some uncertainty regarding the regime on the other side; this is borne out by observations from experiments classified as regimes AB2/B2 or B2/AB2 (Appendix B).

Figure 11. Comparison of the regime maps predicted by the extended and analytical models for (a) $R=5$ and (b) $\psi =1/4$. The regime maps in the extended model, outlined in orange, have a different shape because of the mixing included in the model: $Q_{S}^*$ delays $\zeta _2\rightarrow 1$ and so increases the extent of the three-layer region while $Q_{I}^*$ pushes $\zeta _2$ towards $1$ and so reduces the extent. The regime boundaries in each model match more closely when $Q_{S}^*=Q_{I}^*$, this condition indicated by black dots. Numbered arrows, corresponding to the key, indicate the behaviour at the edge of the solution space.

6.2. Comparisons between models and experiments

Figure 12 compares the regime predictions of the models with the observations, and, overall, the models perform well. The analytical model predicts the correct regime, or either regime of an ambiguous classification, for the majority of the experiments and did not predict an incorrect regime away from the regime boundaries. Accordingly, the missed classifications were in a region of the parameter space where they might be expected. There are eight experiments for which the extended model improves on the analytical model by correctly identifying that a three-layer stratification would prevail (the orange dots labelled B3) or would be disrupted by mixing (the orange dots labelled AB2).

Figure 12. Regime classification of every experiment compared with the predictions of the models. The predictions of the analytical model are shown by the regime map (same colour scheme as figure 4). Black, unlabelled dots denote observations that match the analytical model. Black dots filled with white indicate ambiguous classification, of which the analytical model predicts one aspect. Orange dots indicate that the extended model (but not the analytical model) predicts the observations, with the white-filled dots indicating an ambiguous classification. Red dots indicate that neither model predicted the observed regime.

Figure 13 compares the measured and predicted buoyancy profiles for two example experiments and shows that there is good agreement. The experimental profiles are ‘staircase-like’ although the idealised sharp transition between layers can be smeared in practice by plume-induced mixing. The gradual transition across the intermediate buoyant layer (see $0.4\lesssim z/H_p \lesssim 0.7$ in figure 13a) has been previously described by CL96. The smearing of the upper interface (see $z/H_p\gtrsim 0.8$ in figure 13b) is expected given its close proximity to the impingement region at the top of the box and the surrounding layers have similar buoyancies. If the measured profiles are represented with idealised three-layer stratifications, it is clear that the models capture the interface locations and the layer buoyancy ratio, with the extended model offering improved accuracy.

Figure 13. Comparisons of measured and predicted buoyancy profiles for (a) experiment 39 and (b) experiment 56. The observed profile is shown in red, with the dotted line representing the idealised three-layer stratification based on the measurements; the apparent smearing of the interface is the result of parallax in addition to genuine smearing. The predictions of the analytical and extended models are shown by the dashed and solid lines, respectively.

Figure 14 compares the interface height predictions of both models with the experimental measurements. In general, both models perform well when predicting $\zeta _1$, the largest discrepancies being $0.05H$ and $0.06H$ for the analytical and extended models, respectively. The extended model is better at predicting whether an upper interface exists, for example, see the absence of the interface for experiment 55 (figure 14g at $\phi =0.17$) where the analytical model predicts an interface. The extended model is also better at predicting the height of the upper interface, and the largest discrepancy for $\zeta _2$ is $0.06H$ rather than $0.14H$ in the analytical model (discrepancies could only be evaluated if the model predicts $\zeta _2$ and an upper interface (ambiguous or not) could be identified from the experimental observations). Figure 15 compares the prediction of the buoyancy ratio $g'_2/g'_1$ from each model with the experimental measurements. Incorporating the effects of mixing into the extended model yields smaller predictions for the buoyancy ratio that better match the experimental measurements. While only the measurements from $R=7.6$ are presented, the results are similar for other $R$.

Figure 14. Comparisons of predictions and measurements of interface heights $\zeta _1$ (black) and $\zeta _2$ (red) with varying $\phi$ at the indicated values of $\psi$ and $R$. The predictions of the analytical and extended models are dashed and solid lines, respectively, and are not shown for $\zeta _2$ if a two-layer regime is predicted. Measurements are shown by dots, and a white centre indicates that the existence of the second interface is ambiguous. Error bars indicate the estimated uncertainty, but are not drawn if they overlap with the dots. The background colours show the regime predicted by the analytical model, following the same colour scheme as in figure 4.

Figure 15. Comparisons of predictions and measurements of the buoyancy ratio $g'_2/g'_1$ with varying $\phi$ at the indicated values of $\psi$ and $R$. The axis was set to aid comparison, and the value $g'_2/g'_1=7.36$ (experiment 42, C3/C2) does not appear in (c). The symbols and background colour scheme are otherwise the same as in figure 14.

It seems plausible that a more sophisticated model for $Q_{S}^*$ – for example, one which additionally considers the stabilising buoyancy difference between the upper layers – could improve the predictions of the extended model. However, we have not undertaken this exercise: precise modelling of the turbulent mixing is beyond the scope of this study of emptying–filling boxes, and any refinements would not significantly strengthen the demonstration that accounting for the vertical transport due to mixing is important for predictions. Furthermore, assessing whether a refined model offers significant improvement is difficult as the ambiguity in the regime classifications for many experiments precludes definitive conclusions.

6.3. Extension for the two-layer regimes

Significant mixing can also occur in two-layer regimes as shown by the flow visualisations of figure 16. Both plumes impinge on the top of the box and their outflows mix, preventing all of the outflow of the elevated plume from exiting the box as assumed in the analytical model. Additionally, the interface will potentially be close enough to the impingement regions at the top of the box such that fluid is transferred across the interface at a rate $Q_M$. The mechanism governing $Q_M$ is presumed to be similar to the one that drives $Q_{S}^*$ in the three-layer extended model, although both plumes could contribute to the mixing in the two-layer case. These processes are idealised in the sketch in figure 16(c) where $\sigma$ is the fraction of the flow from the exiting plume that exits (rather than being mixed into the buoyant layer).

Figure 16. (a) Shadowgraph and (b) buoyancy field visualisations from experiment 34 showing that some portion of the outflow of the exiting plume is mixed into the single buoyant layer, and that there is a visible disturbance on the interface. (c) Idealisation for an extended model of regime B2 incorporating mixing. The idealisations for regimes A2 and C2 are similar, but with the first plume as the exiting plume (regime A2) or with the second plume source above the interface (regime C2).

Similar to the extended model for three-layer regimes, solving for $\zeta _1$ and $g'_1$ requires a modification of the equations for the conservation of volume flux and buoyancy flux. We have included the modified governing equations for the two-layer regimes incorporating $Q_M$ and $\sigma$ in Appendix D, which may inspire development of models for $Q_M$ and $\sigma$. However, we have not specified functional forms of these newly introduced quantities herein because of the limited value this would add: the predictions of $\zeta _1$ – which allow calculation of $Q_n$ – in the two-layer regimes using the analytical model are already very good, as shown in figures 14 and 17. Indeed, further developments would only be worthwhile if, for example, a detailed description of the mixing within the buoyant layer was required. Additionally, it seems clear that accurately accounting for mixing between the two plume outflows to model $\sigma$ will depend on the horizontal locations of the plume sources and the upper openings, and is thus beyond the scope of a simplified modelling approach.

Figure 17. Comparisons of predictions of the analytical model and measurements of interface heights $\zeta _1$ (in black) and $\zeta _2$ (in red) with varying $\phi$ at the indicated values of $\psi$ and $R$. The symbols and background colour scheme are the same as in figure 14.

7. Discussion

There are additional refinements or extensions that could be made to our model. Increasing the number of buoyancy sources, while conceptually straightforward in principle, would dramatically increase the complexity of the analysis. Indeed, for three sources of arbitrary strength and elevation, there are 33 different regimes covering stratifications with two, three or four layers. This count is based on the same framework as the analytical model (§ 2); the effects of plume impingement are not included (which would be expected to increase the count), and the number of plumes that exit the box without forming a layer is between zero and two (i.e. one less than the number of plume sources). Non-idealised buoyancy sources could also be considered, specifically finite area sources, which would bridge the parameter space between the idealised point sources considered herein and the large area buoyancy sources considered by Livermore & Woods (2007). There are also questions germane to all wide emptying–filling boxes regarding how the distribution and geometry of openings and the local flow behaviour affects $c_t$, $c_b$ and, thus, $A^*$. However, these further developments are not necessary for the results of the model developed herein to be applied in practice.

While it is clear that our results could be applied to spaces with two buoyancy sources, it is perhaps less obvious that they could be used to assess whether a second source should be intentionally added to a space with a single existing source. As an example, consider a single heat source in a room that is relatively difficult to ventilate ($R=10$) and, consequently, has an interface close to the floor, as sketched in figure 18(a). By adding a second heat source with appropriate conditions, both the height of the lower interface and the flow rate of fresh air can be increased compared with the single plume reference case (figure 18b). Both of these effects could improve the comfort of occupants, particularly if the increase in $\zeta _1$ raises the lower interface above head height. The plots in figure 18(c) show that careful selection of the height of the second source is required: if too high ($\phi \gtrsim 0.70$, in this example), the second source would have no effect (as expected in regime C2) while, if too low ($\phi \lesssim 0.17$), the increase in the flow rate could be outweighed by the lowering of the first interface. The reduction of $\zeta _1$ implies the simultaneous increase of $g'_1$ – consider buoyancy conservation (2.3c) and the variation of the plume volume flux (2.8a) – which could be especially undesirable for occupant comfort; conversely, the increase of $\zeta _1$ results in $g'_1$ decreasing and, accordingly, a cooler layer.

Figure 18. (a) A single plume source in a naturally ventilated room where $R=10$ creates a layer at $\zeta _1=0.364$ and drives the flow rate $Q_{\mathit {ref}}$. (b) Adding a stronger, elevated plume source ($\psi =3$, $\phi =0.3$) creates a new stratification where the first interface is higher and which results in an increased flow rate compared with the single plume scenario. (c) Variation of $\zeta _1$ (black) and $Q_n$ (red) with $\phi$ for $\psi =3$ and $R=10$. The vertical dotted line indicates $\phi =0.3$. The dashed lines indicate the values in the single plume scenario, showing that there is a range of $0.17\lesssim \phi \lesssim 0.71$ where both $\zeta _1$ and $Q_n$ are increased.

Adding a second buoyancy source to enhance the natural ventilation of a space is particularly appealing if the heat would otherwise be wasted. It is intriguing to imagine the heat rejected from an air-conditioning system cooling one part of the building, such as a computer centre, being used to enhance conditions in another room. While the exhaust of the air-conditioning system would not be directly cooling the room (indeed the heat input to the room would have increased), it could nevertheless improve human comfort by increasing the cooling draught in the space, as demonstrated by the preceding example. There is doubtless potential for clever designs to incorporate natural ventilation strategies using two plumes to reduce overall energy consumption of a building and/or to enhance their habitability.

It is also worth highlighting that the model could be applied in cases where both buoyancy sources are at the same physical height, but of a different size. For a range of source conditions, virtual origin corrections could be appropriately used to represent the physical sources as point sources which, in general, would be at different heights.

8. Conclusions

We have addressed the emptying–filling box problem with two localised buoyancy sources of arbitrary elevation and strength, with an emphasis on the previously unconsidered case of sources at distinct elevations. In order to predict the steady stratification created by the turbulent plumes arising from any pair of these sources, we identified a rich array of two- and three-layer stratifications across six possible flow regimes. The analytical model based on simplified governing equations for these regimes reveals how the prevailing regime and stratification properties varies with three key parameters: the relative strengths of the plumes, the height difference between their sources and the resistance of the box to emptying.

The experimental campaign to assess 69 parameter combinations, chosen to ensure that regime boundaries were crossed, broadly confirmed the predictions of the analytical model. However, flow visualisations showed that mixing resulting from plume impingement on either a density interface or on the top of the box could be significant, and, for some combinations of the governing parameters, result in a seventh regime. These two mixing behaviours were thus incorporated into an extended model that provides predictions that more accurately match experimental observations. For certain stratification properties, such as the height of the lowest layer, the predictions of the analytical model are in excellent agreement with measurements, and so, depending on the application, including the effects of mixing will not necessarily be required.

Our analysis revealed two previously undocumented behaviours that, while particularly significant when a source is elevated, can occur even if the sources are at the same elevation. First, two plumes can form a two-layer stratification (rather than the oft-assumed three-layer stratification) if the openings in the box are sufficiently large. We expect this result to generalize such that $n$ plumes in an emptying–filling box could form a number of layers between $2$ and $n+1$. Second, plume impingement on the top of the box can result in significant mixing across a density interface if it is sufficiently close.

The work relies on investigations into a variety of fundamental problems in fluid mechanics, particularly the quantification of entrainment across a density interface. Whilst the impingement of a jet-like flow on a density interface is fairly well understood, the vertical transport across a density interface near to a plume impinging on a horizontal boundary is not. While we demonstrated that our proposed simplified model for this transport is acceptable in the context of improving the predictions of emptying–filling box behaviour, the flow is worthy of further study in its own right. Similarly, an investigation into the interaction between the lateral outflow currents formed by the plumes and the openings at the top of the box would be enlightening.

The analytical model predicted the height of the lower interface to within $5\,\%$ for every experiment, which implies that the flow rate through the box can be determined to within $10\,\%$. From a practical perspective, these quantities will typically be the most important when assessing, for example, human comfort in naturally ventilated buildings, and so the simpler analytical model will likely be suitable for design purposes. We demonstrated one design in an example where a heat input is used to provide a cooling effect and discussed other possible applications. If accounting for the redistribution of fluid between the two buoyant layers is important, then the extended model should be used. This model could be improved or extended ad infinitum – and we have suggested several salient routes – although it is not clear, in the context of the cited application, that further refinement is necessary.