2.1 Apparatus

The experiments were carried out in a deep, square box of internal width $W=0.3$  m and height 0.55 m, with a permeable sintered brass plate as the horizontal bottom boundary for the convection domain (figure 1). Below the permeable plate was a chamber filled with saline solution that was in contact with the permeable plate and slightly pressurised from an external supply reservoir. This chamber also housed an inner manifold which supplied relatively fresh water to small circular areas of the permeable plate through a square array of cylinders, each 25 mm diameter, separating the freshwater from the saline water. In most experiments the inner manifold supplied 16 circular regions, centred a distance $L=75$  mm apart in a $4\times 4$ array, which in turn was symmetrically arranged within the width of the box. A few experiments were carried out with a similarly symmetric $3\times 3$ array of fresh inputs. Peristaltic pumps coupled with flow integrators (to smooth the pump outputs) were used to supply constant volume fluxes (as measured by flowmeters) of de-aerated, filtered saline solution and freshwater to the manifolds from external reservoirs. Both reservoirs were kept well mixed, ensuring constant input salinities throughout a run.

Figure 1.

Schematic of the laboratory apparatus showing the forcing arrangement, with regions of dense saline (blue) and buoyant freshwater (red) input at the permeable base, and vertical cylinders supplying freshwater (at flux  $Q_{f}$ ) from the inner manifold to the permeable plate (black). The sketch shows the water depth for the deep case and the forcing array for the wavelength $L=75$  mm,  $W/L=4$ . The separate planform of the density distribution at the base shows the symmetric arrangement of forcing with respect to the sidewalls.

The approach is the same as used to study ‘pure’ buoyant plumes (Baines & Turner Reference Baines and Turner1969), where a finite but small and dynamically insignificant volume flux is required to deliver the buoyancy. For the present experiments the effect of the small volume input on the bottom boundary layer is evaluated in § 2.3. In order to limit the final water depth $H$ and allow the convection to approach a true statistically steady state, water was also removed from the top by placing the inlet of a Venturi pump on the centreline of the tank at the desired final water depth. This sink withdrew water and salt, but represented zero buoyancy flux ( $g\unicode[STIX]{x1D70C}^{\prime }(Q_{s}+Q_{f})$ ), as the density of solution withdrawn is identical to that of the bulk of the box ( $\unicode[STIX]{x1D70C}^{\prime }=0$ ). The location of the sink at the top also ensures that the withdrawal has negligible influence on the density of the bulk of the box and the individual saline and freshwater buoyancy fluxes at the base. The small volume throughput and associated small mean vertical velocity $\overline{w}\approx (Q_{f}+Q_{s})/W^{2}$ (which will be seen to be around 20 times smaller than the average convection velocities in the interior) also ensure that the sink has no significant influence on the convective velocity field. Adequate quasi-steady states were reached before the water surface reached the sink, as water depth increased very slowly. However, the increasing depth did not suit long measurement periods required to obtain time-averaged flow statistics. Hence the measurements reported here are from the final, true equilibrium states having constant  $H$ . Experiments were carried out with $H=150$  mm ( $H/L=2$ ) and $H=375$  mm ( $H/L=5$ ).

Conductivity measurements were made using a PME MicroScale Conductivity and Temperature probe (Model 125) attached to a traversing mechanism, in turn driven by a Smart Motor and controlled by a computer program. Vertical profiles of conductivity were taken at regular intervals, excepting when high-resolution time records were required, in which case the probe was placed at a fixed location for the desired duration of the record. The location of the probe was varied in one horizontal ( $x$ ) direction in order to obtain profiles above both saline and freshwater input areas. Direct density measurements of samples, withdrawn from various heights through a tube at $(x,y)=(75~\text{mm},112.5~\text{mm})$ , were made using an Anton Paar digital densimeter (Model DMA5000; precision of $\pm 1\times 10^{-6}~\text{g}~\text{cm}^{-3}$ ).

Two-component velocity fields ( $u,w$ ) were obtained by particle tracking velocimetry (PTV). For this the flow was seeded with 125– $250~\unicode[STIX]{x03BC}\text{m}$ diameter particles of Pliolite resin and illuminated from the side in a vertical plane by a light sheet from a 250 W white light projector beam passing through a 5 mm slit at the wall of the box. Input salinities (particularly of the fresher input) and flow rates were chosen such that the bulk density in the convection domain in the final equilibrium state (as predicted theoretically in § 2.3) was approximately $1021~\text{kg}~\text{m}^{-3}$ , matching the density of the Pliolite. Particle motion was recorded, at 25 f.p.s., by a high-definition video camera mounted 4 m from the box, the movie later decomposed into images, particles identified and tracked, and their velocities estimated using the Streams package (Nokes Reference Nokes2014). The particle-centred velocities were then interpolated onto a regular grid at each time step, yielding a time series of two components of the Eulerian velocity field on the measurement plane. In separate runs dedicated to qualitative flow visualisation the same light sheet was used to illuminate fluorescein dye, which was introduced to the freshwater input in trace concentrations.

2.2 Experimental parameters

The governing parameters for these isothermal experiments may be written as the solutal Rayleigh number $Ra_{s}$ (based on the horizontal forcing scale $L$ and solute diffusivity  $\unicode[STIX]{x1D705}_{s}$ ), the Schmidt number $Sc$ , the dimensionless depth and width of the convection domain, a Peclet number $Pe$ based on input flow rate, and the ratio $R_{Q}$ of destabilising to stabilising volume fluxes:

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is the applied density difference between the two input solutions of salinity $S_{f}$ and $S_{s}$ , and $Q_{f}$ , $Q_{s}$ are the corresponding area-integrated imposed volume input rates of the relatively fresh and saline inputs, respectively. The width $W$ was designed to be an integer multiple $n$ of the forcing wavelength  $L$ . Note that the Rayleigh number in the Rossby problem is often based on only one half of the forcing wavelength (heating near one end and cooling near the other end over a box length  $L$ ), whereas $Ra$ is here based on the full wavelength of the forcing and is therefore eight times larger than the more common horizontal Rayleigh number. Given that both the input salinities and the volume fluxes are also imposed, a flux Rayleigh number can also be defined as $Ra_{F}=gQ_{f}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}L^{2}/\unicode[STIX]{x1D70C}_{o}\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}_{s}^{2}=Ra_{s}Pe$ .

The conditions for a series of experiments are summarised in table 1. In this series each experiment commenced with 50 mm of solution above the permeable base at a density greater than the expected final density in the bulk of the domain. The water freshened as the box very slowly filled. Before quasi-equilibrium density was reached the smaller of the input densities drove bottom-to-surface convection. The average salinity of the bulk of the domain eventually approached a constant or equilibrium value, and convection was then maintained by the differential buoyancy inputs. Once the water surface reached the level of the withdrawal sink, the experiment continued with a constant water depth. The water depth could later be increased by raising the sink to a higher level and allowing a further period of equilibration.

Table 1.

Summary of experimental conditions in runs with $n=4$ : imposed density difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f}$ ( $\pm 0.005$ ), volume flow rates $Q_{s}$ ( $\pm 0.3\times 10^{-6}$ ) and $Q_{f}$ ( $\pm 1\times 10^{-6}$ ) of saline and fresher inputs, measured density difference $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }=\bar{\unicode[STIX]{x1D70C}}^{\ast }-\unicode[STIX]{x1D70C}_{f}$ ( $\pm 0.05$ ) between the bulk and the fresher input at equilibrium, buoyancy time scale $t_{b}=[(H/g)(\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast })]^{1/2}$ , Peclet number and Rayleigh number. All runs had $\unicode[STIX]{x1D705}_{s}\approx 1.5\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$ , $Sc\approx 760$ and measurements were generally taken for water depths of $H\approx 2L$ and $H\approx 5L$ . Uncertainty in $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }$ is given by one standard deviation, after steady state was reached.

2.3 Theoretical equilibrium states

For constant volume flow rates of two solutions through a permeable base, into a box initially containing a small volume $V_{0}$ of density $\unicode[STIX]{x1D70C}_{0}$ , the total mass of solution in the box at time $t$ is $M=\unicode[STIX]{x1D70C}_{s}Q_{s}t+\unicode[STIX]{x1D70C}_{f}Q_{f}t+\unicode[STIX]{x1D70C}_{0}V_{0}$ . For a density $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}\,f(S)$ (where $f$ is an arbitrary function, $S$ is salinity and $\unicode[STIX]{x1D70C}_{0}$ is a reference density) the mass of salt divided by the total mass of solution gives a bulk salinity  $\bar{S}$ :

Taking $\unicode[STIX]{x1D70C}_{f}\leqslant \unicode[STIX]{x1D70C}_{0}<\unicode[STIX]{x1D70C}_{s}$ , the solution becomes increasingly saline and approaches the large-time limit (denoted throughout by superscript ^∗)

at times $t\gg t_{eq}=\unicode[STIX]{x1D70C}_{0}S_{0}V_{0}/\unicode[STIX]{x1D70C}_{s}S_{s}Q_{s}$ (assuming $Q_{s}<Q_{f}$ ). In this limit the added mass of salt becomes large compared to the initial mass of salt. Thus the average salinity is constant at large times. If $Q_{s}\ll Q_{f}$ , the equilibrium salinity (2.4) is $\bar{S}^{\ast }\approx \unicode[STIX]{x1D70C}_{s}S_{s}Q_{s}/\unicode[STIX]{x1D70C}_{f}Q_{f}$ and is much closer to the salinity of the fresher input than to that of the saline input. The salinity (2.4) is identical to that expected for a steady state in which the same saline and freshwater are injected continuously into any well-mixed box, with or without a matching volume flux from a sink. Hence the addition of the sink that withdraws the mixture in the convection experiments does not influence the bulk salinity and this, in turn, implies it has no influence on the input buoyancy fluxes.

The buoyancy fluxes into the box are determined by the bulk salinity predicted above. Once the bulk has reached the large-time equilibrium salinity the individual buoyancy fluxes carried by the fresher and more saline inputs become, respectively,

where $\bar{\unicode[STIX]{x1D70C}}^{\ast }=\unicode[STIX]{x1D70C}_{0}\,f(\bar{S}^{\ast })$ is the equilibrium density. Defining $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }=\bar{\unicode[STIX]{x1D70C}}^{\ast }-\unicode[STIX]{x1D70C}_{f}$ , and recalling $\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{f}=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , the destabilising buoyancy input is $B_{f}=g(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }/\unicode[STIX]{x1D70C}_{0})Q_{f}$ and the stabilising input is $B_{s}=g[(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C})/\unicode[STIX]{x1D70C}_{0}]Q_{s}$ . The net boundary buoyancy input is

In the case of a linear equation of state $\bar{\unicode[STIX]{x1D70C}}^{\ast }=(\unicode[STIX]{x1D70C}_{s}Q_{s}+\unicode[STIX]{x1D70C}_{0}Q_{f})/(Q_{f}+Q_{s})$ , $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }/\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=Q_{s}/(Q_{s}+Q_{f})$ and (from (2.6)) the net boundary buoyancy input $B$ in the equilibrium state vanishes. This is a good approximation across the whole range of source salinities used in the present experiments. The accuracy of the approximation is further assured in the bulk of the fluid, where property differences are small compared to those in the dense boundary layer on the base. Thus the net buoyancy flux through any level vanishes at large times $t\gg t_{eq}$ . During the evolution to the steady state the assumption that the bulk remains well mixed requires the equilibration to be slow compared to the convective overturn, thus $t_{eq}\gg t_{b}$ , where $t_{b}=[(H/g)(\unicode[STIX]{x1D70C}_{f}/\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast })]^{1/2}$ is the buoyancy time scale. In the steady state having zero net buoyancy flux through any level, the convective overturn need only be rapid compared to the very slow mean vertical velocity imposed by the finite volume throughput, a condition easily satisfied in the experiments.

For the linear equation of state the individual buoyancy inputs become

and in the limit $Q_{f}\gg Q_{s}$

In this limit the corresponding Nusselt number describing the buoyancy throughput relative to a purely diffusive lateral transport can be written as

where $Ra_{F}^{\ast }$ is a flux Rayleigh number

based on the actual destabilising buoyancy input (rather than on the external variables in (2.1)).

The bottom boundary conditions in the experiments are of mixed kind, as both densities and fluxes are imposed at the permeable base. Furthermore, both diffusive and advective fluxes may be significant at the boundary. The role of imposed fluxes is minimised if the salinity of the saline input is very large, so that the flow is characterised by $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\ast }\ll \unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , which leaves the saline input as an imposed flux condition $Q_{s}S_{s}$ . The effect of a non-zero vertical velocity through the permeable boundary can be assessed by recalling that scaling analysis of the Rossby problem (Rossby Reference Rossby1965, Reference Rossby1998) treats the boundary layer as governed by a balance of vertical advection and diffusion, so that (in the absence of a volume flux through the boundary) the flow has a vertical (downward) velocity $w_{bl}\sim \unicode[STIX]{x1D705}/\unicode[STIX]{x1D6FF}$ , where $\unicode[STIX]{x1D6FF}$ is the boundary layer of thickness. An approximation for the boundary layer thickness can then be found by assuming simple horizontal momentum balances in the boundary layer, either viscous buoyancy (Rossby Reference Rossby1965, Reference Rossby1998) or inertial buoyancy (Gayen et al. Reference Gayen, Griffiths and Hughes2014, equation (5.3)). The viscous balance, expected to give the smallest $w_{bl}$ and valid only at small values of $Ra$ , taken together within continuity, leads to $\unicode[STIX]{x1D6FF}/L\sim Ra^{-1/5}$ , which gives $\unicode[STIX]{x1D6FF}\sim 1$  mm for the Rayleigh numbers in the present compositional convection. The normal velocity imposed at the permeable base in the stabilised areas, $\overline{w_{s}}\approx Q_{s}/W^{2}$ , reduces the role of diffusion in governing boundary layer thickness and is likely to thicken the boundary layer over that given by the Rossby scaling when $\overline{w_{s}}>w_{bl}$ , or $Pe>n^{2}Ra^{1/5}$ . For the experiments this equates to $Pe>10^{6}$ and implies that for the experimental conditions (table 1) the boundary layer is unlikely to be influenced by the volume input. However, the precise nature of the boundary conditions is not important to our conclusions.

3.1 Equilibration of density

Experiments began with a small volume of saline solution in the box and a net destabilising buoyancy flux through the base, causing the density in the bulk of the water to decrease with time. The water depth increased slowly until the surface reached the sink. The system reached an equilibrium state in which the density at any location showed no further secular trend, fluctuating about a mean value (figure 2). The adjustment was approximately exponential with a time scale close to the predicted $t_{eq}$ , and equilibrium was reached by $t\approx 3t_{eq}$ . In the equilibrium state (for $Q_{s}\sim 0.1Q_{f}$ ) the density difference between the bulk and the fresher input, $\unicode[STIX]{x1D6FF}^{\ast }=\bar{\unicode[STIX]{x1D70C}}^{\ast }-\unicode[STIX]{x1D70C}_{f}$ , was 5 %–9 % of the applied density difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ (table 1), and slightly fresher than predicted by (2.4).

Figure 2.

Evolution of densities relative to the large time mean in the interior: (a) shallow-water case (exp. 20), plotted with a linear ordinate; and (b) deep-water case (exp. 31) plotted on a logarithmic ordinate. Symbols show density from samples (taken at $z/L=1$ and $z/L=3.7$ in the deep case and $z/L=0.5$ and $z/L=1.5$ in the shallow case, circles and triangles) and the continuous grey line shows in situ conductivity (at $z=3.7L$ in the deep case and $z/L=1.5$ in the shallow case). Densities are normalised by the overall change from initial to final local interior density and plotted as a function of elapsed time $t$ from the beginning of the experiment. Time is normalised by the equilibration time scale $t_{eq}=(\unicode[STIX]{x1D70C}_{0}S_{0}V_{0})t/(\unicode[STIX]{x1D70C}_{s}S_{s}Q_{s})$ . Fluctuations of $\unicode[STIX]{x1D70C}(t)$ about the mean equilibrium value $\bar{\unicode[STIX]{x1D70C}}^{\ast }$ are both positive and negative and therefore in (b) the absolute value is plotted. An empirical fit of $\unicode[STIX]{x1D6FD}\text{e}^{-t/t_{eq}}$ and the 98 % adjustment level are shown as dot-dashed and dashed lines respectively.

Figure 3.

Vertical profiles of density from two deep-water cases ( $A=5$ ; exps 25, 26; red and black) and two shallow-water experiments ( $A=2$ ; exps 27, 28, cyan and magenta) in the equilibrium states. Density $\unicode[STIX]{x1D70C}^{\ast }$ is shown relative to the density $\unicode[STIX]{x1D70C}_{top}^{\ast }$ at the top of the profile and normalised by the applied density difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ . Profiles cover approximately $2/3$ of the water column to 5 mm from the base.

The bulk of the box was always stably stratified (figure 3). As the box freshened with time the density differences in the box increased, at least in a lower region $z<L$ above areas of saline input. There the bottom boundary layer thickness (defined as the region containing 90 % of the top-to-bottom density difference, assuming a bottom density equal to the saline input $\unicode[STIX]{x1D70C}_{s}$ ) was $\unicode[STIX]{x1D6FF}\approx 0.1L$ , somewhat thicker than estimated for the viscous boundary layer regime. However, the boundary layer was too thin for much of its profile to be measured. The corresponding buoyancy frequency $N=[-(g/\unicode[STIX]{x1D70C}_{0})\,\text{d}\unicode[STIX]{x1D70C}^{\ast }/\text{d}z]^{1/2}$ is over an order of magnitude less than the depth-averaged buoyancy frequency $N_{0}=[(g/\unicode[STIX]{x1D70C}_{0})(\unicode[STIX]{x1D70C}_{s}-\unicode[STIX]{x1D70C}_{top}^{\ast })/H]^{1/2}$ throughout most of the depth ( $1\leqslant z/L\leqslant 3.5$ in the deep cases), and is also consistent with the expectation that the average in the boundary layer will be $(N/N_{0})^{2}\sim H/\unicode[STIX]{x1D6FF}\sim 10^{2}$ .

Figure 4.

Images of fluorescein dye in the shallow case (exp. 31b, $A=2$ , $n=4$ ) using illumination of a vertical plane aligned with a row of freshwater input regions and after the densities had been at equilibrium for at least 80 min. Images were taken 60 s apart, approximately 10 min after the dye was added to the freshwater supply manifold.

Figure 5.

Images of Fluorescein dye injected with the freshwater in the deep case (exp. 18, $A=5$ , $n=4$ ) showing the tracer distribution in the same plane as in figure 4. Images show the flow 140 min after the experiment began and had been in equilibrium for at least 30 min, and approximately 10 min after the dye was added to the freshwater supply manifold. Frames are at intervals of 60 s.

3.2 Flow visualisation

When trace amounts of fluorescein or food dye were added to one of the inputs after the system reached equilibrium, the bottom boundary layer flow could be seen. Over a time much shorter than the equilibration time, tracer in the saline input tended to slowly flood a very thin layer on the permeable plate, of order 2 mm thick. From this sublayer at the bottom of the boundary layer the tracer was carried into the base of plumes as they ascended from the areas of freshwater input. It was not possible to continue the tracer input and retain visible concentration differences for times comparable to the equilibration time, meaning that the full thickness of the boundary layer revealed by the conductivity profiles could not be observed. Addition of fluorescein dye to the fresher input (figures 4 and 5) again showed convergent horizontal flow on the base in the outer parts of each circular area of fresh input, but convective up-flow to heights $z/L\sim 1$ over most of the circular patch. The up-flow formed coherent and somewhat turbulent plumes immediately above each of the regions of fresher input. In all shallow-water cases the plumes extended through the full depth (figure 4). Down-flow between the plumes could be seen once tracer began recycling from the top. In the deep-water cases for $n=4$ (figure 5) the plumes broke down at a height $1<z/L<2$ , above which the tracer showed more disorganised and time-varying motions that included a large box scale of overturning, in which individual plumes on the forcing scale were identifiable only intermittently. No such transition from a fixed array of plumes to large scales was seen for $n=3$ , where the square array of coherent plumes was maintained even in deep water ( $A=5$ ).

Figure 6.

Mean 2D velocity fields and streamfunction at two times for the shallow case $A=2$ (exp. 27) on a vertical plane through a row of the freshwater inputs. Fields are temporal averages over $3t_{b}$ beginning at $t=50$  min (a) and $t=80$  min (b). The density reached equilibrium by $t=30$  min. Velocity scale (colourbar) is in $\text{mm}~\text{s}^{-1}$ ; 2D streamfunction contours are at intervals of $20~\text{mm}^{2}~\text{s}^{-1}$ .

Figure 7.

Mean 2D velocity fields and streamfunction at four times for the deep case $A=5$ (exp. 13), on the same vertical plane as shown in figure 6 and averaged over 100 s (approximately $40t_{b}$ ), beginning at (a) 40 min, (b) 45 min, (c) 80 min and (d) 85 min into the experiment. The density reached steady state before $t=30$  min. The magnitude of velocity vectors is shown by the colour scale in $\text{mm}~\text{s}^{-1}$ ; 2D streamfunction is contoured at intervals of $30~\text{mm}^{2}~\text{s}^{-1}$ .

3.3 Velocity fields and wavenumbers

Velocity measurements confirm the equilibration of the flow. The temporal average $\bar{u}$ , $\bar{w}$ of the two-component velocity field over the vertical planar area $A$ of the light sheet was calculated over three buoyancy time scales $t_{b}$ , and from that the mean flow kinetic energy $\bar{E_{k}}$ ( $\text{J m}^{-1}$ ) was found,

This measure of the kinetic energy is derived only from the resolved scales, noting that the turbulent kinetic energy of the small scales is not evaluated due to the coarse spatial resolution of the velocity field, and the average over a time window effectively provides the energy of the larger scales of motion. The mean kinetic energy is initially large as a result of net input of buoyancy and decays on a time scale of order $10^{3}$  s, which is consistent with the time scale found for density adjustments. At large times the secular trend vanishes. However, the equilibrium continues to show significant variability, reflected in both large standard deviations at time scales less than 300 s and variability between averages over 300–600 s windows, on the order of 50 % of the long-term average kinetic energy. It was confirmed by tracer and PTV results that this variability includes low-frequency reversals of the overturning on the scale of the whole domain.

The velocity fields present a picture of the flow behaviour consistent with that from advection of dye. For the shallow case, as illustrated by 100 s averages in figure 6, the vertical velocity shows maxima in each of the plumes above the regions of destabilising boundary flux, and these maxima are coherent throughout the depth. Between the plumes the flow is predominantly downwards. Although there is some low-frequency fluctuation of the velocity field (as seen over the 55 min between the fields in figure 6) the plume scales of the flow pattern in the shallow case are stationary.

For the deep case (figure 7) velocity maxima and coherent vertical plumes can be associated with only some of the fresher inputs at any one time, only intermittently at any fresh input region, and then only at $z<2L$ . Above this level the velocity fields reveal circulation dominated by the scale of the box. This large scale of overturning frequently involves an asymmetry, with regions of faster, narrower flow that can be either upwards or downwards, and broader return flow. Velocities also vary with time, as illustrated by the fields in figure 7: the maximum velocity doubles in the 5 min elapsed between figures 7(c) and 7(d), and the large-scale motion reverses direction (clockwise to anticlockwise) in the 40 min elapsed between figures 7(a,b) and 7(b,c).

Wavenumber spectra were calculated from a fast Fourier transform of the vertical component $w(x)$ of velocity in the illuminated vertical plane at a given time and at several heights. These velocity–wavenumber spectra were computed from frames at ten second intervals and averaged over a ten minute record, and the results in figure 8 are shown (for mid-depth and upper region) normalised by the maximum in power spectral density across all heights. Only wavenumbers $0\leqslant k\leqslant 25$ , corresponding to length scales larger than 40 mm, are resolved, as the particle velocity information was binned over 10 mm windows. There was a convincing maximum power at the forcing scale for the shallow case at all heights, with larger scales having power comparable to the small scales. In this case the maximum power in the vertical velocity field occurred near mid-height. In the deep case there was instead a maximum power at the box scale at mid-depth and upper levels. The forcing scale in this case appeared as a peak in the wavenumber spectrum only at small heights, where the power density was also a maximum.

Figure 8.

Horizontal wavenumber spectra for (a) shallow (exp. 27) and (b) deep (exp. 30) cases, expressed as the normalised power spectral density of mean kinetic energy. These were found at mid-depth (black line) and at an upper level (grey line); ( $z/L=1.5$ in (a) and $z/L=3.75$ in (b)). The dashed lines correspond to the width of the box ( $k=1/4L$ ) and the forcing scale ( $k=1/L$ ). Dominant peaks are at the forcing scale for the shallow case and at the box scale for the deep case.

3.4 Frequency spectra

In order to investigate the small scales, frequency spectra were obtained from in situ conductivity fluctuations measured at fixed locations. The conductivity was logged at 10 Hz and the time series spanned the duration of an experiment in cases where only one location was sampled. The record was cut into 1 hr segments, the frequency spectra computed for each segment, and then averaged. In other experiments several records of 30–60 min were obtained at different locations. Spectra for shallow and deep cases are shown in figure 9. Both cases have a maximum power density near frequency $U_{0}/(2\unicode[STIX]{x03C0}H)$ (where $U_{0}$ is the measured speed of the box scale of motion and $H/U_{0}$ is approximately the time for advection of tracer or flow structures through the depth  $H$ ). For the deep case (figure 9 b) the low frequencies, corresponding to fluctuations of the large scales, include the observed translations, rotations and reversals of large overturning cells, whereas in the shallow case (figure 9 a) the low frequencies represent only relatively small variability of the array of discrete plumes about a mean state. Higher frequencies, on the other hand, correspond largely to advection of small-scale structure past the measuring probe. For the shallow cases there is only a narrow range (frequencies greater than approximately $H/(2\unicode[STIX]{x03C0}Lt_{b})$ , corresponding to scales a little smaller than  $L$ ) in which the spectral slope is not inconsistent with the $-5/3$ expected for an inertial subrange. This is because experiments achieve only modest Reynolds numbers for the discrete plumes. In contrast, for the deep cases the spectrum is in close agreement with a $-5/3$ slope across more than two orders of magnitude in frequency, corresponding approximately to length scales ranging from somewhat smaller than the forcing scale $L$ up to the domain scale $H$ or  $W$ .

Figure 9.

Measured frequency spectra of density fluctuations in experiments, derived from conductivity for (a) a shallow case ( $A=2$ , at height $z=L$ , exp. 28) and (b) a deep case ( $A=5$ , at $z=4L$ , exp. 31). In (a) the spectrum is from a single long record and in (b) it is the average of four 25 min records from different locations in the domain. Spectra have been smoothed by 3-point binning. The vertical lines show $H/(2\unicode[STIX]{x03C0}Lt_{b})$ (the higher frequency of the two) and $U_{0}/(2\unicode[STIX]{x03C0}H)$ (where $U_{0}$ is the measured velocity of the largest scale of motion), which serve to indicate approximate frequencies corresponding to the advection of features having scale  $L$ and  $H$ , respectively. The power law slope $-5/3$ is shown in (a) for comparison at the frequencies expected of turbulence in discrete plumes; a power law regression fit to frequencies between that of box-scale overturning (vertical line) and $f=1~\text{s}^{-1}$ is shown in (b), and gives power  $-1.71$ .

4.2 DNS results

Simulations of the shallow case ( $A=0.5$ ) with $n=3$ and $n=4$ show a coherent plume ascending above each of the hotspots on the bottom boundary, even at the largest $Ra$ , at which each plume is turbulent (a solution for this case with $Ra=10^{12}$ is shown in figure 1 of Griffiths & Gayen (Reference Griffiths and Gayen2015) and is not reproduced here). The plumes become turbulent at around $Ra\geqslant 10^{9}$ and expand with height as a result of turbulent entrainment of the surrounding fluid. These results are consistent with the saline experiments. For a deep box ( $A=5$ ) and $n=3$ the coherent plumes continue to persist in a stable array, again consistent with observations in the analogue experiments.

Figure 10.

Snapshots of the DNS solutions for $n=4$ and deep water ( $A=5$ ) for (a)  $Ra=10^{7}$ , (b)  $Ra=10^{8}$ , (c)  $Ra=10^{9}$ and (d)  $Ra=10^{12}$ . The three-dimensional renditions show the vertical component of velocity, $w$ , in locations where $w$ is positive (black where $w\leqslant 0$ ). Colour scale indicates velocity normalised by $(g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH)^{1/2}$ . See figure 1 of Griffiths & Gayen (Reference Griffiths and Gayen2015) for the shallow-water case corresponding to (d).

Figure 11.

Kinetic energy power spectra for the deep case ( $A=5$ ), with $n=4$ at four values of $Ra=10^{7}{-}10^{12}$ . The KE is obtained from the three-dimensional DNS velocity field, with velocity normalised by $(g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH)^{1/2}$ , and plotted as a function of horizontal wavenumber $k$ (with wavenumber normalised by the forcing wavenumber $2\unicode[STIX]{x03C0}/L$ of the forcing). Data are averages over 20 distributed locations in the interior of the flow. Vertical dotted lines show the width of the domain at $2\unicode[STIX]{x03C0}/L=0.25$ . A $-5/3$ slope (broken line) describes the spectrum for $Ra=10^{12}$ over the wavenumber range $0.25\leqslant kL/2\unicode[STIX]{x03C0}\leqslant 30$ .

In contrast, the deep case with $n=4$ shows a flow field that changes dramatically with Rayleigh number. At $Ra=10^{7}$ the velocity field (figure 10 a, here normalised by the inertial-buoyancy velocity scale $\sqrt{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH}$ ) shows an array of discrete, almost steady plumes in which the fluid ascends vertically above each hot patch. The plumes remain stable in a square array, although the upper portion of the columnar structures shows a small unsteadiness, the first signs of plume–plume interactions. At this $Ra$ the flow is completely laminar, as is also revealed from the kinetic energy spectrum (figure 11), which shows a sharp decay of energy on each side of a peak at the forcing length scale. At $Ra\sim 10^{8}$ (figure 10 b) the plumes interact more obviously with each other and, as a result, there is significant energy at length scales up to the domain width, with a maximum in the energy spectrum (figure 11) at $4L$ rather than at the forcing scale. This implies an upscale energy transfer. The emergent large-scale circulation tends to organise the smaller scales, acts as a mean ‘wind’ and sweeps plumes into quasi-two-dimensional structures. At $n=4$ and $A=5$ the flow at $Ra=10^{8}$ is unsteady, but there is not a well-developed turbulent energy cascade to small scales.

At $Ra=10^{9}$ (figure 10 c) there is again a spectral peak at the box scale, but more small-scale turbulence kinetic energy is evident. The spectrum (figure 11) is consistent with a $-5/3$ power law expected for the Kolmogorov inertial turbulence subrange, over one order of magnitude in wavenumber, along with a steeper dissipation range at small scales extending to the Kolmogorov scale. At $Ra=10^{12}$ the flow is little further changed apart from evidence of still smaller scales and a significant broadening of the inertial turbulence subrange, which now extends over two orders of magnitude in wavenumber. This change in the smallest scales accurately reflects the predicted dependence of the average Kolmogorov scale on Rayleigh number ( $\unicode[STIX]{x1D702}_{K}\sim L(H/L)^{1/4}Pr^{1/2}Ra^{-1/4}$ ). There is greater intermittency in the warm plumes near the bottom boundary; the large scales involve sweeping motions that tend to carry the warm fluid away from the hotspots of the base, collecting it into one or two updraughts. In contrast to flow in the Rossby problem, where vertical motions are highly asymmetric, the present arrangement shows comparable dimensions and speeds for the large updraughts and downdraughts. This symmetry is also apparent in the saline laboratory experiments reported here.

Figure 12.

(a) Dimensionless heat transport $Nu$ as a function of $Ra$ from DNS for the shallow-water (red circles) and deep-water (green squares) cases. Lines show $Nu=0.195\,Ra^{1/5}$ (red dashed) and $Nu=0.08\,Ra^{1/4}$ (green dashed). The red line is extrapolated to larger values of $Ra$ only as a reference and does not represent a single asymptotic similarity regime. In (b) the heat fluxes from (a) are renormalised as $Nu/Ra^{1/5}$ and lines corresponding to those in (a) are replotted.

The normalised heat flux (Nusselt number) was estimated from

where $\unicode[STIX]{x1D70C}_{0}$ is a reference density, $c_{p}$ is the specific heat, the (dimensional) heat input is $Q=\unicode[STIX]{x1D70C}_{0}c_{p}\unicode[STIX]{x1D705}\int _{a_{+}}\overline{\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}z}|_{z=0}\,\text{d}a$ , the overbar indicates time average and $a_{+}$ is the area over which the vertical temperature gradient on the base is positive (i.e. the cooling area). The results (figure 12) indicate that there are at least two, and arguably three, heat transport regimes. At $Ra\leqslant 10^{8}$ the flow is laminar and the Nusselt numbers for both shallow and deep cases are consistent with $Nu\sim 0.195\,Ra^{1/5}$ , which is the Rossby boundary layer scaling for steady, viscous convection (found at $Ra\leqslant 10^{9}$ in the Rossby problem).

For $Ra\geqslant 10^{9}$ the results depend on water depth. In the shallow case there is a small increase in the heat transport, consistent with a gradual transition, as found in the Rossby problem at $10^{10}\leqslant Ra\leqslant 10^{11}$ (Gayen et al. Reference Gayen, Griffiths and Hughes2014), towards a turbulent inertial regime in which flow becomes governed by turbulent Reynolds stress in the boundary layer. (We do not examine here the effects of radial flow of the boundary layer into each ascending plume as compared with the rectangular geometry in the Rossby problem.) This regime again yields the $1/5$ scaling (but with a larger pre-factor than the viscous regime).

For the deep case the heat transport at $Ra\geqslant 10^{9}$ (figure 12) departs markedly from the $Ra^{1/5}$ scaling. The results are described very well by $Nu\sim 0.08\,Ra^{1/4}$ across four orders of magnitude in $Ra$ . The pre-factor is expected to depend on $n$ , $A$ and $Pr$ . The occurrence of the $1/4$ scaling is coincident with conditions under which scales of motion larger than the forcing and comparable to the width of the box emerge from the turbulence. The larger scales can more rapidly release available potential energy to motion, and transport heat more rapidly, than do the smaller scales. The transition between regimes for $n=4$ and $A=5$ lies at $Ra\approx 10^{8}$ , where we observe the beginning of strong interactions between the discrete plumes (figure 10 b,c).

Closely related to the heat transport is the generation of available potential energy by the boundary heating and cooling, at the (time-averaged) rate $\unicode[STIX]{x1D6F7}_{b2}=\unicode[STIX]{x1D705}g\int _{S}\overline{z^{\ast }\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}z}|_{z=0}\,\text{d}S$ , where $S$ is the surface area of the boundary and $z^{\ast }$ is the height at which each parcel of fluid would lie if the fluid were reorganised adiabatically into the state of gravitational equilibrium (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Hughes et al. Reference Hughes, Hogg and Griffiths2009). A suitable normalisation of this (and other terms in the mechanical energy budget, as discussed below) is a comparison with the rate of viscous dissipation of kinetic energy, the time-averaged box total of which is given by $\unicode[STIX]{x1D700}_{T}=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D708}\int _{V}\overline{\Vert \unicode[STIX]{x1D735}\boldsymbol{u}\Vert ^{2}}\,\text{d}V$ , where $V$ is the total volume of the convecting fluid. The APE generation rate normalised by $\unicode[STIX]{x1D700}_{T}$ , shown in figure 13(a), increases with both $Ra$ and the depth  $H$ . APE generation reaches approximately 40 (shallow case) and 600 times (deep case) the viscous dissipation.

Figure 13.

(a) The rates of generation of available potential energy, $\unicode[STIX]{x1D6F7}_{b2}$ , by boundary buoyancy forcing, normalised by the total viscous dissipation rate integrated over the domain volume, as functions of $Ra$ from DNS for the shallow-water (red circles) and deep-water (green squares) cases. Lines show the power laws $\unicode[STIX]{x1D6F7}_{b2}/\unicode[STIX]{x1D716}_{T}\sim Ra^{1/5}$ (red) and $Ra^{1/4}$ (green). In (b) the values from (a) are renormalised by $ARa^{1/5}$ (where $A=H/L$ ) and lines corresponding to those in (a) are replotted.

The trends in APE generation are better represented by the renormalised form shown in figure 12(b), which is based on the following scaling. Considering the hypothetical re-sorting of fluid elements required to define APE, in which the cold boundary layer fluid stays at the base and the hottest fluid goes to the top, we rewrite the above definition of APE generation by surface buoyancy forcing as $\unicode[STIX]{x1D6F7}_{b2}\approx \unicode[STIX]{x1D705}g\unicode[STIX]{x1D6FC}\langle \overline{\text{d}T/\text{d}z}|_{-}\rangle a_{-}H$ , where the subscript ‘-’ indicates the boundary temperature gradient in the region of negative vertical temperature gradient at the boundary (i.e. the heating area $a_{-}$ ), the overbar represents time average, $\langle \cdot \rangle$ indicates horizontal averaging over the relevant area and we recognise that $|z^{\ast }|\approx H$ . Since heat input is equal to heat withdrawal this can be expressed as $\unicode[STIX]{x1D6F7}_{b2}=\unicode[STIX]{x1D705}g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}T[L\langle \overline{(\text{d}T/\text{d}z)}|_{+}\rangle /\unicode[STIX]{x0394}T](H/L)a_{+}$ . The normalisation of this quantity is based on the viscous dissipation rate in horizontal convection, which is equal to the rate of generation of potential energy $\unicode[STIX]{x1D6F7}_{i}$ ( $\unicode[STIX]{x1D700}_{T}=\unicode[STIX]{x1D6F7}_{i}$ ), which is bounded by $\unicode[STIX]{x1D6F7}_{i}\approx -g\unicode[STIX]{x1D705}(\langle \unicode[STIX]{x1D70C}\rangle _{H}-\langle \unicode[STIX]{x1D70C}\rangle _{0})a_{T}$ , where $a_{T}$ the horizontal area of the domain and horizontal averaging of density is at the top ( $z=H$ ) and bottom ( $z=0$ ) of the domain (Paparella & Young Reference Paparella and Young2002; Hughes et al. Reference Hughes, Hogg and Griffiths2009). The computed total dissipation rate in DNS is not measurably different from this upper bound (Gayen et al. Reference Gayen, Griffiths and Hughes2014; Griffiths & Gayen Reference Griffiths and Gayen2015) and we expect the same expression to hold in the case of the present array forcing. (It was reported in Griffiths & Gayen (Reference Griffiths and Gayen2015) that periodic side boundary conditions gave $\unicode[STIX]{x1D700}_{T}$ larger than $\unicode[STIX]{x1D6F7}_{i}$ . However, this was due to an error in the normalisation of the integrated dissipation rate. The computed total dissipation is actually indistinguishable from $\unicode[STIX]{x1D6F7}_{i}$ for both types of boundary conditions.) Working from $\unicode[STIX]{x1D6F7}_{i}$ , it can be shown that the total viscous dissipation rate is $\unicode[STIX]{x1D700}_{T}\approx g\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}Ta_{+}$ , or equivalently $\unicode[STIX]{x1D700}_{T}\sim (\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D705}^{3}a_{+}/L^{3})(Ra\,Pr)$ . The normalised APE generation rate becomes

which is identical to the previous theoretical scaling proposed for the rectangular Rossby problem (Gayen et al. Reference Gayen, Griffiths, Hughes and Saenz2013a ) and found to be a robust result (Gayen et al. Reference Gayen, Griffiths and Hughes2014; Vreugdenhil et al. Reference Vreugdenhil, Gayen and Griffiths2016). In comparison to the heat transport $Nu$ , APE generation has an additional dependence on water depth $H$ (or the aspect ratio  $A$ ), as a result of gravitational potential. Thus the DNS results are renormalised in figure 12(b) as $\unicode[STIX]{x1D6F7}_{b2}/(\unicode[STIX]{x1D700}_{T}ARa^{1/5})$ and the results closely mirror the trends in renormalised Nusselt number, $Nu/Ra^{1/5}$ (figure 12 b).

The flow velocities are represented in figure 14 as the Reynolds number $Re=UL/\unicode[STIX]{x1D708}$ , which is evaluated from the root-mean-square magnitude of velocities within the volume of the thermal boundary layer (defined, as before, by $z<\unicode[STIX]{x1D6FF}$ , where $\unicode[STIX]{x1D6FF}$ encompasses 95 % of the horizontally averaged bottom-to-top temperature difference). Again these velocities show at least two regimes: the deep-water case at $Ra>10^{8}$ diverges markedly from all other cases and is characterised by larger velocities.

Figure 14.

Reynolds number $Re=UL/\unicode[STIX]{x1D708}$ (a) based on the boundary layer velocity computed from the kinetic energy, $U=\langle (u^{2}+v^{2}+w^{2})/2\rangle ^{1/2}$ , and averaged over the volume of the boundary layer, shown as a function of Rayleigh number, from DNS for the shallow-water (red circles) and deep-water (green squares) cases. Lines show the power laws $Re\sim Ra^{2/5}$ (red) and $Re\sim Ra^{1/2}$ (green). In (b) the values from (a) are normalised by $Ra^{2/5}$ and lines corresponding to those in (a) are replotted.

Figure 15.

Volume-integrated shear production rates $P$ (conversion of mean to turbulent kinetic energy) as functions of Rayleigh number from DNS for the (a) shallow-water and (b) deep-water cases. TKE production rate is integrated over the volume of the thermal boundary layer (open symbols) or the whole domain (filled symbols), and normalised by the total domain-integrated viscous dissipation rate.

The Rayleigh number dependence of the shear production $P=-\unicode[STIX]{x1D70C}_{0}\sum _{i,j=1}^{3}\overline{u_{i}^{\prime }u_{j}^{\prime }}(\unicode[STIX]{x2202}\overline{u}_{i}/\unicode[STIX]{x2202}x_{j})$ (where overbars again represent time averages), when integrated over both the full volume of the domain, $P_{T}=\int _{V}P\,\text{d}V$ , is shown in figure 15, where it is given as a fraction of the total viscous dissipation rate $\unicode[STIX]{x1D716}_{T}$ . For comparison, the proportion of $P_{T}$ found in the boundary layer, $P_{bl}$ , is also shown. Shear production is insignificant in the viscous regime at $Ra<10^{8}$ . Shear production grows to 30 % (shallow case) or 90 % (deep case) of the total kinetic energy production (and dissipation) rate at $Ra=10^{12}$ . At this Rayleigh number the interior accounts for most of the shear production in both cases (approximately 83 % in the shallow case and 89 % in the deep case). Thus the boundary layer accounts for a small fraction of the shear production (and viscous dissipation) under all conditions. The bottom boundary layer at $Ra=10^{12}$ also has a logarithmic profile of the horizontal velocity (Griffiths & Gayen Reference Griffiths and Gayen2015), consistent with shear production within the boundary layer. However, the very small volume of the boundary layer leaves the viscous dissipation almost entirely in the interior.

Figure 16.

Profiles of irreversible mixing rate $\unicode[STIX]{x1D6F7}_{d}$  (a,b) and viscous dissipation rate (c,d) with height for $Ra=10^{9}$ (green line) and $10^{12}$ (red line) from DNS of the (a,c) shallow-water ( $H/L=0.5$ ) and (b,d) deep-water ( $H/L=5$ ) cases. The computed energy dissipations per unit volume are averaged in time and integrated over the horizontal area $W^{2}$ of the domain at each depth, then normalised by the total rate of generation of potential energy $\unicode[STIX]{x1D6F7}_{i}$ divided by the domain depth  $H$ . Horizontal lines show the edge of the thermal boundary layer ( $\unicode[STIX]{x1D6FF}/H\approx 0.002$ for the deep case at $Ra=10^{12}$ ). The boundary layer $\unicode[STIX]{x1D6FF}$ is defined as the layer containing 95 % of the mean temperature difference: at $z=\unicode[STIX]{x1D6FF}$ the horizontally averaged temperature $\langle T\rangle =\langle T\rangle _{H}+0.05(\langle T\rangle _{H}-\langle T\rangle _{0})$ .

In contrast with the viscous dissipation is the behaviour of the other sink in the mechanical energy equations: the irreversible mixing $\unicode[STIX]{x1D6F7}_{d}=-g\unicode[STIX]{x1D705}\int \overline{(\text{d}z^{\ast }/\text{d}\unicode[STIX]{x1D70C})|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}|^{2}}\,\text{d}V$ (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Hughes et al. Reference Hughes, Hogg and Griffiths2009). Figure 16 shows the mixing and viscous dissipation rates as depth profiles of the horizontally integrated $\unicode[STIX]{x1D6F7}_{d}$ and $\unicode[STIX]{x1D700}$ normalised by $\unicode[STIX]{x1D6F7}_{i}/H$ ( $=\unicode[STIX]{x1D700}_{T}/H$ ) for two values of the Rayleigh number. For the shallow case (figure 16 a,c) there is not much difference between $Ra=10^{9}$ and $Ra=10^{12}$ in either the overall (normalised) mixing rate or its distribution. The mixing rates are of $O(1)$ in the interior and $O(10^{4})$ in the thermal boundary layer. Thus for the shallow cases the mixing rates in the interior are similar to the viscous dissipation rates, whereas mixing in the boundary layer is much larger than both the average and local dissipation rates. Most of the mixing is within the volume of the boundary layer (in line with previous results for the Rossby problem, Gayen et al. Reference Gayen, Griffiths and Hughes2014). For the deep case (figure 16 b,d), on the other hand, the normalised mixing rates depend more strongly on Rayleigh number. For $Ra=10^{9}$ mixing is not much greater than in the shallow case, but at $Ra=10^{12}$ the mixing rates are two orders of magnitude greater. The mixing rate is again $10^{3}{-}10^{4}$ times greater in the boundary layer than in the interior, giving a preponderance of mixing in the boundary layer volume, despite the thinness of the boundary layer. Vertical integrations show that the confinement of mixing to the boundary layer is stronger for the deep case with larger $Ra$ . Also in this case, integration over the whole domain gives a mixing rate approximately 600 times the integrated viscous dissipation. The DNS confirms that the APE budget is dominated by an exact balance,

between the rate of energy supply by the surface buoyancy fluxes and its conversion to background potential energy by irreversible mixing, as predicted by Hughes et al. (Reference Hughes, Hogg and Griffiths2009) and Tailleux (Reference Tailleux2009).