Dynamics of mixed convective–stably-stratified fluids

We study the dynamical regimes of a density-stratified fluid confined between isothermal no-slip top and bottom boundaries (at temperatures $T_t$ and $T_b)$ via direct numerical simulation. The thermal expansion coefficient of the fluid is temperature dependent and chosen such that the fluid density is maximum at the inversion temperature $T_b>T_i>T_t$. Thus, the lower layer of the fluid is convectively unstable while the upper layer is stably stratified. We show that the characteristics of the convection change significantly depending on the degree of stratification of the stable layer. For strong stable stratification, the convection zone coincides with the fraction of the fluid that is convectively unstable (i.e., where $T>T_i)$, and convective motions consist of rising and sinking plumes of large density anomaly, as is the case in canonical Rayleigh-Bénard convection; internal gravity waves are generated by turbulent fluctuations in the convective layer and propagate in the upper layer. For weak stable stratification, we demonstrate that a large fraction of the stable fluid (i.e., with temperature $T<T_i)$ is instead destabilized and entrained by buoyant plumes emitted from the bottom boundary. The convection thus mixes cold patches of low density-anomaly fluid with hot upward plumes and the end result is that the $T_i$ isotherm sinks within the bottom boundary layer and that the convection is entrainment dominated. We provide a phenomenological description of the transition between the regimes of plume-dominated and entrainment-dominated convection through analysis of the differences in the heat transfer mechanisms, kinetic energy density spectra, and probability density functions for different stratification strengths. Importantly, we find that the effect of the stable layer on the convection decreases only weakly with increasing stratification strength, meaning that the dynamics of the stable layer and convection should be studied self-consistently in a wide range of applications.

Figure 1

Equation of state for the density anomaly $\rho$ as a function of temperature $T$ $(T$ decreasing upward) for three different stiffness parameters $\mathcal{S}$ (solid lines). For comparison with water, the quadratic equation of state $\rho =-T^2$ is also shown (dashed line); note that the $\mathcal{S}$ parameter that best represents water depends on the temperature range considered.

Figure 2

Transient evolution of (a) and (b) the temperature field and (c) and (d) the normalized vorticity field $\omega {\tau }_c$ at times (a) and (c) $t/{\tau }_c=2.1$ and (b) and (d) $t/{\tau }_c=50.5$. The simulation is started from the static equilibrium of the linear temperature profile and the physical parameters are $\text{Pr}=1,\text{Ra}=8\times {10}^7,\mathcal{S}=1$, and $T_t=-20$. Here ${\tau }_c$ is the reference turnover time scale (see the text).

Figure 3

Simulation snapshots for $\text{Ra}=8\times {10}^7$ and stiffness (a) $\mathcal{S}=2^8$, (b) $\mathcal{S}=1$, and (c) $\mathcal{S}=2^{-8}$ representatives of the system's dynamic at thermal equilibrium, i.e., obtained after several thousand convective turnover times. The density anomaly is shown in the convection zone, below the neutral buoyancy height $Z_{\text{NB}}\left(x\right)$ [i.e., where $\rho (x,Z_{\text{NB}})=-1]$. Above $Z_{\text{NB}}$, the fluid is stably stratified and we plot the vorticity. In (a) convective cells with a top-down symmetry are shown in the lower layer, generating global wave modes in the upper stably stratified layer. In (b) the convection zone is larger, but significant coupling between the two layers is obtained as shown by the strong variations of $Z_{\text{NB}}$ with $x$. In (c) the system changes significantly with rising plumes close to the bottom interacting with large-scale structures (shown by contour lines) higher in the fluid (note that $Z_{\text{NB}}$ is outside the domain such that we show the density field everywhere in this case). We also show the mean temperature and anomaly density profiles by the dashed and solid lines as later discussed and reported in Fig. 6 (see also [30]).

Figure 4

Snapshot of the temperature field corresponding to the density anomaly plot of Fig. 3.

Figure 5

(a) Bulk temperature $T_{\mathrm{bulk}}$ as a function of $\mathcal{S}$ $(\mathcal{S}$ decreases to the right) averaged over the convection zone (see the text for details). The dash-dotted line is a best fit obtained for all data points and is given by Eq. (6). The dashed line highlights the inversion temperature $T_i=0$ where $\alpha$ changes sign and the dotted line represents the bulk temperature expected from classical RBC. (b) Total heat flux $Q$ as a function of $\mathcal{S}$. Results are shown for three different reference Rayleigh numbers $\text{Ra}$.

Figure 6

(a) Mean temperature ${\overline{\langle T\rangle }}_x$ and (b) mean density anomaly ${\overline{\langle \rho \rangle }}_x$ profiles as functions of $z$ for $\text{Ra}=8\times {10}^7$ and for different stiffness parameters $\mathcal{S}$. For strong stratification, the density anomaly in the convection zone is constant within a bulk bounded by two thin layers where the density quickly increases, reminiscent of classical RBC. For weak stratification, the fluid density only increases within the bottom boundary layer and is almost uniform or weakly decreasing above. As in Fig. 5, the vertical dashed line in (a) highlights the inversion temperature $T_i=0$ and the dotted line indicates the bulk temperature expected from classical RBC. The open squares show results from a convection-only simulation with rigid no-slip top and bottom boundaries and a vertical height $h=0.32$, which is an estimate of the depth of the convection zone for $\mathcal{S}=2^8$ [obtained from $h=2L_{db}+L_{cp}+L_{ce}$; see (9)].

Figure 7

Vertical profiles of the heat transfer due to (a) convective plumes $q_{cp}$, (b) entrainment $q_{ce}$, and (c) diffusion $q_d$ for $\text{Ra}=8\times {10}^7$. Different curves correspond to different stiffness $\mathcal{S}$. For high stiffness, we obtain a symmetry in the heat transfer in the convection zone due to plumes $q_{cp}$ in (a), suggesting classical RBC. The contribution from entrainment to the total heat transfer can be seen in (b) to be dominant in the low-stiffness cases and negligible in the high-stiffness regime. The total heat transfer $q=q_{cp}+q_{ce}+q_d$ is constant along the $z$ axis up to 2% maximum relative discrepancy. As in Fig. 6, the open squares show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case.

Figure 8

Volume fractions occupied by the bottom thermal boundary layer, the plume zone, the entrainment zone, and the stably stratified top layer as functions of the stiffness parameter $\mathcal{S}$ and for $\text{Ra}=8\times {10}^7$ (see Fig. 15 in the Appendix for $\text{Ra}=8\times {10}^8)$. The layer thicknesses are computed based on the profiles shown in Fig. 7 and according to Eqs. (9). For each stiffness parameter, the position where the average temperature in the simulations is equal to the inversion temperature is shown by a star. It can be clearly seen that the inversion temperature isotherm is on average close to or within the bottom thermal boundary layer for $\mathcal{S}\le 2^{-2}$. The double-headed arrow on the left-hand side shows the convection layer thickness $h$ prediction based on Eq. (4). The closed diamonds show the empirical prediction of $h$ based on Eq. (14) presented in Sec. 5f.

Figure 9

Spectrogram of the kinetic energy density spectrum $0.5f{\langle {\widehat{u}}^2+{\widehat{w}}^2\rangle }_x$ as a function of frequency and $z$ for each of the three simulations shown in Fig. 3 $(\text{Ra}=8\times {10}^7$; $\mathcal{S}=2^8,2^0,2^{-8}$ from left to right). The frequencies are normalized by the reference turnover frequency $f_c=\sqrt{\text{Ra}\text{Pr}}/2\pi$. The black solid (white dashed) line represents the mean value of the real (imaginary) part of the normalized buoyancy frequency, i.e., $N/f_c$ with $N={\overline{\langle \sqrt{-\text{Ra}\text{Pr}{\partial }_z\rho }\rangle }}_x/2\pi$ based on our notation. The arced dashed lines in (a) highlight global internal-wave modes.

Figure 10

Probability density function of temperature $\mathcal{P}\left(T\right)$ for different stiffness parameters $\mathcal{S}$ and $\text{Ra}=8\times {10}^7$. As $\mathcal{S}$ is increased, the peak of $\mathcal{P}\left(T\right)$ tends to 0.5 (dotted line), the bulk temperature expected for classical RBC. As $\mathcal{S}$ is decreased, the peak of $\mathcal{P}\left(T\right)$ drops below the inversion temperature (dashed line). The profile of $\mathcal{P}\left(T\right)$ is symmetric about the peak for hight stiffness (the black solid line shows exponential fit), but asymmetric for low stiffness. Small local maxima of $\mathcal{P}\left(T\right)$ in the range $T\in [0.5,1]$ are expected to be representative of the plumes' temperature. As in Fig. 6, the open squares show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case. Note that $\mathcal{P}\left(T\right)$ for the convection-only results is normalized such that the integral value is equal to the integral value of $\mathcal{P}\left(T\right)$ for $\mathcal{S}=2^8$ over the range $T\in [0,1]$.

Figure 11

(a) Probability density function $\mathcal{P}$ of velocity $w$ from the convective region for different stiffness parameters $\mathcal{S}$ and for $\text{Ra}=8\times {10}^7$. Note that the convective region encompasses both the entrainment zone and the plume zone (cf. Fig. 8). (b) Same as (a) but normalized by the standard deviation $\sigma$ and mean $\mu$ of the distributions. For high stiffness (red solid line), the distribution looks exponential for slow speeds (black dashed line). For low stiffness (blue solid line) the distribution looks Gaussian (black dotted line) for slow speeds with asymmetric exponential tails. As in Fig. 6, the open squares show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case.

Figure 12

Joint probability density function of temperature and vertical velocity for $\text{Ra}=8\times {10}^7$ and (a) $\mathcal{S}=2^8$, (b) $\mathcal{S}=2^0$, (c) $\mathcal{S}=2^{-8}$. Note that the oscillations in (c) are most likely due to statistics that are not yet fully converged in the stable layer for the highest-stiffness case.

Figure 13

Variations of $\text{Nu}/{\text{Ra}}^{0.28}$ with $\mathcal{S}$ where $\text{Nu}=-Q/(T_t-1)$ for three reference Rayleigh numbers. The dashed line shows the best-fit power-law regression of $\text{Nu}$ in terms of $\mathcal{S}$ for $\text{Ra}=8\times {10}^7,8\times {10}^8$.

Figure 14

Rayleigh-Nusselt scaling based on ${\text{Nu}}_C=Qh$ and effective Rayleigh number ${\text{Ra}}_{\mathrm{eff}}=\text{Ra}{h}^3$, where $h$ is an estimate of the convective layer depth [see Eq. (12) and the text for details]. The circles, triangles, and stars correspond to reference Rayleigh number $\text{Ra}=8\times {10}^6,8\times {10}^7,8\times {10}^8$, respectively. The dash-dotted line has slope 0.28, while the dashed line has slope 1. Results of five convection-only simulations are shown by closed squares; note that the lowest ${\text{Ra}}_{\mathrm{eff}}$ convection-only simulation is stationary in time; also, the scaling exponent for the convection-only results is 0.28 (discarding the lowest ${\text{Ra}}_{\mathrm{eff}}$ case).

Figure 15

Same as Fig. 8 but for the higher-Rayleigh-number $\text{Ra}=8\times {10}^8$ case.