3.1. Flow organization

In this subsection, a qualitative description of the main features of the two-layer Rayleigh–Bénard convection is presented through a series of different flow visualizations. First, figure 3 presents snapshots of the temperature volumetric rendering, including the isosurface representing the location of the two-fluid interface. The activity of the thermal plumes is depicted clearly, especially at the top half of the cavity, where hotter plumes rise from the interface along with colder plumes descending from the top cooled wall. The plume activity typically feeds the large-scale circulation structures that sweep the boundary layers at their periphery. The bottom half of the cavity appears to be less active, with thicker thermal structures that occupy a significant portion of the available volume. As observed in single-phase thermal convection, the thermal structures become finer as the Rayleigh number increases, intensifying the thermal transport from the boundary layers to the large-scale circulation. On the contrary, based on the various snapshots in figure 3, the Weber number does not appear to have any obvious effect on the temperature field.

Figure 3.

Snapshots of volume representations of the temperature field (semi-transparent coloured isosurfaces) and the interface (light brown surface at approximately mid-height) for the six different DNS cases. Red denotes a higher temperature compared to blue.

Focusing on the two-fluid interface as depicted in figure 3, its deformation is mostly visible for the ${Ra}=10^6$ cases. For larger Rayleigh numbers, this deformation becomes less pronounced. More specifically, the r.m.s. values of the interface elevation $\eta _{rms}$ reduce by 69 % from ${Ra}=10^6$ to ${Ra}=10^8$ considering Weber number 100, and by 66 % considering Weber number 1000. The effects of the Weber number on the interface deformation are more clearly depicted in figure 4, where cases with larger Weber numbers exhibit a more noticeable deformation. As concerns the velocity fields, we note that the appearance of progressively smaller structures as the Rayleigh number increases well correlates with the ejected thermal plumes. This observation is more evident in the top half of the cavity, as the bottom half exhibits smaller temperature differences. The more qualitative observations discussed in this subsection will be verified in the quantitative statistical analysis presented in the following subsections.

Figure 4.

Snapshots of $x$–$z$ planar representations of the temperature field, including streamlines (black lines with arrows) and the interface (white dashed line) for the six different DNS cases. Red denotes a higher temperature compared to blue.

3.2. Global properties

The analysis of global properties focuses on the interface temperature ($\varTheta _\varGamma$) and the Nusselt and Reynolds numbers, which represent, in non-dimensional form, heat transfer and turbulence intensity, respectively. In this subsection, results extracted from the present DNS are compared against recently developed scaling laws with respect to the Rayleigh and Weber numbers. We note here that Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022) have argued that the proposed scaling laws can hold for different Rayleigh numbers. Nonetheless, this hypothesis has not been assessed so far. The scaling laws for $\varTheta _\varGamma$, $Nu$ and $Re$ are based on the Grossmann–Lohse (GL) theory, which consists of two nonlinear equations, originating from the relations for the kinetic and thermal energy dissipation rates. For a multiphase problem, this set of equations is written for both layers separately (Liu et al. Reference Liu, Chong, Yang, Verzicco and Lohse2022):

(3.1a)\begin{equation} ({Nu}_j-1)\,{Ra}_j\,{Pr}_j^{{-}2} = \dfrac{c_1\,Re_j^2}{g\left(\sqrt{\dfrac{Re_c}{Re_j}}\,\right)} + c_2\,Re_j^3, \end{equation}

(3.1b) \begin{align} ({Nu}_j-1) &= c_3\,Re_j^{1/2}\,Pr_j^{1/2}\sqrt{f\left[\dfrac{2c_0\,{Nu}_j}{\sqrt{Re_j}}\,g \left(\sqrt{\dfrac{Re_c}{Re_j}}\,\right)\right]} \nonumber\\ &\quad +c_4\,Pr_j\,Re_j\,f\left[\dfrac{2c_0\,{Nu}_j}{\sqrt{Re_j}}\,g\left(\sqrt{\dfrac{Re_c}{Re_j}}\,\right)\right], \end{align}

where indices $j=1,2$ denote the corresponding parameters in the heavier and lighter fluids, respectively. The coefficients $c_{i}$, with $i=0,\ldots,4$, are the prefactors that are chosen as $c_{i}=0.922,8.05,1.38,0.487,0.0252$, as suggested in Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013), together with $Re_c=(2c_0)^2$. The crossover functions $g$ and $f$ in (3.1) are given by $g(x)=(1+x^n)^{1/n}$ and $f(x)=x^n(1+x^n)^{1/n}$ with $n=4$ (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013). The Prandtl number in (3.1) is computed as ${Pr}=\hat {\mu }_j\hat {c}_{p,j}/\hat {\zeta }_j$, while the Rayleigh numbers of each phase, ${Ra}_j$, read

(3.2a)\begin{gather} {Ra}_1 = C_{tot}^3(1/2-\varTheta_\varGamma)\,{Ra}, \end{gather}

(3.2b)\begin{gather}{Ra}_2=(1-C_{tot})^3(1/2+\varTheta_\varGamma)\,{Ra}\,\varLambda_\mu\varLambda_\zeta/ (\varLambda_\rho^2\varLambda_{c_p}\varLambda_\alpha), \end{gather}

where $C_{tot}$ is the overall volume fraction of phase 1. Note that in Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022), the properties ratios $\varLambda _\mu$, $\varLambda _\rho$, $\varLambda _{c_p}$ and $\varLambda _\alpha$ are omitted in (3.2) since they are considered equal to $1$, but are kept here for completeness. Therefore, the Nusselt numbers are defined as

(3.3a)\begin{gather} {Nu}_1=\frac{\hat{Q}_1C_{tot}\hat{l}_z}{\hat{\zeta}_1(1/2-\varTheta_\varGamma)\,\Delta\hat{\varTheta}}, \end{gather}

(3.3b)\begin{gather}{Nu}_2=\frac{\hat{Q}_2(1-C_{tot})\hat{l}_z}{\hat{\zeta}_2(1/2+\varTheta_\varGamma)\,\Delta\hat{\varTheta}}, \end{gather}

where the heat fluxes $\hat {Q}_{1},\hat {Q}_{2}$ at the bottom and top walls are considered equal once a statistical equilibrium is reached, i.e. $\langle \hat {Q}_1\rangle _t=\langle \hat {Q}_2\rangle _t$. Equations (3.1) are applied to each layer separately, coupled to (3.2) and (3.3). The solution of this system of equations provides the values of $\varTheta _\varGamma$, ${Nu}_j$ and $Re_j$. Since the GL theory is implicit in ${Nu}_j$ and $Re_j$, an iterative procedure is required. Once ${Nu}_j$, $Re_j$ and $\varTheta _\varGamma$ are known, the global Nusselt number can be obtained readily from ${Nu}_j$ at one of the interface sides. Taking as reference the top wall and using the relation for $Nu_2$ in (3.3), the global $Nu$ reads

(3.4)\begin{equation} {Nu} = \dfrac{{\hat{Q}_2}\hat{l}_z}{\hat{\zeta}_2\,\Delta\hat{\varTheta}} = {Nu}_2\,\dfrac{1/2+\varTheta_\varGamma}{1-C_{tot}}. \end{equation}

Similarly, we can define two global Reynolds numbers for the top ($Re_{top}$) and bottom ($Re_{bot}$) halves of the cavity. These quantities can be related to the turbulence in each phase, with $Re_{bot}=Re_1\,(\hat {\nu }_1/\hat {\nu }_{ref})/C_{tot}$ and $Re_{top}=Re_2c(\hat {\nu }_2/\hat {\nu }_{ref})/(1-C_{tot})$. Note that in the present set-up, $\hat {\nu }_{1}=\hat {\nu }_{2}=\hat {\nu }_{ref}$ and $C_{tot}=0.5$, therefore $Re_{bot}=2\,Re_{1}$ and $Re_{top}=2\,Re_{2}$.

The Nusselt and Reynolds numbers predictions are shown in figure 5, exhibiting excellent agreement against the corresponding scaling laws. This confirms the observation from Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022) that the weak dependence on the Weber number is expected to hold as long as the interface does not break up and two separated layers can be clearly identified. Furthermore, the Reynolds number is approximately two times higher in the lighter fluid compared to the denser fluid, for all the Weber and Rayleigh numbers explored. This indicates significantly higher turbulence levels in the top half of the cavity.

Figure 5.

(a) Nusselt number and (b) Reynolds number predictions as functions of the Rayleigh number, for different Weber numbers. The GL theory predictions are calculated from (3.1)–(3.3).

Furthermore, as an alternative to (3.1), which requires an iterative solution, an approximation of $\varTheta _\varGamma$ can be obtained explicitly. This approach has been proposed in Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022) for multiphase thermal convection, and in Scapin et al. (Reference Scapin, Demou and Brandt2023) for evaporating Rayleigh–Bénard convection. The main idea is to employ a simplified scaling of the form ${Nu}_{j}=A_j\,{Ra}_j^{\gamma _j}\,{Pr}_j^{m_j}$ and to assume that the Rayleigh and Prandtl numbers of both phases are sufficiently similar to fall inside the same scaling regime of the GL theory, so that $A_1=A_2=A$, $\gamma _1=\gamma _2=\gamma$ and $m_1=m_2=m$. Then, by using (3.3) to express the ${Nu}_j$ and (3.2) to express the ${Ra}_j$, the following approximate relation for $\varTheta _\varGamma$ emerges:

(3.5)\begin{equation} \varTheta_\varGamma ={-}\frac{1}{2}+\left(1+\left(\frac{1-C_{tot}}{C_{tot}}\right)^{{(1-3\gamma)}/{(1+\gamma)}} \left(\frac{\varLambda_\rho^2\varLambda_{c_p}\varLambda_{\alpha}}{\varLambda_\mu}\right)^{{\gamma}/{(1+\gamma)}} \varLambda_\zeta^{{(1-\gamma)}/{(1+\gamma)}}\right)^{{-}1}. \end{equation}

Note that the Prandtl number dependence is typically omitted since the corresponding scaling exponent becomes $0< m\ll 1$ for ${Pr}>0.5$ (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013). In contrast to the implicit relations of the GL theory, (3.5) depends only on the various property ratios $\varLambda$ and the chosen scaling exponent $\gamma$. Figure 6 shows the interface temperature obtained from the DNS, compared against the predictions of the GL theory (3.1)–(3.3) and the predictions of the simplified scaling (3.5). The results of the simplified scaling assume the typical scaling exponent range $1/4\leq \gamma \leq 1/3$. As observed in figure 6, the DNS results from the low Weber number cases follow the GL theory closely. The higher Weber number cases deviate from the GL theory as the lower surface tension induces a stronger interface deformation, which influences the resulting scaling. Noteworthy is the fact that the simplified scaling provides a reliable estimation of $\varTheta _\varGamma$ for the different values of the Rayleigh number considered in the present study, with a maximum deviation of less than 2 % from the GL theory. Furthermore, despite its apparent simplicity, the simplified scaling directly highlights some important features of the equilibrium interface temperature $\varTheta _\varGamma$. First, we highlight that the weak dependence of $\varTheta _\varGamma$ on the Rayleigh number (absent in (3.5)) is confirmed by the DNS results, since the value of $\varTheta _\varGamma$ decreases by less than 3 % as $Ra$ increases from $10^6$ to $10^8$. Second, (3.5) enables us to quantify the role of each thermophysical property in modulating $\varTheta _\varGamma$, because each property ratio $\varLambda$ has its own scaling exponent. In particular, the density and thermal conductivity ratios have the largest exponents, representing the dominant sources of influence on the interface temperature $\varTheta _\varGamma$.

Figure 6.

Average interface temperature $\varTheta _\varGamma$ as a function of the Rayleigh number, for different Weber numbers. The GL theory predictions are calculated from (3.1)–(3.3), while the predictions of the simplified scaling are obtained from (3.5). The green region represents possible $\varTheta _\varGamma$ values for scalings $1/4<\gamma <1/3$, as per the simplified scaling.

We conclude this subsection by examining the thicknesses of the scaled thermal boundary layers at different locations in the cavity, which are shown in figure 7. The single-phase scaling of ${Ra}^{-1/4}$ is used, as this scaling was proven to agree reasonably well with the interface temperature predictions in figure 6. Due to the presence of two fluid layers, the effective Rayleigh numbers ${Ra}_{1,2}$ defined in (3.2) are used to scale the thermal boundary layer thickness at the bottom and top halves of the cavity. The effective Rayleigh numbers consider the height and the temperature difference of the corresponding fluid layer, and for the parameters of the present study, can be simplified as ${Ra}_1 = {Ra}\,(0.5-\varTheta _\varGamma )/8$ at the bottom layer, and ${Ra}_2 = {Ra}\,(\varTheta _\varGamma +0.5)/8$ at the top layer. As shown in figure 7, this scaling is successful in collapsing the thicknesses of the thermal boundary layers that form on both the solid surfaces and the fluid interface. Note that this is true despite the variation of the Weber number, which, however, does not have a clear impact on the scaled thermal boundary layer thicknesses. As shown in figure 6, the Weber number has a minor influence on the interface temperature, which is not clearly noticeable when examining the thermal boundary layer thicknesses in figure 7. Indeed, the value of the interface temperature, $\varTheta _\varGamma$ appears in the definitions of both $h^{\theta }_{1,2}$ (2.10) and ${Ra}_{1,2}$ (3.2).

Figure 7.

Scaled thermal boundary layer thickness at (a) the solid walls and (b) the interface. The scalings used are ${Ra}_{1}^{-1/4}$ at the bottom and ${Ra}_{2}^{-1/4}$ at the top halves of the cavity, where ${Ra}_{1,2}$ are defined in (3.2).

3.3. Two-phase statistics

The vertical distributions of the average and r.m.s. phase indicator functions are shown in figures 8(a,b), respectively, with a focus on the region adjacent to the interface. As the Rayleigh number increases, the transition from one fluid to the other at $z=0.5$ becomes sharper, as shown by the average distribution of the indicator function. This indicates the reduced interface deformation as the Rayleigh number increases. The same effect can be deduced from the r.m.s. values of the indicator function, where the fluctuations are shown to extend to a progressively smaller region as the Rayleigh number increases. This observation can be linked to the size of the thermal plumes that cause the interface deformation. In fact, with increasing Rayleigh number, the thermal plumes become progressively smaller (Zhou & Xia Reference Zhou and Xia2010), while the flow in the bulk becomes increasingly more homogeneous (Zhou, Sun & Xia Reference Zhou, Sun and Xia2008). Therefore, the stronger interface deformations at low Rayleigh numbers result from the larger thermal plumes in the cavity. Furthermore, increasing the Weber number leads to smoother profiles in the average distribution in figure 8(a), and larger fluctuating regions in figure 8(b). These trends are due to the decreased surface tension, leading to larger interface deformation. The interface deformation will be analysed further and characterized in § 3.4.

Figure 8.

Vertical distributions of (a) average and (b) r.m.s. values of the indicator function for different values of the Rayleigh and Weber numbers.

The vertical profiles of the temperature field statistics are depicted in figure 9. The data in the figure suggest several considerations. First, the average temperature field in figure 9(a) reveals an increasingly sharper temperature profile as the Rayleigh number increases. Specifically, at the bottom and top walls of the cavity, the temperature gradients become larger (in absolute value) with increasing Rayleigh number, in accordance with the ${Ra}^{-1/4}$ scaling of the thermal boundary layers in figure 7. The same characteristic is also observed in the temperature distribution next to the interface. In the bulk of the bottom half (approximately $z=0.2\unicode{x2013}0.4$) and the top half (approximately $z=0.6\unicode{x2013}0.8$), a nearly uniform temperature profile is observed, with a very weak dependence on the Rayleigh number. The temperature values in these regions are approximately 0.42 at the bottom half and $-0.05$ at the top half, approximately constant in all the cases. This significant asymmetry is mainly attributed to the different thermal conductivities on each side of the interface. Since the thermal conductivity of the denser fluid is an order of magnitude larger than that of the lighter fluid, the bulk temperature at the bottom half is much closer to the temperature at the bottom heated wall than that of the top half of the cavity. A final observation from figure 9(a) is that the average temperature distributions show negligible sensitivity to the Weber number.

Figure 9.

Vertical distribution of (a) average and (b) r.m.s. values of the temperature field, for different values of the Rayleigh and Weber numbers. The line styles are the same as in figure 8.

The vertical distribution of the temperature r.m.s. fields is depicted in figure 9(b). The most distinct feature is the significantly smaller r.m.s. values at the bottom of the cavity compared to the top, hinting at a weaker turbulent state in the denser fluid, in line with the observations made from the instantaneous fields in § 3.1 and the Reynolds number predictions discussed in § 3.2. At the top of the cavity, two maxima are observed: a smaller one next to the top wall, and a larger one next to the two-phase interface, reinforced by the interface oscillation. For both peaks, the maximum values become smaller, and their location moves closer to the boundaries as the Rayleigh number increases. The same behaviour of the temperature r.m.s. field was observed in other single-phase studies, such as du Puits et al. (Reference du Puits, Resagk, Tilgner, Busse and Thess2007) for the turbulent convection in a cylindrical cell, and Demou & Grigoriadis (Reference Demou and Grigoriadis2019) in a cuboid cavity. Furthermore, increasing the Weber number leads to a noticeable increase in the temperature r.m.s. maximum values next to the interface. More specifically, the Weber number effects are stronger as the Rayleigh number increases. On the other hand, not surprisingly, the temperature r.m.s. values next to the top wall remain relatively unaffected with increasing Weber numbers.

Moving on to the statistics of the velocity field, we recall that since $\langle u \rangle _{x,y,t} = \langle v \rangle _{x,y,t} = \langle w \rangle _{x,y,t} = 0$, the velocity r.m.s. reduces to $u_{rms}=\sqrt {\langle u^2 \rangle _{x,y,t}}$, $v_{rms}=\sqrt {\langle v^2 \rangle _{x,y,t}}$ and $w_{rms}=\sqrt {\langle w^2 \rangle _{x,y,t}}$, which corresponds to the square root of the single-phase kinetic energy per unit mass. Therefore, the velocity r.m.s. values are associated with the average kinetic energy per unit mass of the large-scale vortical structures in each fluid layer. Figures 10(a,b) show the vertical distributions of the vertical and horizontal components of the velocity r.m.s. field, respectively. As already mentioned, the turbulence at the bottom fluid layer is weaker compared to the top fluid layer, something that is clearly depicted in the r.m.s. profiles. Still, the rotation of the large-scale vortical structures is adequate in retaining relatively large velocity r.m.s. values in the denser fluid, in comparison to the temperature r.m.s. values in figure 9(b).

Figure 10.

Vertical distribution of (a) the vertical and (b) the horizontal components of the velocity r.m.s. field, for different values of the Rayleigh and Weber numbers. The line styles are the same as in figure 8.

Focusing on the vertical velocity r.m.s. component, shown in figure 10(a), the different profiles exhibit a concave shape in each fluid layer, with a maximum at $z=0.25$ and $z=0.75$, i.e. in the middle of the two layers. The curvature at the maxima decreases with increasing Rayleigh number, indicating a change in the circulation of the large-scale vortical structures. The actual maximum values behave non-monotonically, with larger maxima for the $Ra=10^7$ case. Moreover, the effects of the Weber number are more visible in the lighter fluid. They are much more pronounced for the lower Rayleigh number cases, as can be expected considering the stronger interface deformation in these cases. More specifically, the maximum vertical velocity r.m.s. component in the lighter fluid for $Ra=10^6$ increases noticeably with increasing Weber number, while there is no significant shift in the higher Rayleigh number cases. On the other hand, the horizontal velocity r.m.s. components, shown in figure 10(b), exhibit a similar shape to the temperature r.m.s. profiles in figure 9(b), with two maxima in each layer; one next to the solid wall, and one next to the two-fluid interface. The difference compared to the temperature profiles is that the maximum velocity r.m.s. values are approximately symmetric in each layer. As concerns variations of interface deformability, we observe trends similar to those observed for the vertical velocity r.m.s. distribution, with a more pronounced effect of the Weber number for smaller Rayleigh numbers.

3.4. Surface displacement

The dynamics of a mechanically perturbed interface is thoroughly discussed in the literature, usually referred to as wave-turbulence (Falcon & Mordant Reference Falcon and Mordant2022). By applying a statistically stationary large-scale perturbation (through either wind or mechanical excitation), the interface deforms into waves, whose elevation $\eta (x,y,t)$ may vary significantly in time and space. Energy is redistributed towards small (dissipative) scales through nonlinear interactions, leading to a cascade process that is described thoroughly in the literature (Zakharov & Filonenko Reference Zakharov and Filonenko1967; Nazarenko Reference Nazarenko2011). Larger waves are generated through gravity forces, while smaller ones are governed by capillarity. The dispersion relation for the linear (small-amplitude) gravity-capillary waves reads (see Lamb Reference Lamb1993)

(3.6)\begin{equation} \hat{\omega}^2 = \left(\frac{\hat{\rho}_1-\hat{\rho}_2}{\hat{\rho}_1+\hat{\rho}_2}\,\hat{g} + \dfrac{\hat{\sigma}}{\hat{\rho}_1+\hat{\rho}_2}\,\hat{k}^2\right)\hat{k}, \end{equation}

where $\hat {\omega }=2{\rm \pi} \hat {f}$ is the angular frequency (corresponding to time period $\hat {T}_f=1/\hat {f}$), and $\hat {k}=2{\rm \pi} /\hat {\lambda }$ is the wavenumber, corresponding to wavelength $\hat {\lambda }$. The first term on the right-hand side of (3.6) represents the contribution due to the gravity forces, while the second is the modulation due to capillary stresses. These two contributions are equal at the capillary wavenumber $\hat {k}_c=\sqrt {\hat {g}(\hat {\rho }_1-\hat {\rho }_2)/\hat {\sigma }}$, which, using (3.6), gives the capillary frequency $\hat {f}_c=\hat {g}^{3/4}(\hat {\rho }_1-\hat {\rho }_2)^{3/4}(\hat {\rho }_1+\hat {\rho }_2)^{-1/2}\hat {\sigma }^{-1/4}/(\sqrt {2}{\rm \pi} )$. The crossover scale is the capillary wavelength $\hat {\lambda }_c=2{\rm \pi} \sqrt {\hat {\sigma }/((\hat {\rho }_1-\hat {\rho }_2)\hat {g})}$: waves longer than this are driven by gravity, while smaller waves are driven by capillarity. For both the capillary ($\hat {\lambda }>\hat {\lambda }_c$) and gravity ($\hat {\lambda }<\hat {\lambda }_c$) regimes, theoretical predictions are available for spectra computed in both time, $\hat {S}_\omega =\langle |\tilde {\eta }(x,y,\omega )|^2 \rangle _{x,y}$ and space $\hat {S}_k=\langle |\tilde {\eta }(k_x,k_y,t)|^2 \rangle _{t}$, with $(\tilde {\cdot })$ being the Fourier transform. For gravity waves, Zakharov & Filonenko (Reference Zakharov and Filonenko1966) found that $\hat {S}_\omega \sim \hat {\omega }^{-4}$ and $\hat {S}_k\sim \hat {k}^{-5/2}$, while for capillary waves, the scaling laws $\hat {S}_\omega \sim \hat {\omega }^{-17/6}$ and $\hat {S}_k\sim \hat {k}^{-15/4}$ were obtained by Zakharov & Filonenko (Reference Zakharov and Filonenko1967).

The dispersion relation (3.6) is obtained by solving the linearized system composed of the Euler equations (in Fourier space) for an incompressible, irrotational and inviscid velocity field, and the equation for the interface elevation dynamics. However, the present system is quite different, with gravity effects not only producing waves at the interface (as in the aforementioned references) but also affecting the velocity field through buoyancy due to the dependency of density on the temperature. In fact, the full system of equations for the present configuration should comprise the Navier–Stokes equations, including not only viscosity but, more importantly, the effect of density variations that lead to the production of vorticity by the baroclinic torque, making (3.6) indicative only for the present study. Furthermore, the linearization employed to obtain (3.6) assumes that the nonlinearities of the problem provide only weak contributions. If nonlinearities become stronger, then the dispersion relation widens (Aubourg & Mordant Reference Aubourg and Mordant2016) and eventually deviates from the theoretical prediction in (3.6) (Herbert, Mordant & Falcon Reference Herbert, Mordant and Falcon2010).

Given this short overview, we move on to examine the results of the present simulations. Figure 11 shows the temporal spectrum $\hat {S}_\omega$ of the wave elevation for each case, together with the scaling exponent predicted for gravity and capillary waves. When surface tension forces are higher, i.e. low $We$ (figures 11a–c), we observe that the capillary regimes ($\hat {\omega }>\hat {\omega }_c$) can be observed over an extended range of frequencies, especially for the highest Rayleigh numbers under investigation here (see figures 11b,c). For the larger Weber number (figures 11d–f), we observe a region with a flatter spectrum around $\hat {\omega }\sim \hat {\omega }_c$, with a secondary peak and a scaling behaviour close to the capillary regime, $\hat {\omega }^{-17/6}$. The secondary peak becomes less pronounced with increasing Rayleigh number, while no capillary scaling can be observed at large wavenumbers. While it is difficult to explain the origin of the secondary peak, this is likely not connected to gravity waves. In fact, by increasing the Weber number, we do not find any consistent trend, with either the capillary or gravity power laws. This is indeed counter-intuitive, because increasing the Weber number should lead to the formation of larger waves, as the lower surface tension forces cannot withstand the local velocity fluctuations. Such waves should, in principle, behave as gravity waves, as $\hat {\omega } < \hat {\omega }_c$, but this is not observed in the spectra.

Figure 11.

Temporal spectrum of the wave elevation $\eta (x,y,t)$. In each plot, we show the spectrum $\hat {S}_\omega$ (black lines), the $-4$ power law for gravity waves (red dashed lines), and the $-17/6$ power law for capillary waves (blue dash-dotted lines), for (a–c) all cases at ${We}=100$, and (d–f) all cases at ${We}=1000$. The $x$-axes are normalized by the capillary frequency $\hat {\omega }_c=2{\rm \pi} \hat {f}_c$.

A similar behaviour is also observed consistently in the spatial spectrum, shown in figure 12. Again, at low $We$ and high $Ra$ (figures 12b,c), the spatial spectra show a clean capillary wave range, which is observed only partially at low $Ra$ (figure 12a), and not observed at high $We$ ((figures 12d–f) except from the highest Rayleigh number. As shown in the temporal analysis, the spatial spectra also do not display a clean gravity waves range. In contrast, the secondary peak is clearly observed here for both Weber number values considered, although it is less visible at higher Rayleigh numbers. As evidenced, a pronounced secondary peak is associated with the disruption of the theoretical power-law behaviour of the gravity and capillary regimes.

Figure 12.

Spatial spectrum of the wave elevation $\hat {\eta }(x,y,t)$. In each plot, we show the spectrum $\hat {S}_k$ (black lines), the $-5/2$ power law for gravity waves (red dashed lines), and the $-15/4$ power law for capillary waves (blue dash-dotted lines), for (a–c) all cases at ${We}=100$, and (d–f) all cases at ${We}=1000$.

Next, we discuss the existence of the secondary peak. First, it is worth recalling that this is observed for frequencies (and wavenumbers) larger than the capillary scale. Similar deviations from the theoretical power laws are not unusual in wave turbulence (e.g. see Berhanu & Falcon Reference Berhanu and Falcon2013; Zonta, Soldati & Onorato Reference Zonta, Soldati and Onorato2015; Aubourg & Mordant Reference Aubourg and Mordant2016), and are typically attributed to dissipative effects or strong nonlinearities. In the present study, the most striking result is that the secondary peak in the spectrum attenuates at high Rayleigh number and low Weber number (a less pronounced peak appears in the spectrum), when a neater capillary wave regime can be observed. This observation can be linked to the interface fluctuations expressed through the r.m.s. values of the indicator function in figure 8(b), revealing smaller r.m.s. values at higher Rayleigh numbers. Considering as well the smaller values of the temperature r.m.s. distribution at higher Rayleigh numbers exhibited in figure 9(b) in combination with a less deformable interface at lower Weber numbers, the influence of nonlinear effects on the interface dynamics weakens.

To further investigate the interface dynamics, we display in figure 13 the space–time spectrum $E_\eta (k,\omega )$ for all flows under investigation, and compare those with the dispersion relation for capillary and gravity-capillary waves. In all cases, we observe that neither the capillary nor the gravity-capillary dispersion relations provide good predictions for the maxima of the spectra, which could be expected in the current system. Nevertheless, most of the measured maxima lie between the two theoretical dispersion relations, revealing an intermediate regime created most likely by the combination of nonlinear effects due to buoyancy and viscous effects at small scales. For all the studied cases, we note a widening of the dispersion relation, meaning that nonlinearities are a critical aspect of these flows. Such widening may also affect some of the maxima, which therefore lie outside the theoretical dispersion relation for the capillary and gravity-capillary laws. Furthermore, there is a clear range of low frequencies where the footprints of the forcing can be observed (see Aubourg & Mordant Reference Aubourg and Mordant2016). On the other hand, it is important to notice that as $Ra$ increases, nonlinearities seem to decrease, and the maxima converge towards a line. This suggests that the energy transport mechanism at high Rayleigh and low Weber numbers is characterized by reduced nonlinearities, which may help to predict the effective interface area across which heat/mass is transferred.

Figure 13.

Space–time spectrum $\hat {E}_\eta (k,\omega )$ of the wave elevation $\hat {\eta }(x,y,t)$. In each plot, we show the spectrum contour, and the theoretical dispersion relation for capillary waves $\hat {\omega }_c^2 = [\hat {\sigma }/(\hat {g}(\hat {\rho }_1+\hat {\rho }_2))]\hat {k}^3$ is shown with a black dashed line, while the gravity-capillary dispersion relation from (3.6) is shown with a red dashed curve. Vertical and horizontal dotted black lines mark the capillary angular frequency and wavenumber. Black crosses indicate the maximum of $\hat {E}_\eta$ for each angular frequency. The configuration of the plots follows figures 11 and 12, i.e. (a–c) all cases at ${We}=100$, and (d–f) all cases at ${We}=1000$. Increasing Rayleigh numbers are shown from left to right: (a,d) ${Ra}=10^6$, (b,e) ${Ra}=10^7$, and (c,f) ${Ra}=10^8$.