Two-Scalar Turbulent Rayleigh-Benard Convection: Numerical Simulations and Unifying Theory

We conduct direct numerical simulations for turbulent Rayleigh-B\'{e}nard (RB) convection, driven simultaneously by two scalar components (say, temperature and salt concentration) with different molecular diffusivities, and measure the respective fluxes and the Reynolds number. To account for the results, we generalize the Grossmann-Lohse theory for traditional RB convections~(Grossmann and Lohse, J. Fluid Mech., 407, 27-56; Phys. Rev. Lett., 86, 3316-3319; Stevens et al., J. Fluid Mech., 730, 295-308) to this two-scalar turbulent convection. Our numerical results suggest that the generalized theory can successfully predict the overall trends for the fluxes of two scalars and the Reynolds number. In fact, for most of the parameters explored here, the theory can even predict the absolute values of the fluxes and the Reynolds number with good accuracy. The current study extends the generality of the Grossmann-Lohse theory in the area of the buoyancy-driven convection flows.

above, and is subject to an external gravitational field. Such system has been extensively studied in recent years, e.g. see the reviews of Ahlers, Grossmann & Lohse (2009), Lohse & Xia (2010) and Chillà & Schumacher (2012). One key question is how the scalar flux and flow velocity depend on the control parameters. The unifying theory for the flux and flow velocity ('GL theory ' Grossmann & Lohse 2000, 2002, 2004, which are measured respectively by the Nusselt number Nu and the Reynolds number Re, has achieved great success for RB flows (Ahlers et al. 2009;Stevens et al. 2013), and now has predictive power for the absolute values of Nu and Re for given control parameters, i.e. the Rayleigh number Ra and the Prandtl number Pr.
However, in many cases the convection flow can be more complex than in the idealized RB system, e.g. as in the case of external rotation (e.g. King et al. 2009;Horn & Shishkina 2015;Wei, Weiss & Ahlers 2015), inhomogeneities of the top and bottom boundaries (e.g. Wang, Huang & Xia 2017;Bakhuis et al. 2018), or wall roughness (e.g. Roche et al. 2001;Shishkina & Wagner 2011;Salort et al. 2014;Wagner & Shishkina 2015;Xie & Xia 2017;Zhu et al. 2017;Jiang et al. 2018). In this study we will investigate another type of complexity, i.e. RB convection driven by two different scalar components. Multiple-component convection is commonly encountered in natural environments. For instance, the density of seawater is mainly determined by temperature and salinity, and chemical reaction flows usually have more than one species. In the ocean, the vertical convection flow driven by both temperature and salinity gradients is usually referred to as double diffusive convection (DDC) (Turner 1985;Radko 2013). Our previous study on DDC was confined in the so-called fingering regime, where the fluid layer experiences an unstable salinity gradient and a stable temperature gradient Yang, Verzicco & Lohse 2016). DDC can also happen in an electrodeposition cell with heat and ion concentration as two scalars (Hage & Tilgner 2010;Kellner & Tilgner 2014).
Here we will focus on convection flow driven simultaneously by two scalar components with different molecular diffusivities. Our previous study showed that the original GL model can be used to describe the salinity transfer in fingering DDC flow (Yang et al. , 2016. The theory has also been applied to DDC in the diffusive regime, in which the fluid layer is subjected to an unstable temperature gradient and a stable salinity gradient (Hieronymus & Carpenter 2016). In the present study the RB convection is driven by two scalar components which are both unstably stratified. Recall that the key idea of the GL theory is to divide both the momentum and thermal fields into their own boundary layer and bulk regions. Then scaling relations are developed for the two respective contributions to the respective dissipation rates, leading to Nu(Ra, Pr) and Re(Ra, Pr). It is straightforward to generalize this type of argument to multiple scalar fields. We will validate this generalization of the GL theory by direct numerical simulations. We stress that not a single new fitting parameter is introduced; we simply take those determined in Stevens et al. (2013).
The paper is organized as follows. In § 2 we will provide the governing dynamical equations of the system and the details of our simulations. In § 3 we will present the generalization and application of the GL theory to the two-scalar RB convection. Finally § 4 concludes the paper.

Governing equations and numerical simulations
The flow under consideration is incompressible and the density ρ is determined by two scalar components, say temperature θ (x, t) and concentration field s(x, t).
The Oberbeck-Boussinesq approximation is adopted, i.e. the fluid density depends linearly on both scalars as ρ(θ, s) = ρ 0 [1 − β θ θ + β s s]. ρ 0 is a reference density, while θ and s are the temperature and concentration relative to their respective reference values. β ζ is the positive expansion coefficient respectively for temperature (ζ = θ ) and concentration (ζ = s). From now on the subscript ζ = θ or s stands for a quantity associated with scalar ζ . The signs before the two terms indicate that the density is bigger for either lower temperature or higher concentration, which are the usual cases in practice.
The governing equations consist of the momentum equation and the advectiondiffusion equations for two scalars, which read Here u i with i = 1, 2, 3 are the three velocity components, p is the kinematic pressure, ν is the kinematic viscosity, g is the gravitational acceleration and κ ζ is the molecular diffusivity, respectively. The dynamic system is further constrained by the continuity equation ∂ i u i = 0. The fluid layer is between two parallel plates which are perpendicular to gravity and separated by a height H. At each plate both the temperature and concentration are kept constant. The scalar difference between two plates is denoted by ∆ ζ . The top plate has lower temperature and higher concentration, thus the flow is driven by both scalars. The flow has four control parameters, namely two Prandtl numbers and two Rayleigh numbers, Another useful parameter, which is borrowed from the DDC community, is the density ratio Λ = β θ ∆ θ β s ∆ s = Pr s Ra θ Pr θ Ra s . (2.3) Λ measures the relative strength of the buoyancy force induced by the temperature difference to that induced by the concentration difference. Λ < 1 indicates that the buoyancy force of the concentration difference is stronger than that of the temperature field, which we refer to as the 'concentration-dominant' (CD) regime. Accordingly, Λ > 1 is referred to as the 'temperature-dominant' (TD) regime. The three key responses of the system are the scalar fluxes and the flow velocity, which are measured by the two Nusselt numbers and the Reynolds number Here · h represents the average over time and a horizontal plane at arbitrary height.
In this work we calculate the two Nusselt numbers by the scalar gradients on the two plates. u rms is the root-mean-square value of the velocity magnitude calculated over the entire domain. The governing equations (2.1) are numerically solved by a highly efficient code developed in our group (Ostilla-Mónico et al. 2015), which has been used  (1) 192 (1)  For all cases Pr θ = 1. Columns from left to right: Prandtl number of concentration field, Rayleigh numbers of temperature and concentration, density ratio, aspect ratio of domain (horizontal length over height), horizontal resolution of base mesh (refinement factor of refined mesh), vertical resolution of base mesh (refinement factor of refined mesh), two Nusselt numbers and Reynolds number. Note that some cases appear repeatedly for the completeness of each group.
intensively in our previous DDC studies (Yang et al. , 2016. The code employs a multiple-grid method, which solves different flow quantities on either a base mesh or a refined mesh. Specifically, momentum and scalars with Pr 1 are always solved on the base mesh. Scalars with Pr > 1 usually require higher resolution and are solved on the refined mesh. The flow quantities are non-dimensionalized by the height H, the free-fall velocity defined by concentration difference U s = √ gβ s ∆ s H and the scalar differences ∆ ζ , respectively. At two plates no-slip boundary conditions are applied and both scalars are kept constant. In the two horizontal directions the periodic boundary conditions are employed. We fix Pr θ = 1 and change Ra θ , Ra s and Pr s . For Pr s < Pr θ (Pr s > Pr θ ) concentration diffuses faster (slower) than temperature.
The simulated cases are summarized in table 1. All cases are sorted into five groups, as shown in table 1. Within each group we only vary one control parameter and fix all Three-dimensional volume rendering of the temperature (a,c,e) and the concentration (b,d,f ) for three runs with different control parameters. The control parameters (Pr s , Pr θ , Ra s , Ra θ , Λ) are (10, 1, 10 7 , 10 5 , 0.1) for the (a,b), (10, 1, 10 7 , 10 7 , 10) for the (c,d) and (0.1, 1, 10 5 , 10 7 , 10) for the (e,f ), respectively. Note that given Pr s , Pr θ , Ra s and Ra θ , giving the value for Λ (by 2.3) is redundant information, but useful for the interpretation of the figures.
the others. To ensure the resolution is adequate, the mesh size is chosen to meet the criteria proposed in Stevens, Verzicco & Lohse (2010), i.e. for momentum the mesh size is smaller than the Kolmogorov length scale and for each scalar smaller than the Batchelor length scale, respectively. Furthermore, for the case with Pr s = 10, Ra s = Ra θ = 10 7 we run another simulation with a higher resolution n x (m x ) = 512(4) and n z (m z ) = 384(2) as a benchmark. Two different resolutions generate similar Nusselt and Reynolds numbers, with differences smaller than 1 %. In the table we also give the three responses of the flow, i.e. Nu s , Nu θ and Re. The uncertainty is measured by the standard deviation of temporal fluctuation.
In figure 1 we show the scalar fields for three runs with very different flow morphologies. The case shown in figure 1(a,b) has (Pr s , Ra s , Ra θ ) = (10, 10 7 , 10 5 ), which corresponds to Λ = 0.1. Thus, the buoyancy force is dominated by the concentration difference. However, for this case the molecular diffusivity of temperature is ten times bigger than that of concentration, and the typical size of the temperature plumes is much larger than the concentration ones due to the fast horizontal diffusion. In figure 1(c,d) we show the case with (Pr s , Ra s , Ra θ ) = (10, 10 7 , 10 7 ) and Λ = 10. Now the temperature difference dominates the buoyancy force, and the temperature plumes are very active. The concentration plumes become thin filaments embedded within temperature plumes due to the slow diffusion. Figure 1(e,f ) displays the case with (Pr s , Ra s , Ra θ ) = (0.1, 10 5 , 10 7 ), or Λ = 10. Although the temperature difference dominates the buoyancy force as is the case in figure 1(c,d), the concentration field has larger molecular diffusivity and therefore the concentration plumes have bigger horizontal size than the temperature plumes.   figure 2(b,d), Pr s = 10 and Ra s = 10 7 . Both scalar boundary layers are always nested inside the momentum one. As Ra θ increases from 10 5 in figure 2(b) to 10 7 in figure 2(d), all three boundary-layer thicknesses decrease due to the stronger buoyancy force. Interestingly, the concentration boundary layer also becomes thinner even though Ra s is the same for the two cases. The cases of panels figure 2(d,f ) have the same Ra θ = 10 7 and Λ = 10, and temperature difference dominates the buoyancy force. The temperature and momentum boundary layers have very similar heights for the two cases. However, among the three fields the thickness of the concentration boundary layer is the smallest for Pr s = 10 (figure 2d) and the largest for Pr s = 0.1 (figure 2f ).

Key idea of the original Grossmann-Lohse theory
The GL theory was originally developed for RB flow (Grossmann & Lohse 2000, 2002, 2004 to provide an unifying theory to account for the dependences Nu(Ra, Pr) and Re(Ra, Pr), and indeed successfully does so for most of the existing experimental and numerical results (Ahlers et al. 2009;Stevens et al. 2013). Here we will briefly describe the theory for completeness, and then discuss the extension to the current two-scalar convection flow. For full details the readers are referred to the original references.
The GL theory is built upon the exact relations between the dissipation rates and the global fluxes. For convection flow driven by an active scalar (e.g. temperature), these exact relations read (Shraiman & Siggia 1990) Here · represents the average over time and the entire domain. The key idea of the theory is to split the flow domain into boundary-layer and bulk regions and thus The individual contributions bl and bulk can be modelled for both momentum and thermal fields. For the bulk region, by using the Kolmogorov's energy-cascade picture, bulk can be estimated as the scaling laws of Re defined by a characteristic velocity of the bulk flow. For the boundary-layer region, bl can be estimated by assuming the Prandtl-Blasius-Pohlhausen profiles in the boundary layers and by using the scaling laws for the boundary-layer thicknesses. Then a cross-over function f is introduced to model the transition from the regime where the kinetic boundary layer is nested in the thermal one to the regime where the thermal boundary layer is thinner than the kinetic one. Another cross-over function g is employed to account for the saturation limit at small Re when the kinetic boundary-layer thickness is of the order of the domain height. By combining all the modelling, one obtains the original GL theory The model has five free coefficients a and c i with i = 1, 2, 3, 4. Re c can be calculated as 4a 2 . Stevens et al. (2013) explained in detail how these coefficients were fitted based on only four (experimental or numerical) data points Nu(Ra θ,j , Pr θ,j ), j = 1, 2, 3, 4 and one measure of Reynolds number, giving c 1 = 8.05, c 2 = 1.38, c 3 = 0.487, c 4 = 0.0252, a = 0.922.
(3.4a−e) Stevens et al. (2013) also showed that the theory successfully accounts for most of the existing RB results. Since for different flow configurations the Reynolds number can be defined by different characteristic velocities, to correctly predict the Reynolds number a transformation coefficient α needs to be determined and the procedure is described in Grossmann & Lohse (2002) and in Stevens et al. (2013). Such transformation does not change the theoretical prediction of the Nusselt number and merely accommodates the specific definition of Reynolds number in the flow investigated.
3.2. Extension to two-scalar convection and rational of the coefficients Extend the GL theory to the current two-scalar RB flow is straightforward. First of all, exact relations similar to (3.1) still exist and read They can be readily obtained from the dynamic equations of s 2 , θ 2 and u 2 /2 − gβ θ zθ + gβ s zs by the procedure given in Shraiman & Siggia (1990). We then notice that for the current flow the momentum and both scalar fields can still be divided into boundarylayer and bulk regions, as supported by figures 1 and 2. Thus for each field the dissipation rate can be modelled separated for the boundary-layer and bulk regions as in the original theory. Intuitively we assume that in the current flow the same scaling relations still hold for the dissipation rates in each region. By following the same arguments as in the original theory, one readily obtains the GL theory for the turbulent two-scalar convective flow, namely Equations (3.6a) without the first term and (3.6b) form the GL theory for the singlescalar RB flow, e.g. see (3.3). Comparing to the original theory, the first terms of (3.6a) and (3.6c) are introduced due to the existence of the second scalar. Note that (3.6b) and (3.6c) have the same coefficients. The physical reason for this necessity is that, as explained in Yang et al. (2015), when either of the two scalar differences decreases to zero, the flow must still reduce to the traditional RB flow with one scalar and the theory must degenerate to the original form. Thus the generalized theory still has the very same five coefficients c i with i = 1, 2, 3, 4 and a, and not two extra ones for the second scalar field, as one may naively (but erroneously) expect. Furthermore, the values of these coefficients can directly be taken from Stevens et al. (2013), i.e. those given in (3.4).

Application to the current numerical results
We now apply the theory to the current numerical results. As we explained before, a transformation coefficient α needs to be determined to correctly predict the Reynolds number for the current flow. Here we fix α = 1.453 by using the Reynolds number of the case with Pr s = 10, Ra θ = Ra s = 10 7 . Figure 3 shows both the theoretical predictions and the numerical measurements for the five groups of the cases listed in table 1. The overall trends of all three global responses, namely Nu θ , Nu s and Re, are very well captured by the theory. Moreover, the theory can even predict the absolute values for most of the cases with good accuracy. Specifically, the relative error of Nu θ prediction is always smaller than 20 %. For the Re prediction over 90 % of all cases have a relative error smaller than 20 %. The theory has the least accuracy for the Nu s prediction, but still the relative error for over 50 % of all cases are smaller than 20 %. The largest error is always smaller than 50 % for both Nu s and Re.
Among the three global responses, the theoretical prediction for the concentration flux Nu s exhibits the biggest deviation from the numerical results, especially in the deep TD regime (with large Λ) for Pr s > 1 (see the left panel of figure 3a-c,e). In this regime the flow is mainly driven by the temperature difference, such as the case shown in figure 1(c,d). All the concentration plumes are very thin and stay in the core regions of the temperature plumes, therefore the buoyancy anomaly associated with concentration and its effects on the momentum field are confined inside the temperature plumes due to the slow sideward diffusion. The morphology and dynamics of such thin concentration plumes are very different from the traditional RB plumes, thus the original scaling arguments of the GL theory become less accurate.
For the opposite reason in the deep CD regime (with small Λ), although the buoyancy force is mainly generated by the concentration component, the temperature anomaly diffuses faster in the lateral direction and is not confined inside the concentration plumes. The temperature anomaly can directly interact with the momentum field and thus more similar to the situation in the RB flows. Therefore the GL theory performs better in the CD regime than in the TD regime, as shown in figure 3(a-c).

Conclusions and discussion
In summary, we conducted direct numerical simulations of the RB convection driven by two scalar components. The flow morphology changes for different ratios of the buoyancy forces associated with the two scalar differences. We have generalized the GL theory for the single-scalar RB convection to the current problem. The results show that the theory captures the overall trends of the dependences of scalar fluxes and flow velocity on the control parameters. For most of the cases the theory predicts the absolute values with good accuracy. This comparison demonstrates the applicability of the GL theory to multiple-component convection flows. The accuracy of the theory decreases when temperature difference dominates the buoyancy force, i.e. in the TD regime, and when the concentration has large Prandtl number. We argue that in this regime, the structures of the concentration plumes are different from those in the traditional RB flows, and the argument in the original GL theory becomes less accurate. A more refined generalization of the GL theory should assume that the coefficients are not constant but some functions of the density ratio Λ, and they recover the original GL values when the system degenerates to a single-scalar RB flows. 10 6 10 7 10 8 10 9 10 5 10 3 10 2 10 6 10 7 10 8 10 9 10 5 10 6 10 7 10 8 10 5 10 4 10 6 10 7 10 8 10 5 10 4 10 6 10 7 10 8 10 0 10 -1 10 1 10 2 10 0 10 -1 10 1 10 2 10 0 10 -1 10 1 10 2 10 5 10 4 10 6 10 7 10 8 10 5 10 4 10 6 10 7 10 8 10 5 10 4 10 6 10 7 10 8 10 5 10 4 10 6 10 7 10 8 10 9 10 5 10 6 10 7 10 8 10 9 10 5 10 6 10 7 10 8 10 9