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).

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 Reference Hunt and Coffey2010).

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. Reference Linden, Lane-Serff and Smeed1990) 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 Reference Hunt and Holford2000; Radomski Reference Radomski2009; Vauquelin et al. Reference Vauquelin, Koutaiba, Blanchard and Fromy2017). 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 Reference Hunt and Linden2001). 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 (Reference Hunt and Kaye2001) 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 Reference Morton and Middleton1973). 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 Reference Cenedese and Dalziel1998; Allgayer & Hunt Reference Allgayer and Hunt2012). 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 Reference Cooper and LindenCL96 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.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 Reference Baines1975; Kumagai Reference Kumagai1984; Baines, Corriveau & Reedman Reference Baines, Corriveau and Reedman1993). Shrinivas & Hunt (Reference Shrinivas and Hunt2014b) 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. (Reference Morton, Taylor and Turner1956) 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 (Reference Shrinivas and Hunt2014b) 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 Reference Hopfinger and Toly1976) and wind shear over a stratified water column (e.g. Mellor & Durbin Reference Mellor and Durbin1975) where the turbulence generated away from a density interface causes vertical transport across it; Fernando (Reference Fernando1991), Sullivan & McWilliams (Reference Sullivan and McWilliams2010) and Caulfield (Reference Caulfield2021) 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 Reference Kaye and Hunt2007), 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 (Reference Ezzamel2011), $\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 Reference Cooper and LindenCL96, 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 Reference Cooper and LindenCL96. 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.