1. Introduction
Evaporation on a horizontal gas–liquid interface plays a pivotal role in different contexts, from geophysical processes such as moisture convection and vapour distribution within the atmosphere (Colman & Soden Reference Colman and Soden2021), to industrial applications such as the cooling of fuel rods in the spent-fuel pools of nuclear reactors (Hay & Papalexandris Reference Hay and Papalexandris2020). The presence of vapour in the gas changes the local thermophysical properties, in particular, the density and heat capacity (Colman & Soden Reference Colman and Soden2021), and thus modifies the global heat transfer in the system, quantified by the Nusselt number, 
$Nu$ (Schumacher & Pauluis Reference Schumacher and Pauluis2010). The mean vapour content in the gas phase depends on its value at the interface, which in turn is a function of the partial pressure and the interface temperature in an exponential fashion (e.g. Clausius–Clapeyron law, Span–Wagner relation). Thus, small changes in the interface temperature significantly impact the amount of vapour in the gas and the total heat transfer. A conceptually simple set-up to study these flows in a precise and controlled manner is the multiphase Rayleigh–Bénard (RB) configuration: two infinitely extended fluid layers confined by two horizontal walls at a fixed temperature, heated from below and cooled from above. Inspired from the classical single-phase counterpart used to model turbulent convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Chillà & Schumacher Reference Chillà and Schumacher2012), this configuration has been the object of numerical (Nataf, Moreno & Cardin Reference Nataf, Moreno and Cardin1988; Prakash & Koster Reference Prakash and Koster1994) and experimental (Xie & Xia Reference Xie and Xia2013; Zhang, Chong & Xia Reference Zhang, Chong and Xia2019) studies. In particular, the multiphase RB set-up has been recently adopted to study: (i) the interface breakup in the presence of buoyancy (Liu et al. Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a), (ii) the heat transfer enhancement due to the manipulation of the wall wettability (Liu et al. Reference Liu, Chong, Ng, Verzicco and Lohse2022a) and (iii) the modulation of heat transfer and interface temperature induced by the variation of the liquid layer height and thermal conductivity of the two phases (Liu et al. Reference Liu, Chong, Yang, Verzicco and Lohse2022b). All these numerical studies did not consider phase change and assumed constant thermophysical properties within the Oberbeck–Boussinesq (OB) approximation, with the exception of Biferale et al. (Reference Biferale, Perlekar, Sbragaglia and Toschi2012) for boiling flows and of Favier, Purseed & Duchemin (Reference Favier, Purseed and Duchemin2019) for ice melting.
Here, we include evaporation at the two-phase interface and relax the assumption of constant and uniform thermophysical properties in the gas phase while keeping those of the liquid uniform and constant. In particular, we extend the theory for the interface temperature proposed in Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022b) to account for (i) phase change at the interface and (ii) non-OB (NOB) effects in the gas phase induced by variations of the thermodynamic pressure, temperature and composition. We propose analytical scaling laws for predicting the interface temperature and the heat transfer modulation with respect to a RB system without evaporation. The resulting expressions are compared against high-fidelity direct numerical simulations (DNS), which are performed using a weakly compressible multiphase formulation with phase change, covering a substantial region of the 
$Ra$–
$\varepsilon$ parameter space.
This paper is organized as follows. In § 2, we introduce the main assumptions and derive analytical expressions for the interface temperature and heat transfer modulation in the evaporating RB system. In § 3, we describe the mathematical and numerical model employed to validate the analytical scaling laws. In § 4, we present the validation of the theory and an assessment of the assumptions behind the model. The main findings and conclusions are summarized in § 5.
2. Interface temperature and global heat transfer modulation
We consider a cavity partially filled with an evaporating single-component liquid and an initially dry gas, as shown in figure 1. The domain is laterally unbounded and confined by two horizontal walls separated by a distance 
$\hat {l}_z$ (
$\hat {\cdot }$ indicates a dimensional quantity). Constant temperatures 
$\hat {T}_b=(1+\varepsilon )\hat {T}_r$ and 
$\hat {T}_t=(1-\varepsilon )\hat {T}_r$ are imposed on the top heated and bottom cooled walls, where 
$\varepsilon =(\hat {T}_b-\hat {T}_t)/(2\hat {T}_r)= \widehat {{\rm \Delta} T}/(2\hat {T}_r)$ is the dimensionless temperature differential and 
$\hat {T}_r=(\hat {T}_t+\hat {T}_b)/2$ is the mean temperature, taken, hereinafter, as reference value. Under these conditions, the system eventually reaches a statistically stationary condition, with the vapour saturating the gas layer, leading to a dynamic balance between evaporation and condensation at the interface 
$\varGamma$. The presence of vapour dramatically changes the statistically stationary state: it modifies the gas thermophysical properties and hence the heat transfer inside the cavity, while reducing the height of the liquid layer. Here, we characterize these variations as a function of the amount of vapour inside the cavity. For this purpose, we write an expression for the mean interfacial vapour mass concentration 
$\bar {Y}_{l,\varGamma }^v$ (deduced from Raoult's law) and the Span–Wagner model for the interfacial vapour pressure 
$p_{s,\varGamma }$:
\begin{equation} \left.\begin{array}{c@{}} \bar{Y}_{l,\varGamma}^v =\dfrac{\lambda_Mp_{s,\varGamma}}{\lambda_Mp_{s, \varGamma}+(p_{th}-p_{s,\varGamma})},\\ p_{s,\varGamma} =\varPi_P^{{-}1}\exp[(B_1\eta_{sw}+B_2\eta_{sw}^{1.5} +B_3\eta_{sw}^{2.5}+B_4\eta_{sw}^{5})(1 -\eta_{sw})^{{-}1}].\end{array}\right\} \end{equation}
In (2.1), 
$\lambda _M=\hat {M}_l/\hat {M}_{g,r}$ is the molar mass ratio between liquid and inert gas and 
$\varPi _P=\hat {p}_{th,r}/\hat {p}_{cr}$ with 
$\hat {p}_{th,r}$ the reference thermodynamic pressure and 
$\hat {p}_{cr}$ the critical pressure. The quantity 
$\eta _{sw}=1-\hat {T}_{\varGamma }/\hat {T}_{cr}$ is the Span–Wagner parameter and the coefficients 
$B_{i=1,4}$ depend on the substance under consideration. By introducing the dimensionless interface temperature 
$\varTheta _{\varGamma } = (\hat {T}_{\varGamma }-\hat {T}_r)/\widehat {{\rm \Delta} T}$, 
$\eta _{sw}$ becomes
\begin{equation} \eta_{sw} = 1-\frac{\hat{T}_{\varGamma}}{\hat{T}_{cr}} =1-\varPi_T(1+2\varepsilon\varTheta_{\varGamma}),\end{equation}
where 
$\varPi _T=\hat {T}_r/\hat {T}_{cr}$. Equations (2.1) and (2.2) show that fixing the type of substance, 
$\bar {Y}_{l,\varGamma }^v$ is determined by five quantities: (i) the ratio between the mean temperature and the critical temperature, 
$\varPi _T$, (ii) the ratio between the reference thermodynamic pressure and the critical pressure, 
$\varPi _P$, (iii) the temperature differential 
$\varepsilon$, (iv) the interface temperature 
$\varTheta _\varGamma$ and (v) the thermodynamic pressure 
$p_{th}$. The first three quantities depend on the ambient conditions (typically given or measured) and the type of substance, whereas 
$\varTheta _{\varGamma }$ and 
$p_{th}$ depend on the flow in the two phases and are determined below.
Multiphase RB convection for 
$Ra=10^8$ and 
$\varepsilon =0.20$ (taken from the numerical simulations). Snapshots of the temperature distribution 
$\varTheta$ (a) at the start of evaporation (
$t=0$) and (b) after a statistically stationary state is achieved (
$t>300$). Snapshots of the vapour mass fraction, 
$Y_l^v$, distribution in the gas region at (c) 
$t\approx 60$, (d) 
$t\approx 122$ and (e) steady state, 
$t>300$ (
$t$ is scaled with the free-fall time 
$\hat {t}_{ff}=\hat {l}_z//(2\varepsilon | \hat {\boldsymbol {g}}|\hat {l}_z)^{0.5}$).
Our derivation relies on three central assumptions: (i) the OB approximation can be applied to the liquid phase, while the gas thermophysical properties are a generic function of thermodynamic pressure, local temperature and vapour composition, (ii) the Grossmann–Lohse (GL) theory for thermal convection can be applied to the liquid and gas layers separately and (iii) the vapour content in the gas can be taken as the mean value at the gas–liquid interface. The validity of these assumptions is assessed and discussed in § 4.
2.1. Interface temperature
The first step is to account for the variation of the liquid height induced by phase change and for the NOB effects (relevant for 
$\varepsilon \geq 0.05$; see Chillà & Schumacher Reference Chillà and Schumacher2012; Wan et al. Reference Wan, Wang, Wang, Xia, Zhou and Sun2020) due to variations of the local temperature, thermodynamic pressure and composition. Note that based on the first assumption of our derivation, NOB effects manifest only in the gas phase, whereas the liquid phase is described under the OB approximation. For simplicity of notation, a generic quantity 
$\hat {\xi }_g$ is expressed as
\begin{equation} \hat{\xi}_g =\hat{\xi}_{g,r}\left(1+\frac{{\rm \Delta} \hat{\xi}_g}{\hat{\xi}_{g,r}}\right)= \hat{\xi}_{g,r}f_{g,\xi},\end{equation}
where 
${\rm \Delta} \hat {\xi }_g=\hat {\xi }_g-\hat {\xi }_{g,r}$ and 
$f_{g,\xi }$ is the mean normalized variation of 
$\hat {\xi }_g$ with respect to the same quantity evaluated at reference condition, 
$\hat {\xi }_{g,r}$. Both 
$\hat {\xi }_{g,r}$ and 
$f_{g,\xi }$ refer to the layer pertaining to the gas phase. Following the approach proposed by Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022b) and using (2.3), we define a Rayleigh number in each phase:
\begin{equation} \left.\begin{array}{c@{}} Ra_l = \dfrac{\hat{\beta}_l(\hat{T}_b-\hat{T}_{\varGamma}) |\hat{\boldsymbol{g}}|\hat{h}_l^3\hat{\rho}_l^2 \hat{c}_{pl}}{\hat{\mu}_l\hat{k}_l}=\alpha_0^3(1/2- \varTheta_{\varGamma})Ra\dfrac{\lambda_{\rho}^2 \lambda_{cp}\lambda_{\beta}}{\lambda_{\mu}\lambda_k}f_{l,h}^3, \\ Ra_g = \dfrac{\hat{\beta}_g(\hat{T}_{\varGamma}-\hat{T}_t) |\hat{\boldsymbol{g}}|\hat{h}_g^3\hat{\rho}_g^2\hat{c}_{pg}}{\hat{\mu}_g \hat{k}_g}=(1-\alpha_0)^3(1/2+\varTheta_{\varGamma})Ra \dfrac{f_{g,h}^3f_{g,\rho}^2f_{g,cp}}{f_{g,\mu}f_{g,k}}, \end{array}\right\} \end{equation}
where 
$\hat {h}_i$ is the height of each layer, 
$\hat {\beta }_{i}$ is the thermal expansion coefficient, 
$\hat {\rho }_{i}$, 
$\hat {\mu }_{i}$, 
$\hat {k}_{i}$ and 
$\hat {c}_{p,i}$ are the fluid density, dynamic viscosity, conductivity and specific heat capacity and 
$\lambda _{\xi }$ is the associated property ratio scaled with respect to the reference gas property evaluated at 
$\hat {T}_r$, 
$\hat {p}_{th,r}$ and in a dry condition, i.e. 
$Y_l^v=0$. Further, 
$|\hat {\boldsymbol {g}}|$ is the gravitational acceleration, 
$\alpha _0$ is the initial liquid volume fraction and 
$Ra=\hat {\beta }_g\widehat {{\rm \Delta} T}|\hat {\boldsymbol {g}}|\hat {l}_z^3\hat {\rho }_{g,r}^2\hat {c}_{pg,r}/(\hat {\mu }_{g,r}\hat {k}_{g,r})$ is a ‘fictitious’ Rayleigh number based on the reference gas thermophysical properties, the height of the cavity and the temperature difference between top and bottom walls. Given the NOB effects in the gas, we define 
$\hat {\beta }_g=1/\hat {T}_{c,g}$. In contrast, 
$\hat {\beta }_l$ is taken as constant and independent of 
$\varepsilon$ in the liquid, where we impose the OB approximation. Temperature 
$\hat {T}_{c,g}$ is the central temperature in the gas region and its estimation is given later in this section.
Next, we define two separate Nusselt numbers, 
$Nu_l$ and 
$Nu_g$:
\begin{equation} \left.\begin{array}{c@{}} Nu_l =\dfrac{\hat{Q}_{\varGamma,l}\hat{h}_l}{\hat{k}_l(\hat{T}_b -\hat{T}_\varGamma)}=\dfrac{\hat{Q}_{\varGamma,l}\alpha_0 \hat{l}_zf_{l,h}}{\hat{k}_{l}(1/2-\varTheta_{\varGamma})\widehat{{\rm \Delta} T}}, \\ Nu_g =\dfrac{\hat{Q}_{\varGamma,g}\hat{h}_g}{\hat{k}_g(\hat{T}_\varGamma- \hat{T}_t)}=\dfrac{\hat{Q}_{\varGamma,g}(1-\alpha_0) \hat{l}_zf_{g,h}}{\hat{k}_{g,r}f_{g,k}(1/2+\varTheta_{\varGamma}) \widehat{{\rm \Delta} T}}, \end{array}\right\} \end{equation}
where 
$\hat {Q}_{\varGamma,l}$ and 
$\hat {Q}_{\varGamma,g}$ are the heat fluxes on the liquid and gas side of the interface. In the absence of phase change and when evaporation and condensation events are statistically balanced, these two quantities are on average equal. Employing the GL theory, the Nusselt number in both layers is then related to the corresponding Rayleigh number. Note that the complete GL theory is a system of equations that provides the value of the Nusselt 
$Nu$ and the Reynolds 
$Re$ numbers for given Rayleigh 
$Ra$ and Prandtl 
$Pr$ numbers. The implicit nature of this system prevents obtaining an explicit expression for 
$\varTheta _\varGamma$ and, therefore, the following simplified scaling laws are here considered (Weiss et al. Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018):
\begin{equation} Nu_l = A_lRa_l^{\gamma_l} Pr_l^{m_l}, \quad Nu_g = A_gRa_g^{\gamma_g} Pr_g^{m_g}. \end{equation}
In this work, we consider 
$A=A_l=A_g$, 
$\gamma =\gamma _l=\gamma _g$ and 
$m=m_l=m_g$. As remarked in Weiss et al. (Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018), employing the simplified GL theory in (2.6a,b) is valid as long as 
$Ra_l$, 
$Ra_g$ and 
$Pr_l$, 
$Pr_g$ are sufficiently similar to fall inside the same scaling regime (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001) so that the same 
$\gamma$ and 
$m$ can be used for both layers. Taking the ratio 
$Nu_l/Nu_g$ yields
\begin{equation} \frac{Nu_l}{Nu_g} =\left(\frac{Ra_l}{Ra_g}\right)^\gamma \left(\frac{Pr_l}{Pr_g}\right)^m.\end{equation}
As suggested in Weiss et al. (Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018), if 
$Pr>0.5$ the GL theory suggests a scaling exponent 
$m$ very close to zero and the Prandtl dependence in (2.7) can be omitted. To obtain an explicit relation for 
$\varTheta _\varGamma$, we compute the ratio 
$Ra_l/Ra_g$ from (2.4):
\begin{equation} \frac{Ra_l}{Ra_g} =\frac{\alpha_0^3}{(1-\alpha_0)^3} \frac{f_{l,h}^3}{f_{g,h}^3}\frac{(1/2-\varTheta_\varGamma)}{(1/2+ \varTheta_\varGamma)}\frac{\mathcal{F}_\lambda}{f_{g,\rho}^2 \mathcal{F}_g}\frac{\lambda_\beta}{\lambda_k}f_{g,k},\end{equation}
where 
$\mathcal {F}_\lambda =\lambda _{\rho }^2\lambda _{cp}/\lambda _\mu$ and 
$\mathcal {F}_g=f_{g,cp}/f_{g,\mu }$. Likewise, we express 
$Nu_l/Nu_g$ using (2.5) as
\begin{equation} \frac{Nu_l}{Nu_g} =\frac{\alpha_0}{(1-\alpha_0)} \frac{f_{l,h}}{f_{g,h}}\frac{1}{\lambda_k} \frac{(1/2+\varTheta_\varGamma)}{(1/2 -\varTheta_\varGamma)}f_{g,k}.\end{equation}
Note that to derive (2.9), the heat fluxes on the liquid and gas side are taken equal 
$\hat {Q}_{\varGamma,l}=\hat {Q}_{\varGamma,g}$. Once more, this is a valid assumption when evaporation and condensation balance at the interface and the gas layer is at saturation. By employing the scaling relations (2.6a,b) and (2.8) and (2.9), we get
\begin{equation} \frac{\alpha_0}{(1-\alpha_0)}\frac{f_{l,h}}{f_{g,h}} \frac{1}{\lambda_k}f_{g,k}\frac{(1/2+\varTheta_\varGamma)}{(1/2- \varTheta_\varGamma)}=\left[\frac{\alpha_0^3}{(1-\alpha_0)^3} \frac{f_{l,h}^3}{f_{g,h}^3}\frac{(1/2-\varTheta_\varGamma)}{(1/2+ \varTheta_\varGamma)}\frac{\mathcal{F}_\lambda}{f_{g,\rho}^2 \mathcal{F}_g}\frac{\lambda_\beta}{\lambda_k}f_{g,k}\right]^\gamma, \end{equation}which, after some manipulation, reads
\begin{equation} \frac{1}{\varTheta_\varGamma^*} = 1+\left(\frac{\alpha_0}{1-\alpha_0} \frac{f_{l,h}}{f_{g,h}}\right)^{({1-3\gamma})/({1+\gamma})} \left(\frac{f_{g,\rho}^2\mathcal{F}_g}{\mathcal{F}_\lambda \lambda_\beta}\right)^{{\gamma}/({1+\gamma})} \left(\frac{f_{g,k}}{\lambda_k}\right)^{({1-\gamma})/({1+\gamma})}.\end{equation}
Note that in the derivation of (2.11), we have performed a change of variable, 
$\varTheta _\varGamma ^*=\varTheta _\varGamma +1/2$. Equation (2.11) can be then rearranged in terms of 
$\varTheta _\varGamma ^*$ and, finally, of 
$\varTheta _\varGamma$:
\begin{equation} \varTheta_{\varGamma} ={-}\frac{1}{2}+\left(1+\left(\frac{\alpha_0}{1- \alpha_0}\frac{f_{l,h}}{f_{g,h}}\right)^{({1-3\gamma})/({1+\gamma})} \left(\frac{f_{g,\rho}^2\mathcal{F}_g}{\mathcal{F}_{\lambda} \lambda_{\beta}}\right)^{{\gamma}/({1+\gamma})} \left(\frac{f_{g,k}}{\lambda_k}\right)^{({1-\gamma})/({1+\gamma})} \right)^{{-}1}\!,\end{equation}
where 
$\lambda _\beta$ is defined as 
$\lambda _\beta =\hat {\beta }_l/\hat {\beta }_g$ with 
$\hat {\beta }_g=1/\hat {T}_{c,g}$. Accordingly, 
$\lambda _\beta$ can be finally expressed as
\begin{equation} \lambda_\beta = \hat{\beta}_l\hat{T}_{c,g}=\hat{\beta}_l\widehat{{\rm \Delta} T}\left(\frac{\hat{T}_{c,g}-\hat{T}_r}{\widehat{{\rm \Delta} T}}+\frac{\hat{T}_r}{\widehat{{\rm \Delta} T}}\right)=\hat{\beta}_l\hat{T}_r(1+2\varepsilon\varTheta_c). \end{equation}
In (2.12), 
$f_{g,k}$ and 
$\mathcal {F}_{g}$ are functions of temperature and composition. The dependence on the thermodynamic pressure is typically important for the density, i.e. 
$f_{g,\rho }$, while it can be omitted for the specific heat capacity, viscosity and thermal conductivity. Therefore, to evaluate 
$f_{g,k}$ and 
$\mathcal {F}_{g}$ in (2.12), we need to specify a reference temperature and reference vapour concentration. This aspect has already been discussed in Weiss et al. (Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018), where the authors derive an expression for the central temperature for the gas region 
$\hat {T}_{c,g}$ and show that it represents the temperature at which the thermophysical properties should be evaluated for a RB cell under strong NOB effects. Note that the derivation is based on the same assumptions that lead to (2.6a,b) and (2.7) and, therefore, no additional hypotheses are introduced here. By incorporating this approach in our model, 
$\hat {T}_{c,g}$ reads as
\begin{equation} \hat{T}_{c,g} =\frac{\hat{k}_{+,g}^{3/4}\hat{\eta}_{+,g}^{1/4} \hat{T}_\varGamma+\hat{k}_{-,g}^{3/4}\hat{\eta}_{-,g}^{1/4} \hat{T}_t}{\hat{k}_{+,g}^{3/4}\hat{\eta}_{+,g}^{1/4}+ \hat{k}_{-,g}^{3/4}\hat{\eta}_{-,g}^{1/4}}.\end{equation}
Note that when the gas thermophysical properties are uniform, 
$\hat {T}_{c,g}$ reduces to 
$(\hat {T}_\varGamma +\hat {T}_t)/2$ and, therefore, (2.14) can be interpreted as a more general choice of the central temperature than the arithmetic mean. By introducing 
$\varTheta _{c,g}=(\hat {T}_{c,g}-\hat {T}_r)/\widehat {{\rm \Delta} T}$, (2.14) can be written in dimensionless form as
\begin{equation} \varTheta_{c,g} =\frac{\varTheta_\varGamma-\dfrac{1}{2} \left(\dfrac{\hat{k}_{-,g}^3\hat{\eta}_{-,g}}{k_{+,g}^3 \hat{\eta}_{+,g}}\right)^{1/4}}{1+\left(\dfrac{\hat{k}_{-,g}^3 \hat{\eta}_{-,g}}{\hat{k}_{+,g}^3\hat{\eta}_{+,g}}\right)^{1/4}}.\end{equation}
Note that in (2.14) and (2.15), the group 
$\hat {\eta }_{\pm,g}$ corresponds to 
$(\hat {\beta }\hat {\rho }^2/\hat {k}\hat {\mu })_{\pm,g}$. The properties 
$\hat {k}_{+,g}$, 
$\hat {\eta }_{+,g}$ and 
$\hat {k}_{-,g}$, 
$\hat {\eta }_{-,g}$ are evaluated at the crossover points, which are located at the transition points between the boundary layer and the bulk region of the cell. Following once more the procedure in Weiss et al. (Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018), the crossover temperature 
$\hat {T}_{\pm,g}$ at which these properties should be evaluated is determined as a linear combination between 
$\hat {T}_{\varGamma }$ and 
$\hat {T}_t$. In particular, for the gas region we have
\begin{equation} \hat{T}_{+,g} = \delta \hat{T}_\varGamma + (1-\delta)\hat{T}_{c,g}, \quad \hat{T}_{-,g} = \delta \hat{T}_t +(1-\delta)\hat{T}_{c,g}. \end{equation}
The last parameter to be specified is 
$\delta$, which depends exponentially on the aspect ratio 
$\mathcal {A}$, i.e. 
$\delta = \delta _0 + (1-\delta _0)\exp {(-B\mathcal {A})}$. By experimental fitting, Weiss et al. (Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018) suggest employing 
$\delta _0=0.235$ and 
$B=1.14$, which provide an accurate estimation of 
$\varTheta _{c,i}$ for a wide range of 
$Ra$ and 
$\mathcal {A}$.
For the terms 
$f_{g,\rho }$, 
$f_{l,h}$ and 
$f_{g,h}$, more analysis is needed. Since evaporation changes the mass of the gas and its local density, and it decreases the volume of the liquid, we introduce the mass ratio 
$G_g=\hat {G}_g/\hat {G}_{g,r}$ and the volume ratio 
$V_g=\hat {V}_g/\hat {V}_{g,r}$. These ratios are defined as the values of the quantities in the evaporating regime divided by the corresponding values for the flow without evaporation. Indicating with 
$\hat {G}_{l,r}$ the initial liquid mass, 
$\hat {V}_{l,r}$ the initial liquid volume and 
$\hat {G}_l^e$ the mass of the evaporated liquid, we get
\begin{equation} \left.\begin{array}{c@{}} \displaystyle \dfrac{\hat{G}_g}{\hat{G}_{g,r}} = \dfrac{\hat{G}_{g,r}+\hat{G}_{l,r}-\hat{G}_l}{\hat{G}_{g,r}} = 1 + \dfrac{\hat{G}^e_{l}}{\hat{G}_{g,r}} = 1 + \dfrac{1}{\hat{G}_{g,r}}\displaystyle{\int_{\hat{V}_g} \hat{\rho}_gY_l^v\widehat{{\rm d} V}_g},\\ \displaystyle \dfrac{\hat{V}_g}{\hat{V}_{g,r}} = \dfrac{\hat{V}_{g,r}+\hat{V}_{l,r}-\hat{V}_l}{\hat{V}_{g,r}} = 1 + \dfrac{\hat{V}_l^e}{\hat{V}_{g,r}} = 1 + \dfrac{1}{\lambda_{\rho}\hat{G}_{g,r}} \displaystyle{\int_{\hat{V}_g}\hat{\rho}_gY_l^v\widehat{{\rm d} V}_g}, \end{array}\right\}\end{equation}
where 
$\hat {V}_l^e$ is the volume of the incompressible liquid that turns into vapour. Note that 
$\hat {G}_l^e$, given by the integrals in (2.17), requires knowledge of the vapour distribution. Approximating 
$Y_l^v$ with 
$\bar {Y}_{l,\varGamma }^v$ (second assumption of our derivation) allows us to estimate the mass of the evaporated liquid as 
$\hat {G}_l^e\approx \bar {Y}_{l,\varGamma }^v\hat {G}_g$. Accordingly, 
$G_g=1/(1-\bar {Y}_{l,\varGamma }^v)$ and 
$V_g=1+\bar {Y}_{l,\varGamma }^v/(\lambda _{\rho }(1-\bar {Y}_{l,\varGamma }^v))$, which provide the following estimation of 
$f_{g,\rho }$:
\begin{equation} f_{g,\rho} = \frac{G_g}{V_g} = \frac{\lambda_{\rho}}{\bar{Y}_{l,\varGamma}^v +(1-\bar{Y}_{l,\varGamma}^v)\lambda_{\rho}}.\end{equation}
Since 
$f_{i,h}\approx \hat {G}_i/(\hat {\rho }_i\hat {V}_{i,r})$, the relations (2.17) allow also the estimation of 
$f_{l,h}$ and 
$f_{g,h}$.
To proceed, it is worth noticing that 
$\bar {Y}_{l,\varGamma }^v=g_1(\varPi _T,\varPi _P,\varepsilon, \varTheta _{\varGamma },p_{th})$ from (2.1) and 
$\varTheta _{\varGamma }=g_2(\varPi _T,\varPi _P,\varepsilon,p_{th})$ as in (4.2) and, therefore, a relation for the thermodynamic pressure (supposed uniform; see Chillà & Schumacher Reference Chillà and Schumacher2012) is required. To derive it, we integrate over the gas region the equation of state for the local gas density, i.e. 
$p_{th}M_{m}/(1+2\varepsilon \varTheta _g)$, and express the result in terms of 
$p_{th}$:
\begin{equation} p_{th} ={G_g}\left({{\int_{V_g}M_{m}\varTheta_g^idV_g}}\right)^{{-}1} \approx \frac{f_{g,\rho}}{\bar{M}_{m}}\varTheta_c^i, \end{equation}
where 
$\varTheta _g^i=(1+2\varepsilon \varTheta _g)^{-1}$ and 
$\bar {M}_{m}=\lambda _M/(\bar {Y}_{l,\varGamma }^v +(1-\bar {Y}_{l,\varGamma }^v)\lambda _M)$ is the mean molar mass of the mixture computed using the harmonic average between 
$\hat {M}_l$ and 
$\hat {M}_{g,r}$ (Scapin et al. Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022). Note that in (2.19), we employ the second hypothesis and we approximate the volume integral of 
$\varTheta _g^i$ as 
$V_g\bar {\varTheta }_g^i$ using 
$\bar {\varTheta }_g^i=\varTheta _c^i$ from (2.15).
The model to estimate 
$\varTheta _\varGamma$ is based on (2.12), (2.15), (2.18) and (2.19), coupled with appropriate equations of state for specific heat capacity, thermal conductivity and viscosity, as detailed in Appendix B. The system is not linear; however, a simple iterative procedure can be used to obtain 
$\varTheta _\varGamma$, 
$\bar {Y}_l^v$ and 
$p_{th}$ together with 
$\varTheta _{c,g}$. We remark here that the proposed model for 
$\varTheta _\varGamma$ is more general than that presented in Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022b) in three aspects: (i) we account for phase change, (ii) we account for NOB effects in the gas phase and (iii) we include the density, viscosity, specific heat capacity and thermal expansion ratios in (2.12). It is worth mentioning that in the absence of evaporation, with uniform bulk properties (i.e. 
$\mathcal {F}_g=f_{g,\xi }=1$) and for 
$\lambda _\rho =\lambda _\mu =\lambda _{cp}=\lambda _{\beta }=1$, the general expression for 
$\varTheta _\varGamma$ in (2.12) reduces to the estimate by Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022b).
2.2. Nusselt number
With the estimated interface temperature 
$\varTheta _{\varGamma }$, we can derive a scaling law for the ratio between the global Nusselt number with and without evaporation, i.e. 
$Nu^e=\hat {Q}_t^e\hat {l}_z/(\hat {k}_{g,r}\widehat {{\rm \Delta} T})$ and 
$Nu=\hat {Q}_t\hat {l}_z/(\hat {k}_{g,r}\widehat {{\rm \Delta} T})$, where 
$\hat {Q}_t^e$ and 
$\hat {Q}_t$ are the global heat fluxes with and without evaporation measured at the top boundary. Taking the ratio between the two, we immediately see that 
$Nu^e/Nu = \hat {Q}_t^e/\hat {Q}_t$. To compute 
$\hat {Q}_t^e/\hat {Q}_t$, we first define the Nusselt number on the gas side of the interface, with and without evaporation:
\begin{equation} Nu_g^e = \frac{\hat{Q}_{\varGamma}^e\hat{h}_g^e}{\hat{k}_g^e (\varTheta_{\varGamma}^e+1/2)\widehat{{\rm \Delta} T}}, \quad Nu_{g} = \frac{\hat{Q}_{\varGamma}\hat{h}_g}{\hat{k}_{g} (\varTheta_{\varGamma}+1/2)\widehat{{\rm \Delta} T}},\end{equation}
where 
$\hat {Q}_{\varGamma }^e$ and 
$\hat {Q}_{\varGamma }$ are the heat fluxes exchanged at the interface. Taking the ratio of the two Nusselt numbers in (2.20a,b) and applying again the GL theory (i.e. 
$Nu_g^e/Nu_g=(Ra_g^e/Ra_g)^{\gamma }$) yields
\begin{equation} \frac{\hat{Q}_{\varGamma}^e}{\hat{Q}_{\varGamma}} = \frac{(\varTheta_{\varGamma}^e+1/2)}{(\varTheta_{\varGamma}+1/2)} \left(\frac{Ra_{g}^e}{Ra_g}\right)^{\gamma} \frac{f_{g,k}^e}{f_{g,k}}\frac{1}{f_{g,h}^e}.\end{equation}
Note that the variations of the thermophysical properties induced by evaporation are denoted with a superscript ‘
$e$’. We then use (2.4) to compute the ratio
\begin{equation} \frac{Ra_g^e}{Ra_{g}} =\frac{(1+2\varepsilon\varTheta_c^e)}{(1 +2\varepsilon\varTheta_c)}\frac{(\varTheta_{\varGamma}^e +1/2)}{(\varTheta_{\varGamma}+1/2)}f_{g,h}^{3,e}f_{g,\rho}^{2,e} \frac{f_{g,cp}^{e}}{f_{g,\mu}^{e}f_{g,k}^{e}} \frac{f_{g,\mu}f_{g,k}}{f_{g,cp}}.\end{equation}
Since the global heat flux is equal to the heat flux at the interface, i.e. 
$\hat {Q}_{\varGamma }^e=\hat {Q}_{t}^e$ and 
$\hat {Q}_{\varGamma }=\hat {Q}_{t}$, we can finally combine equations (2.21) and (2.22) into
\begin{equation} \frac{Nu^e}{Nu} = \left(\frac{1+2\varepsilon\varTheta_c}{1 +2\varepsilon\varTheta_c^e}\right)^{\gamma}\left(\frac{\varTheta_{\varGamma}^e +1/2}{\varTheta_{\varGamma}+1/2}\right)^{1+\gamma} \frac{f_{g,\rho}^{2\gamma,e}}{f_{g,h}^{1-3\gamma,e}} \frac{f_{g,cp}^{\gamma,e}}{f_{g,cp}^{\gamma}} \frac{f_{g,\mu}^{\gamma}}{f_{g,\mu}^{e,\gamma}} \frac{f_{g,k}^{1-\gamma,e}}{f_{g,k}^{1-\gamma}}.\end{equation}
Solving iteratively the system composed of (2.12) and (2.19), together with the relations (2.18) and (2.1), provides the values of 
$\varTheta _\varGamma$ and 
$p_{th}$. Once these are known, we can directly estimate the global heat transfer modulation with (2.23). This expression predicts 
$Nu$ for an evaporating system with respect to the configuration without phase change, described by the GL theory. Inspection of (2.23) allows us to draw some initial conclusions concerning the role of phase change in the RB system. First, the change in liquid height influences the heat transfer modulation only weakly, since its exponent is 
$1-3\gamma \approx 0$. Second, higher gas density decreases the interface temperature, as suggested by (4.2). Nevertheless, despite that in (2.23) the exponent of 
$\varTheta _\varGamma$ is larger than the exponent of 
$f_{g,\rho }$, this last term is expected to be dominant given its stronger dependence on 
$\varepsilon$, as is clearly shown in the next section. Last, a non-negligible effect is present due to the variation of 
$c_p$, 
$\mu$ and 
$k$ whose contribution to the heat transfer modulation scales with 
$\gamma$ and 
$1-\gamma$.
3. Numerical methodology
3.1. Governing equations
The validation of the model previously described is performed with the in-house code for phase-changing flows extensively described in Scapin, Costa & Brandt (Reference Scapin, Costa and Brandt2020) and Scapin et al. (Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022) and, therefore, we briefly mention here only the main features. First, to distinguish between the phases, an indicator function 
$H$ is introduced, defined equal to 
$1$ in the liquid phase and 
$0$ in the gas phase. Function 
$H$ is governed by the following transport equation:
\begin{equation} \frac{\partial H}{\partial t} + \boldsymbol{u}_\varGamma\boldsymbol{\cdot}\boldsymbol{\nabla} H = 0,\end{equation}
where 
$\boldsymbol {u}_\varGamma$ is the interface velocity computed as the sum of an extended liquid velocity, 
$\boldsymbol {u}_l^e$, and a term due to phase change, 
$\dot {m}_\varGamma /\rho _l\boldsymbol {n}_\varGamma$, with 
$\dot {m}_\varGamma$ the mass flux and 
$\boldsymbol {n}_\varGamma$ the unit normal vector at the interface. Note that 
$\boldsymbol {u}_l^e$ is computed as described in Scapin et al. (Reference Scapin, Costa and Brandt2020). Next, the numerical code solves the governing equations assuming that the liquid phase has constant properties and can be treated within the OB approximation, while the gas phase manifests compressible effects that can be described within the low-Mach-number formulation. Accordingly, the dimensionless conservation equations for momentum, vaporized species 
$Y_l^v$, temperature 
$\varTheta$ and mass flux 
$\dot {m}_{\varGamma }$ across the interface read (Scapin et al. Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022)
\begin{gather} \rho\frac{{\rm D}\boldsymbol{u}}{{\rm D}t} ={-}\boldsymbol{\nabla} p+\sqrt{\frac{Pr}{Ra}}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}+ \frac{\boldsymbol{f}_{\sigma}}{We}+\left[\lambda_{\rho} (1-\varPi_{\beta}2\varepsilon\varTheta)H+ \frac{p_{th}M_{m}}{(1+2\varepsilon\varTheta)}(1-H)\right] \frac{\boldsymbol{e}_z}{2\varepsilon}, \end{gather}
\begin{gather}\rho_g\frac{DY_{l}^{v}}{{\rm D}t} = \frac{\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v})}{\sqrt{RaSc}} , \end{gather}
\begin{gather}\rho c_p\frac{D\varTheta}{{\rm D}t} = \frac{\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla}\varTheta)}{\sqrt{RaPr}} + \varPi_{R}\frac{{\rm d} p_{th}}{{\rm d} t}(1-H) - \frac{(\dot{m}_{\varGamma}\delta_{\varGamma})}{2\varepsilon Ste}, \end{gather}
\begin{gather}\dot{m}_{\varGamma}=\frac{1}{\sqrt{RaSc}}\frac{\rho_{g, \varGamma}D_{lg,\varGamma}}{1-Y_{l,\varGamma}^v} \nabla_{\varGamma} Y_{l}^{v}\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}. \end{gather}
In (3.2), 
$\boldsymbol {u}$ is the velocity, 
$p$ is the hydrodynamic pressure, 
$\tau$ is the viscous stress tensor for compressible Newtonian flows and 
$\boldsymbol {f}_{\sigma }=\kappa _{\varGamma }\delta _{\varGamma } \boldsymbol {n}_\varGamma$ with 
$\kappa _{\varGamma }$ the interfacial curvature, 
$\delta _\varGamma$ a regularized Dirac-delta function (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999) and 
$\boldsymbol {n}_\varGamma$ the normal vector. The unit vector 
$\boldsymbol {e}_z$ points in the gravity direction, i.e. 
$\boldsymbol {e}_z=(0,0,-1)$.
The generic thermophysical property 
$\xi$ (density 
$\rho$, dynamic viscosity 
$\mu$, thermal conductivity 
$k$ or specific heat capacity 
$c_p$) is computed with an arithmetic average, i.e. 
$\xi =1+(\lambda _{\xi }-1)H$, where 
$\lambda _{\xi }=\xi _l/\xi _{g,r}$ with 
$\xi _{g,r}$ evaluated at reference condition (i.e. 
$\hat {T}_r$, 
$\hat {p}_{th,r}$ and 
$Y_l^v=0$). Since 
$\xi _l$ is kept constant and uniform, no further modelling is needed, while the generic gas property 
$\xi _{g}$ is computed with the appropriate equation of state. For example, the gas density is computed with the ideal gas law, 
$\rho _g=p_{th}M_{m}/(1+2\varepsilon \varTheta )$, and the vapour diffusion coefficient with the Wilke–Lee correlation (Reid, Prausnitz & Poling Reference Reid, Prausnitz and Poling1987), 
$D_{lg}=(1+2\varepsilon \varTheta )^{3/2}/p_{th}$ (Wan et al. Reference Wan, Wang, Wang, Xia, Zhou and Sun2020). The remaining gas thermophysical properties are computed as detailed in Appendix B. Note that in the OB limit, (i.e. 
$\varepsilon \rightarrow 0$, 
$p_{th}=1$ and 
$M_{m}=1$), the gravity term active in the gas region reduces to 
$1-2\varepsilon \varTheta$ with 
$\hat {\beta }_g=1/\hat {T}_r$ and hence matches that employed in previous works (Liu et al. Reference Liu, Ng, Chong, Lohse and Verzicco2021b,Reference Liu, Chong, Wang, Ng, Verzicco and Lohsea, Reference Liu, Chong, Ng, Verzicco and Lohse2022a,Reference Liu, Chong, Yang, Verzicco and Lohseb).
Equations (3.2)–(3.5) are written in dimensionless form by introducing the free-fall velocity scale 
$\hat {u}_r=(2\varepsilon |\hat {\boldsymbol {g}}|\hat {l}_z)^{1/2}$ and the free-fall time scale 
$\hat {t}_r=\hat {l}_z/\hat {u}_r$. Accordingly, we define the Weber number 
$We=\hat {\rho }_{g,r}2\varepsilon |\hat {\boldsymbol {g}}| \hat {l}_{z}^2/\hat {\sigma }$ with 
$\hat {\sigma }$ the surface tension, 
$\varPi _{\beta }=\widehat {\beta _l}\widehat {T_r}$; 
$Sc=\hat {\mu }_{g,r}/(\hat {\rho }_{g,r}\hat {D}_{lg,r})$ and 
$Pr=\hat {\mu }_{g,r}\hat {c}_{pg,r}/\hat {k}_{g,r}$ are the Schmidt and the Prandtl numbers. Based on the chosen reference quantities, the Rayleigh number 
$Ra=2\varepsilon |\hat {\boldsymbol {g}}|\hat {l}_z^3 \hat {\rho }_{g,r}^2\hat {c}_{pg,r}/(\hat {\mu }_{g,r}\hat {k}_{g,r})$ corresponds to the fictitious one defined in (2.4). Note that the temperature equation (3.4) requires the definition of the Stefan number 
$Ste=\hat {c}_{pg,r}\hat {T}_r/\hat {{\rm \Delta} h}_{lv}$, where 
$\hat {{\rm \Delta} h}_{lv}$ is the latent heat, and of the dimensionless group 
$\varPi _{R}=\hat {R}_u/(\hat {c}_{pg,r}\hat {M}_{g,r})$, where 
$\hat {M}_{g,r}$ is the molar mass of the gas phase and 
$\hat {R}_u$ the universal gas constant. In (3.5), the vapour mass fraction at the interface 
$Y_{l,\varGamma }^v$ is computed using (2.1). These are Raoult's law (Reid et al. Reference Reid, Prausnitz and Poling1987) and the Span–Wagner equation of state for the vapour pressure 
$\hat {p}_{s,\varGamma }$ at the gas–liquid interface. The employed coefficients are those for pentane and taken equal to 
$B_{i=1,4}=[-7.327140,+1.823650,-2.272744,-2.711929]$. Note that we prefer to employ the slightly more elaborated Span–Wagner model over ‘simpler’ equations for the partial pressure, e.g. the Clausius–Clapeyron law or Antoine's law (Reid et al. Reference Reid, Prausnitz and Poling1987). The motivation behind our choice is twofold. First, the Span–Wagner equation is based on critical quantities, 
$\hat {T}_{cr}$ and 
$\hat {p}_{cr}$, which are intrinsic properties of the substance and thus independent of the local ambient conditions (e.g. 
$\hat {p}_{th,r}$). Next, the Span–Wagner model provides accurate and reliable results for most of the substances over a wide range of temperature and thermodynamic pressure, well below the critical point as well as near it.
To form a closed set of equations, one needs a relation for the velocity divergence and for the thermodynamic pressure:
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} = f_{\varGamma} + \left[\frac{p_{th}f_{Y} + f_\varTheta}{p_{th}} - \left(1-\frac{\varPi_{R}}{c_pM_{m}}\right)\frac{1}{p_{th}} \frac{{\rm d} p_{th}}{{\rm d} t}\right](1-H), \end{gather}
\begin{gather}\frac{1}{p_{th}}\frac{{\rm d} p_{th}}{{\rm d} t}\int_{V_g}\left(1- \frac{\varPi_{R}}{c_pM_{m}}\right) {\rm d} V_g= \int_V\left[f_{\varGamma}+\frac{f_\varTheta+p_{th} f_Y}{p_{th}}(1-H)\right]\, {\rm d}V . \end{gather}
The local divergence constraint, (3.6), is derived by applying the divergence operator to the one-fluid velocity defined as 
$\boldsymbol {u}=H\boldsymbol {u}_l+(1-H)\boldsymbol {u}_g$ and by applying the continuity equation of both phases. Equation (3.7) is derived by integrating equation (3.6) over the total domain 
$V$, sum of the liquid and gas domains, and by imposing the volume conservation over 
$V$, i.e. 
$\int _V\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {u}\, \textrm {d}V=0$. In (3.6) and (3.7) the functions 
$f_{\varGamma }$, 
$f_{Y}$ and 
$f_{\varTheta }$ represent the different contributions to the total velocity divergence from the phase change (
$f_{\varGamma }$) and the change of the gas density due to composition (
$f_Y$) and temperature (
$f_{\varTheta }$) (see again Scapin et al. (Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022) for details):
\begin{gather} f_{\varGamma} = \dot{m}_{\varGamma}\left(\frac{1}{\rho_{g,\varGamma}}- \frac{1}{\lambda_{\rho}}\right)\delta_{\varGamma}, \end{gather}
\begin{gather}f_{Y}=\frac{1}{\sqrt{RaSc}}\frac{M_{m}}{\rho_g} \left(\frac{1}{\lambda_M}-1\right)\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v}), \end{gather}
\begin{gather}f_{\varTheta} =\frac{\varPi_{R}}{c_pM_{m}}\frac{1}{\sqrt{RaPr}}\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla}\varTheta). \end{gather}
We remark that a weakly compressible formulation is still required to model an evaporating RB cell for 
$\varepsilon \rightarrow 0$. Indeed, the term 
$f_\varGamma$ is not negligible in the case of evaporation and contributes to expansion and contraction at the two-phase interface.
The governing equations (3.2)–(3.5) are solved on a uniform Cartesian grid with a standard marker and cell method (Harlow & Welch Reference Harlow and Welch1965). The interface dynamics is captured with an algebraic volume-of-fluid method, VoF-MTHINC (Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012; Rosti, De Vita & Brandt Reference Rosti, De Vita and Brandt2019), which ensures excellent conservation properties of the liquid volume both with and without phase change, provided (3.6) is correctly imposed on 
$\boldsymbol {u}$. Thus, a pressure correction method is employed, and the associated Poisson equation is first factorized into a constant-coefficient one (Dodd & Ferrante Reference Dodd and Ferrante2014) and then solved with the eigenexpansion technique (Schumann & Sweet Reference Schumann and Sweet1988), as implemented in the open-source code CaNS (Costa Reference Costa2018). Further details and validations are provided for phase-change problems with constant and variable properties in Scapin et al. (Reference Scapin, Costa and Brandt2020), Dalla Barba et al. (Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021) and Scapin et al. (Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022). Appendix A reports an additional validation against the multiphase RB convection case in Liu et al. (Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a). Note that the baseline two-phase code, including the VoF-MTHINC and heat transfer effects, FluTAS, is described in Crialesi-Esposito et al. (Reference Crialesi-Esposito, Scapin, Demou, Rosti, Costa, Spiga and Brandt2023) and is released as open-source software.
3.2. Computational set-up
For the validation of the model, we consider three values of the Rayleigh number, 
$10^6$, 
$10^7$ and 
$10^8$, and four values for the temperature differential 
$\varepsilon =0.05$, 
$0.10$, 
$0.15$ and 
$0.20$. At fixed dimensionless mean temperature 
$\varPi _T$ and pressure 
$\varPi _P$, the temperature differential is the only parameter affecting 
$\varTheta _\varGamma$, 
$Y_l^v$ and 
$p_{th}$. The remaining dimensionless parameters are not varied; we choose the Prandtl, Schmidt and Stefan numbers equal to unity, i.e. 
$Pr=Sc=1$ and 
$Ste=1$, and the property ratios as 
$\lambda _{\rho }=\lambda _{\mu }=\lambda _k=20$, 
$\lambda _{cp}=1$ and 
$\lambda _M=2.58$. Moreover, we set 
$\varPi _T=0.8$, 
$\varPi _{P}=1.65$, 
$\varPi _{\beta }=0.6$ and 
$\varPi _{R}=0.18$. This choice corresponds to a light hydrocarbon (e.g. pentane) at high temperature (below the critical value) and high pressure, which can be used as a coolant in industrial applications. Finally, we set the Weber number 
$We=5$ to limit interface deformation for all the investigated values of the Rayleigh number (Liu et al. Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a). All the cases are first simulated without evaporation until a statistically stationary condition. Following the procedure proposed for NOB flows (Demou & Grigoriadis Reference Demou and Grigoriadis2019), temporal convergence is assessed by comparing the Nusselt number at the bottom and top wall, ensuring that the relative difference is lower than 
$1\,\%$. Next, evaporation is activated until a new statistically stationary regime is reached. Also in this case, temporal convergence is assessed by comparing the top and bottom values of 
$Nu$. First- and second-order statistics of the generic variable 
$g$ (denoted as 
$\langle g\rangle _x$ and 
$\langle g_{rms}\rangle _x$ with 
$x$ the periodic direction) are collected for a sampling period sufficient to ensure their being independent of the size of the sample. Note that we use both Favre and Reynolds averaging. Unless otherwise stated, only the latter is employed in the current work, as we found the difference between the two negligible. In table 1, we report the sample size and the time step employed for each case with and without evaporation.
Time window for statistical sampling (
$T_{avg}$) and the fixed time step 
${\rm \Delta} t_{avg}$ employed to collect the statistics, for both the cases with (WT) and without (WO) evaporation (EV). Time is reported in units of free-fall time 
$\hat {t}_{ff}$. The cases where the gas density and the gas–liquid diffusion coefficient are the only variable properties are conducted at 
$Ra=10^6$, 
$10^7$ and 
$10^8$, and for 
$\varepsilon =0.05$, 
$0.10$, 
$0.15$ and 
$0.20$. The cases where all the gas thermophysical properties are varied are conducted at 
$Ra=10^6$ and 
$10^8$, and for 
$\varepsilon =0.05$, 
$0.10$, 
$0.15$ and 
$0.20$.
The computational domain is a two-dimensional cavity with aspect ratio 
$\mathcal {A}=2$, periodic in the horizontal direction and with two walls at the bottom and top. Here, a no-slip condition is applied for the velocity, a Dirichlet condition for the temperature, while zero flux is imposed on the remaining quantities. In all the cases, we employ the same uniform grid with 
$N_x=1024$ and 
$N_z=512$ along the periodic and the wall-normal directions. This resolution fulfils the two requirements proposed in Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). Taking as a reference the case at 
$Ra=10^8$ and 
$\varepsilon =0.20$ (the most demanding among our cases in terms of grid resolution), we perform two checks summarized below.
First, we ensure that the grid size 
$\hat {{\rm \Delta} }=\hat {l}_z/N_z$ is smaller than the Kolmogorov length scale:
\begin{equation} \hat{{\rm \Delta}}_b \leq {\rm \pi}\hat{\eta} \approx {\rm \pi}\hat{l}_z\sqrt[4]{\frac{Pr^2}{RaNu}}. \end{equation}
Using the a posteriori estimation of the Nusselt number 
$Nu\approx 40$ at 
$Ra=10^8$ and 
$\varepsilon =0.20$, we get 
$\hat {{\rm \Delta} }/({\rm \pi} \hat {\eta })\approx 0.20$; thus the requirement is well met. When 
$\mu$, 
$\rho$ and 
$c_p$ are not uniform, 
$Pr_g$ slightly increases above unity up to 
$1.53$. Therefore, even though the Batchelor scale is smaller than the Kolmogorov scale, the limited 
$Pr$ makes the employed resolution suitable for our configuration.
Next, we estimate the minimum number of grid points required to fully resolve the thermal and hydrodynamic boundary layers of the gas phase 
$N_{gp}$, which represents the most stringent condition. Following once more Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010), this requirement reads as
\begin{equation} N_{gp} = \sqrt{2}aNu_g^{1/2}Pr^{0.3215+0.011\log(Pr)},\end{equation}
with 
$a=0.482$ in (3.10) and the Nusselt number on the gas side 
$Nu_g$ estimated as
\begin{equation} Nu_g =\frac{\hat{Q}_g\hat{h}_g}{\hat{k}_{g,r}(\hat{T}_\varGamma- \hat{T}_t)}=Nu\frac{\alpha_0f_{g,h}}{(\varTheta_\varGamma+1/2)}. \end{equation}
Taking 
$\alpha _0=0.5$, 
$f_{g,h}=1.05$ and 
$\varTheta _\varGamma =0.36$ in (3.10), 
$N_{gp}$ is equal to 
$5$, significantly less than the value of 
$10$ employed in the current set-up. Finally, given the spatial variation of 
$\rho _g$ and 
$D_{lg}$ with the state variables, the local 
$Sc$ may differ from unity. Taking once more the case at 
$Ra=10^8$ and 
$\varepsilon =0.2$, the mean gas density and the vapour diffusion coefficient become 
$1.8$ and 
$0.9$ times the corresponding reference values, corresponding to an effective 
$Sc_{eff}\approx 0.6$. Since 
$Sc_{eff}< Pr$, the velocity and thermal fields impose the stricter resolution requirement.
We have performed a mesh convergence study for the case at 
$Ra=10^8$ and 
$\varepsilon =0.2$ by doubling the grid in both directions (i.e. 
$2048\times 1024$) and comparing the result with the corresponding coarser simulation. As shown in figure 2, the chosen resolution (
$1024\times 512$) guarantees excellent spatial convergence for first- and second-order temperature statistics, confirming once more that it can be considered as adequate for the current study.
Grid convergence studies for the case at 
$Ra=10^8$ and 
$\varepsilon =0.20$: (a) mean vertical profile and (b) root mean square of temperature using 
$1024\times 512$ and 
$2048\times 1024$ grid points.
4. Results and discussion
4.1. Temporal evolution of 
$Nu$ and $\varTheta _\varGamma$
All the cases are first simulated without evaporation until a statistically stationary equilibrium is reached. Once this condition is met, evaporation is activated at the gas–liquid interface. During the phase-change process, the liquid height is reduced, and the gas region changes its mean temperature, composition and pressure due to increased vapour content. Moreover, phase change introduces another heat transfer component, i.e. latent heat, which is responsible for the mismatch between the Nusselt number values measured at the top and the bottom of the cell, 
$Nu_t$ and 
$Nu_b$, during the transient phase. This condition is clearly displayed in figure 3 for 
$Ra=10^6$ and 
$10^8$, and 
$\varepsilon =0.05$ and 
$0.20$. In particular, 
$Nu_b$ rapidly increases since evaporation is an endothermic process, while 
$Nu_t$ rapidly reduces since less heat is transported from the interface to the upper wall. Note that, irrespective of 
$Ra$, the transient between the two statistical equilibria is faster for larger values of 
$\varepsilon$ since more vapour is released in the gas layer. Eventually, the gas layer saturates with a balance of evaporation and condensation at the gas–liquid interface, and a new statistical equilibrium is achieved. In this new condition, the time-averaged mass flux is zero, i.e.
\begin{equation} \int_{t_{eq}}^{t_{eq}+T}\left(\int_\varGamma \dot{m}_\varGamma \, {\rm d} S\right)\, {\rm d} t=0, \end{equation}
where 
$t_{eq}$ is the physical time at which a statistical equilibrium is reached and 
$T$ is the time window employed for statistical sampling. Accordingly, the latent heat exchanged at the interface is statically zero, only sensible heat is exchanged and, thus, the Nusselt number values measured at the top and bottom walls converge to the same statistical mean. This final condition is also displayed in figure 3 for 
$Ra=10^6$ and 
$10^8$, and 
$\varepsilon =0.05$ and 
$0.20$. Higher 
$\varepsilon$ leads to a larger 
$\bar {Y}_l^v$ and, thus, to an enhancement of the heat transfer in the cell. Moreover, large values of 
$\bar {Y}_l^v$ increase the fluctuations of the Nusselt number around its mean value, an effect more pronounced at lower 
$Ra$.
Temporal evolution of the Nusselt number 
$Nu$ for (a) 
$Ra=10^6$ and (b) 
$10^8$, considering 
$\varepsilon =0.05$ and 
$0.20$, measured at the bottom wall (BW) and top wall (TW). The time instant 
$t_{eq}=0$ refers to the instant when evaporation at the gas–liquid interface is activated.
Figure 4 displays the time evolution of the interface temperature 
$\varTheta _\varGamma$ for the same cases. When evaporation is active, the interface cools and 
$\varTheta _\varGamma$ suddenly drops due to the latent heat. The final statistically stationary condition is reached when condition (4.1) is met. Similarly to the temporal behaviour of 
$Nu$, larger values of 
$\varepsilon$ cause more significant fluctuations of the interface temperature 
$\varTheta _\varGamma$ and reduce the transient before the saturation condition is met. It is worth mentioning that a different choice of Stefan number would affect the system's behaviour during the transient phase. On the other hand, different values of 
$Ste$ would lead to the same statistical equilibrium, i.e. 
$Nu^e/Nu$ and 
$\varTheta _\varGamma ^e$ when the condition (4.1) is met.
Temporal evolution of the interface temperature 
$\varTheta _\varGamma$ for (a) 
$Ra=10^6$ and (b) 
$10^8$, considering 
$\varepsilon =0.05$ and 
$0.20$. The instant 
$t_{eq}=0$ refers to the time that evaporation at the gas–liquid interface is activated.
4.2. Validation of the model
The model described in § 2 is here validated against two-dimensional interface-resolved DNS of the evaporating system. We first study the variation of 
$\rho$, 
$c_p$, 
$\mu$ and 
$k$, together with the molar mass 
$\bar {M}_m$, with the temperature differential 
$\varepsilon$. The results, displayed in figure 5, are computed by combining the analytical model described in § 2.1 together with the equations of state reported in Appendix B. For this reason, a non-negligible dependence on the scaling exponent 
$\gamma$ is observed, especially for the gas density and the molar mass. The current set-up considers a mixture of air and light hydrocarbon with molar mass ratio 
$\lambda _M$ larger than 1. Therefore 
$\bar {M}_m$ increases with 
$\varepsilon$. Also 
$f_{g,\rho }$, computed with (2.18), increases with the temperature differential 
$\varepsilon$, but, as shown in figure 5, the variation of 
$f_{g,\rho }$ is dominant over 
$\bar {M}_m$. Figure 5 displays the dependence of dynamic viscosity, thermal conductivity and specific heat capacity with 
$\varepsilon$. Since the specific heat capacity of the gas, 
$c_p$, increases with the mean vapour content, its normalized variation, 
$f_{g,cp}>1$, and increases with 
$\varepsilon$. The opposite occurs for 
$\mu _g$ and 
$k_g$; therefore, their normalized variations 
$f_{g,\mu }$ and 
$f_{g,k}$ are lower than unity with a decreasing trend with 
$\varepsilon$.
(a) Variation of the normalized mean density 
$f_{g,\rho }$ and molar mass 
$\bar {M}_m$ as a function of 
$\varepsilon$. (b) Variation of the normalized mean dynamic viscosity 
$f_{g,\mu }$, thermal conductivity 
$f_{g,k}$ and specific heat capacity 
$f_{g,cp}$. Parameters 
$f_{g,i}$ with 
$i=\rho$, 
$\mu$, 
$k$ and 
$c_p$ are computed using the definition, i.e. (2.3), and the equations of state detailed in Appendix B.
To better disentangle the effects of evaporation in the RB cell, we start the discussion by considering a subset of the general model, where the gas density and the liquid–gas diffusion coefficient are the only thermophysical properties that vary. This simplification allows us to isolate the effect of the density variations and omit, without further approximations, the Prandtl number dependence in the scaling 
$Nu_g\sim Ra_g^\gamma$, since the Prandtl number is not a function of 
$\rho _g$. In this simplified setting, 
$f_{g,\mu }=f_{g,cp}=f_{g,k}=1$ and (2.12) for the interface temperature 
$\varTheta _{\varGamma }$ becomes
\begin{equation} \varTheta_{\varGamma}
={-}\frac{1}{2}+\left(1+\left(\frac{\alpha_0}{1-
\alpha_0}\frac{f_{l,h}}{f_{g,h}}\right)^{({1-3\gamma})/({1+\gamma})}
\left(\frac{f_{g,\rho}^2\lambda_{\mu}}{\lambda_{\rho}^2
\lambda_{cp}\lambda_{\beta}}\right)^{{\gamma}/({1+\gamma})}
\left(\frac{1}{\lambda_k}\right)^{({1-\gamma})/({1+
\gamma})}\right)^{{-}1}.\end{equation}
Similarly, the Nusselt number ratio 
$Nu^e/Nu$ reads
\begin{equation} \frac{Nu^e}{Nu} = \left(\frac{1+2\varepsilon\varTheta_c}{1+ 2\varepsilon\varTheta_c^e}\right)^{\gamma}\left(\frac{\varTheta_{\varGamma}^e +1/2}{\varTheta_{\varGamma}+1/2}\right)^{1+\gamma} \frac{f_{g,\rho}^{2\gamma,e}}{f_{g,h}^{1-3\gamma,e}}.\end{equation}
Figure 6 compares the analytical prediction of 
$\varTheta _{\varGamma }$ and 
$Nu^e/Nu$ against the values extracted from the DNS. For both quantities, the predicted results agree very well with the theory, for a choice of scaling exponent 
$1/4\leq \gamma \leq 1/3$, as suggested by the GL theory (grey region).
Comparison between the analytical predictions (with 
$\gamma \in [1/4,1/3]$) and the numerical simulations for (a) the interface temperature 
$\varTheta _{\varGamma }$ and (b) the ratio 
$Nu^e/Nu$ as a function 
$\varepsilon$ for different values of 
$Ra$, when only the gas density is varied. Note that in (a), the dotted lines correspond to the prediction of 
$\varTheta _\varGamma$ without evaporation.
The data also reveal different aspects of the role of the temperature differential 
$\varepsilon$ in 
$\varTheta _\varGamma$. First, 
$\varTheta _{\varGamma }$ shows no or little dependence on 
$\varepsilon$ in a dry environment (i.e. no evaporation). This suggests that strong density variations have a negligible impact on 
$\varTheta _\varGamma$ and 
$Nu$. Second, when evaporation is active, the gas mixture has higher density and 
$\varTheta _{\varGamma }$ decreases with respect to the case without evaporation, exhibiting a stronger-than-linear dependence on 
$\varepsilon$. Based on the proposed model, the increase of the gas density is expected to occur regardless of the value of 
$\lambda _M$ (equation (2.18)). Nonetheless, 
$\lambda _M$ affects the statistically steady value of 
$p_{th}$. For 
$\lambda _M<1$, typical of mixtures of water vapour and inert gas, 
$\bar {M}_{m}<1$ and, therefore, variations of density and molar mass, 
$f_{g,\rho }$ and 
$\bar {M}_{m}$, contribute to the increase of 
$p_{th}$ (equation (2.19)). The opposite occurs for 
$\lambda _M>1$ as in mixtures of hydrocarbons and inert gases, where there is a competitive effect between 
$f_{g,\rho }$ and 
$\bar {M}_{m}$, with the former dominating the latter in the current set-up (see figure 5). Finally, note that the increase in the gas density with 
$\varepsilon$, i.e. 
$f_{g,\rho }>1$, compensates the drop in the interface temperature and, thus, contributes to an enhancement of the global heat transfer, i.e. 
$Nu^e/Nu>1$. We now consider the general case when all the gas thermophysical properties are varied with temperature and composition. Note that the dependence on the thermodynamic pressure is assumed only for the gas density. Based on the parameters and the equations of state we employ for 
$\mu$, 
$c_p$ and 
$k$ (reported in Appendix B), the temperature dependence is weaker than the effect of the composition. The variation of 
$\bar {Y}_l^v$ drives the change of 
$c_p$, 
$\mu$ and 
$k$ compared with the dry case. More specifically and as anticipated in figure 5, 
$c_{pg}$ increases with 
$\varepsilon$, i.e. 
$f_{g,cp}>1$, while 
$k$ and 
$\mu$ decrease with 
$\varepsilon$, i.e. 
$f_{g,\mu }<1$ and 
$f_{g,k}<1$. Note that from (2.12), 
$f_{g,cp}>1$ and 
$f_{g,\mu }<1$ promote a reduction of the interface temperature in the case of evaporation. Conversely, 
$f_{g,k}<1$ promotes an increase in 
$\varTheta _\varGamma$. For the current choice of parameters, the changes in viscosity and heat capacity have a stronger effect compared with the thermal conductivity. This factor leads to a decrease of 
$\varTheta _\varGamma$ with 
$\varepsilon$. Figure 6 displays the analytical predictions for these two quantities and the results from DNS conducted for 
$Ra=10^6$ and 
$10^8$, and 
$\varepsilon =0.05$, 
$0.10$, 
$0.15$ and 
$0.20$. As for the other cases, the scaling exponent 
$\gamma$ is chosen in the interval 
$[1/4, 1/3]$. Compared with the simplified setting where the gas density is the only variable thermophysical property, 
$\varTheta _\varGamma$ exhibits slightly lower values (see figure 6a). This behaviour is attributed to the larger sensitivity of the interface temperature to 
$c_p$ and 
$\mu$ than to variations of 
$k$. By comparing figure 6(b) with figure 7(b), we note that accounting for the variability of 
$\mu$, 
$c_p$ and 
$k$ leads to a more significant increase of the heat transfer, 
$Nu^e/Nu$. As shown in (2.23), the increase of 
$Nu^e/Nu$ with 
$\varepsilon$ is driven by the decrease in the heat capacity and viscosity ratio. This increase is only partially compensated by the increase in the thermal conductivity ratio.
Comparison between the analytical predictions (with 
$\gamma \in [1/4,1/3]$) and the numerical simulations for (a) the interface temperature 
$\varTheta _{\varGamma }$ and (b) the ratio 
$Nu^e/Nu$ as a function 
$\varepsilon$ for different values of 
$Ra$, when all the gas thermophysical properties are varied. Note that in (a), the dotted lines correspond to the prediction of 
$\varTheta _\varGamma$ without evaporation.
It is worth noticing that in all the cases and irrespective of which gas thermophysical properties vary, the solution does not converge to the case without evaporation for 
$\varepsilon \rightarrow 0$, i.e. 
$\varTheta _\varGamma ^e/\varTheta _\varGamma \nrightarrow 1$ and 
$Nu^e/Nu \nrightarrow 1$ for 
$\varepsilon \rightarrow 0$. This can be explained by considering the Span–Wagner equation of state as in (2.1) and, in particular, the parameter 
$\eta _{sw}=1-\varPi _T(1+2\varepsilon )$. When 
$\varepsilon$ is reduced, 
$\eta _{sw}$ approaches 
$1-\varPi _T$ and, therefore, some vapour still reaches the gas region, i.e. 
$\bar {Y}_{l,\varGamma }^v>0$. The presence of vapour changes the local composition, modifies the bulk thermophysical properties and affects both 
$\varTheta _\varGamma$ and 
$Nu^e$ even for 
$\varepsilon \rightarrow 0$.
4.3. Assessment of the hypotheses
We finally assess the validity of the assumptions invoked at the beginning. With regards to the first assumption, it is reasonable to neglect the variation of the liquid thermophysical properties. This assumption can be easily relaxed by using the mathematical framework proposed here to include appropriate equations of the state for the liquid thermophysical properties.
Moving to the second assumption, previous studies in single-phase and multiphase thermal convection have already proven that the GL theory accurately predicts the Nusselt number in the case of NOB effects (Weiss et al. Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018) and two phases (Liu et al. Reference Liu, Chong, Yang, Verzicco and Lohse2022b). In this work, we have employed a simplified scaling of the form 
$Nu=A Ra^\gamma Pr^m$, rather than the complete GL theory. Since the simplified scaling is an explicit relation for 
$Nu$ as a function of 
$(Ra, Pr)$, we could also derive explicit laws for 
$\varTheta _\varGamma$ and 
$Nu^e/Nu$ as shown in § 2. This would not be possible if adopting the complete GL theory, which provides an implicit relation for 
$(Nu,Re)$ as a function of 
$(Ra,Pr)$ (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001). Even though it is possible to extend the present model with the complete GL theory without conceptual modifications, the simplified scaling is still a valid approximation for the present set-up since the liquid Rayleigh and Prandtl numbers, 
$Ra_l$ and 
$Pr_l$, are similar to those in the gas phase, 
$Ra_g$ and 
$Pr_g$. It is worth emphasizing that the GL theory ceases to be valid when the interface breaks and significant topological changes occur, as shown in Liu et al. (Reference Liu, Chong, Yang, Verzicco and Lohse2022b).
To confirm the validity of the last assumption, we employ the results of the DNS. Figure 8 displays the mean vertical profile of 
$Y_{l}^v$, normalized by 
$\bar {Y}_{l,\varGamma }^v$, for the different flow configurations under investigation. In all cases, we observe a small positive deviation from the interface values, less than 
$1\,\%$ for the highest 
$\varepsilon$, which confirms the validity of approximating 
$Y_{l}^v$ with 
$\bar {Y}_{l,\varGamma }^v$.
Vertical distribution of the vapour mass fraction 
$\langle Y_{l}^v\rangle _x$ (normalized by 
$\bar {Y}_{l,\varGamma }^v$) for 
$Ra=10^6$, 
$10^7$ and 
$10^8$, and 
$\varepsilon =0.05$ (dot-dashed lines) and 
$\varepsilon =0.2$ (solid lines).
Finally, it is worth mentioning that despite the overall model being assessed with two-dimensional simulations, we believe its validity is going to be confirmed also in a three-dimensional configurations without any apparent modification. As discussed in Van Der Poel, Stevens & Lohse (Reference Van Der Poel, Stevens and Lohse2013), the simplified scaling 
$Nu\sim Ra^\gamma$ is valid in both two and three dimensions and none of the three assumptions set restrictions on the dimensionality of the problem.
5. Conclusions
We propose a model for the analytical estimation of the interface temperature 
$\varTheta _{\varGamma }$ and the heat transfer modulation, quantified by the Nusselt number, for an evaporating two-layer RB configuration at a statistically stationary state. The model is based on three assumptions: (i) the OB approximation can be applied to the liquid phase, while the gas thermophysical properties are generic functions of the thermodynamic pressure, local temperature and vapour composition, (ii) the GL theory for thermal convection can be applied to the liquid and gas layers separately and (iii) the vapour content in the gas can be taken as the mean value at the gas–liquid interface. The model provides a quantitative prediction of the ratio 
$Nu^e/Nu$, enabling us to predict the global heat transfer of the evaporating system once the value for the same system in dry conditions, i.e. without phase change, is known. We validate the analytical predictions using DNS in the low-Mach-number regime in a parameter space defined by 
$10^6\leq Ra\leq 10^8$ and 
$0.05\leq \varepsilon \leq 0.20$. Simulations are performed in two settings: (i) assuming the gas density and liquid-diffusion coefficient are the only variable property and (ii) in the general case where all the gas thermophysical properties depend on the state variables. Irrespective of the setting, a very good agreement between the model predictions and the numerical simulations is found by just adopting reasonable values for the scaling exponent of the GL theory, i.e. 
$1/4\leq \gamma \leq 1/3$. Finally, we assess the basic assumptions on which the entire model is built and conclude that they are generally valid unless the interface undergoes large deformation and breakup, which would make the GL theory no longer valid (Liu et al. Reference Liu, Chong, Yang, Verzicco and Lohse2022b). We believe that the proposed model and further extensions may find applications where accurate predictions of 
$\varTheta _\varGamma$ and 
$Nu$ are required and to improve single-phase models of turbulent convection in the presence of evaporation (Schumacher & Pauluis Reference Schumacher and Pauluis2010; Hay & Papalexandris Reference Hay and Papalexandris2020).
Acknowledgements
Dr P. Costa is acknowledged for useful discussions, and Mr W. Gross is thanked for pointing out some typos in an early version of the manuscript.
Funding
N.S., A.D.D. and L.B. acknowledge support from the Swedish Research Council via the multidisciplinary research environment INTERFACE (Pr. 2016-06119). Computer time was provided by the Swedish National Infrastructure for Computing (SNIC) and by the Norwegian Research Infrastructure Services (NRIS, Pr. NN9561K).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation
For completeness, we validate our code using the data from the two-layer RB convection study in Liu et al. (Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a). Here, the OB approximation is assumed for the gas and liquid layers, and phase change is absent. To reproduce this condition, we ‘switch off’ phase change and solve only (3.1), (3.2) and (3.4) setting 
$\varepsilon =0.005$. Moreover, we consider 
$Ra=10^8$, 
$\lambda _{\rho }=3.33$ and Weber number 
$We=5$, while the remaining thermophysical property ratio 
$\lambda _\xi =1$. Simulations are conducted in a two-dimensional domain discretized with 
$N_x=1024$ and 
$N_z=512$. Figure 9 displays the mean vertical temperature profile. The excellent agreement between our simulations and the reference data in Liu et al. (Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a) validates our numerical algorithm for a two-layer RB configuration in the OB limit.
Mean temperature profiles obtained with the numerical code employed in the present study (continuous line) and the results in Liu et al. (Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021a) (circles) for a two-dimensional two-fluid RB flow with 
$Ra=10^8$, 
$\lambda _\rho =3.33$ and 
$We=5$.
Appendix B. Equations of state
We detail the equations of state employed to compute the gas density, specific heat capacity, dynamic viscosity and thermal conductivity of the gas phase. Note that the gas density is a function of temperature, thermodynamic pressure and vapour composition. The dependence on the thermodynamic pressure is usually negligible for the remaining thermophysical properties (Reid et al. Reference Reid, Prausnitz and Poling1987) and, therefore, it is omitted in the present work.
B.1. Gas density
The gas density 
$\hat {\rho }_g$ is evaluated with the equation of state for an ideal gas:
\begin{equation} \hat{\rho}_g = \frac{\hat{p}_{th}\hat{M}_m}{\hat{R}_u\hat{T}_g}, \end{equation}
where 
$\hat {p}_{th}$ is the thermodynamic pressure, 
$\hat {R}_u$ is the ideal gas constant, 
$\hat {T}_g$ is the gas temperature and 
$\hat {M}_m$ is the molar mass, computed using the harmonic average between that of the liquid and that of the gas, i.e. 
$\hat {M}_m=(Y_l^v/\hat {M}_l+(1-Y_l^v)/\hat {M}_g)^{-1}$. By introducing a reference thermodynamic pressure 
$\hat {p}_{th,r}$, a reference temperature difference 
$\widehat {{\rm \Delta} T}=2\varepsilon \hat {T}_r$, a reference molar mass 
$\hat {M}_{g,r}=\hat {M}_{g}$ and using 
$\varTheta _g=(\hat {T}_g-\hat {T}_r)/\widehat {{\rm \Delta} T}$, (B1) can be written in dimensionless form as
\begin{equation} \rho_g = \frac{p_{th}M_m}{1+2\varepsilon\varTheta_g}. \end{equation}
Note that in (B2) the reference density 
$\hat {\rho }_{g,r} = \hat {p}_{th,r}\hat {M}_g/(\hat {R}_u\hat {T}_r)$.
B.2. Specific heat capacity
The specific heat capacity 
$\hat {c}_{pg}$ of the gas phase is computed as a linear combination between that of the vapour, 
$\hat {c}_{p,v}$, supposed to be independent of temperature, and that of the dry gas, 
$\hat {c}_{p,dg}$, as
\begin{equation} \hat{c}_{pg} = \hat{c}_{p,v}Y_l^v + \hat{c}_{pg,d}(1-Y_l^v). \end{equation}
The gas heat capacity in dry conditions, 
$\hat {c}_{pg,d}$, is a function of temperature and it is computed as
\begin{equation} \hat{c}_{pg,d} = \hat{C}_1+\hat{C}_2((\hat{C}_3/\hat{T}_g)/ (\sinh(\hat{C}_3/\hat{T}_g)))^2 + \hat{C}_4((\hat{C}_5/\hat{T}_g)/ (\sinh(\hat{C}_5/\hat{T}_g)))^2,\end{equation}
where 
$\hat {C}_{i=1,5}$ are semi-empirical constants for the dry gas taken equal to 
$\hat {C}_{i=1,5}=[4.13 \times 10^{4}$ J (
$\textrm {kmol}\, \textrm {K})^{-1}$, 
$1.34 \times 10^{4}$ J (
$\textrm {kmol}\,\textrm {K})^{-1}$, 
$3012$ K, 
$1.08 \times 10^{4}$ J (
$\textrm {kmol} \, \textrm {K})^{-1}$, 
$1484$ K], as suggested in Reid et al. (Reference Reid, Prausnitz and Poling1987).
Equation (B3) can be rewritten in dimensionless form using the reference specific heat capacity evaluated at 
$\hat {T}_r$ from (B4):
\begin{equation} c_{pg} = \varPi_{cp}Y_l^v + c_{pg,d}(1-Y_l^v),\end{equation}
where 
$\varPi _{cp}=\hat {c}_{p,v}/\hat {c}_{pg,r}$, in the present work equal to 
$2.0216$. Parameter 
$c_{pg,d}$ is computed from the dimensionless expression of (B4). By using 
$\widehat {{\rm \Delta} T}=2\varepsilon \hat {T}_r$, (B4) can be written in dimensionless form as
\begin{equation} c_{pg,d} = C_1+C_2((C_3/\varTheta_g^*)/(\sinh(C_3/\varTheta_g^*)))^2 + C_4((C_5/\varTheta_g^*)/(\sinh(C_5/\varTheta_g^*)))^2,\end{equation}
where 
$\varTheta _g^*=1+2\varepsilon \varTheta _g$ and 
$C_{i=1,5}= [0.6966,0.2259,8.0175,0.1824,3.9502]$.
B.3. Dynamic viscosity
The dynamic viscosity of the gas phase 
$\hat {\mu }_g$ is computed as a combination between that of the vapour, 
$\hat {\mu }_v$, supposed to be independent of temperature, and that of the dry gas, 
$\hat {\mu }_{dg}$. Differently from the specific heat capacity, the mixture rule is typically nonlinear. In the present work, we employ the Wilke–Lee mixture rule as suggested in Reid et al. (Reference Reid, Prausnitz and Poling1987). This reads as
\begin{equation} \hat{\mu}_g =\frac{Y_{l,m}^v\hat{\mu}_v}{Y_{l,m}^v +(1-Y_{l,m}^v)\phi_{vg}}+\frac{(1-Y_{l,m}^v)\hat{\mu}_{dg}}{Y_{l,m}^v \phi_{gv}+(1-Y_{l,m}^v)}.\end{equation}
Equation (B7) requires to first evaluate 
$Y_{l,m}^v$ as
\begin{equation} Y_{l,m}^v = \frac{Y_l^v}{\hat{M}_l} + \frac{1-Y_l^v}{\hat{M}_g}.\end{equation}
Next, the weighting coefficients 
$\phi _{vg}$ and 
$\phi _{gv}$, which are functions of the molar mass of the vapour and the dry gas, are evaluated as
\begin{equation} \phi_{vg}=\frac{\left(1+\sqrt{\dfrac{\hat{\mu}_v}{\hat{\mu}_{g,d}}} \lambda_M^{{-}1/4}\right)^2}{\sqrt{8(1+\lambda_M)}}, \quad \phi_{gv}=\frac{\left(1+\sqrt{\dfrac{\hat{\mu}_{g,d}}{\hat{\mu}_{v}}} \lambda_M^{1/4}\right)^2}{\sqrt{8(1 +\lambda_M^{{-}1})}}.\end{equation}The gas viscosity can be evaluated with the simplified Sutherland's law:
\begin{equation} \hat{\mu}_{g,d} = \hat{\mu}_{g,r}\left(\frac{\hat{T}_g}{\hat{T}_r} \right)^{2/3}.\end{equation}
Equation (B7) can be rewritten in dimensionless form using the reference viscosity evaluated at 
$\hat {T}_r$ from (B10):
\begin{equation} \mu_g =\frac{Y_{l,m}^v\varPi_\mu}{Y_{l,m}^v+(1-Y_{l,m}^v) \phi_{vg}}+\frac{(1-Y_{l,m}^v)\mu_{g,d}}{Y_{l,m}^v \phi_{gv}+(1-Y_{l,m}^v)},\end{equation}
where 
$\varPi _{\mu }=\hat {\mu }_v/\hat {\mu }_{g,r}=0.3321$. Moreover, by using 
$\widehat {{\rm \Delta} T}=2\varepsilon \hat {T}_r$, (B10) can be written in dimensionless form as
\begin{equation} \mu_{g,d} = (1+2\varepsilon\varTheta_g)^{2/3}.\end{equation}B.4. Thermal conductivity
The thermal conductivity 
$\hat {k}_g$ of the gas phase is computed similarly to the gas viscosity, i.e. using the nonlinear Wilke–Lee mixture rule with suitable modifications:
\begin{equation} \hat{k}_g = \frac{Y_{l,m}^v\hat{k}_v}{Y_{l,m}^v+(1-Y_{l,m}^v) \phi_{vg}}+\frac{(1-Y_{l,m}^v)\hat{k}_{dg}}{Y_{l,m}^v \phi_{gv}+(1-Y_{l,m}^v)},\end{equation}
where 
$Y_{l,m}^v$ is evaluated with (B8). Next, the weighting coefficients 
$\phi _{vg}$ and 
$\phi _{gv}$, which are functions of the molar mass of the vapour and the dry gas, are evaluated using expressions similar to (B9a,b):
\begin{equation} \phi_{vg}=\frac{\left(1+\sqrt{\dfrac{\hat{k}_v}{\hat{k}_{g,d}}} \lambda_M^{{-}1/4}\right)^2}{\sqrt{8(1+\lambda_M)}}, \quad \phi_{gv}=\frac{\left(1+\sqrt{\dfrac{\hat{k}_{g,d}}{\hat{k}_{v}}} \lambda_M^{1/4}\right)^2}{\sqrt{8(1+\lambda_M^{{-}1})}}. \end{equation}The gas thermal conductivity in dry condition can be evaluated with the following expression:
\begin{equation} \hat{k}_{g,d} =\frac{\hat{C}_1\hat{T}_g^{\hat{C}_2}}{1+ \dfrac{\hat{C}_3}{\hat{T}_g}+\dfrac{\hat{C}_4}{\hat{T}_g^2}},\end{equation}
where 
$\hat {C}_{i=1,4}$ are semi-empirical constants for the dry gas taken equal to 
$\hat {C}_{i=1,4}=[3.8889 \times 10^{-4}$ W 
$(\textrm {K}\, \textrm {m})^{-1}$, 
$0.7786$, 
$-0.7716$ K, 
$2121.7\, {\rm K}^{2}$], as suggested in Reid et al. (Reference Reid, Prausnitz and Poling1987). Equation (B13) can be rewritten in dimensionless form using the reference thermal conductivity evaluated at 
$\hat {T}_r$ from (B10):
\begin{equation} k_g =\frac{Y_{l,m}^v\varPi_k}{Y_{l,m}^v+(1-Y_{l,m}^v)\phi_{vg}}+ \frac{(1-Y_{l,m}^v)k_{g,d}}{Y_{l,m}^v\phi_{gv}+(1-Y_{l,m}^v)},\end{equation}
where 
$\varPi _{k}=\hat {k}_v/\hat {k}_{g,r}=0.5824$. Parameter 
$k_{g,d}$ is computed from the dimensionless expression of (B15), obtained defining appropriate reference quantities as for the previous thermophysical properties:
\begin{equation} k_{g,d} =\frac{C_1\varTheta_g^{*,C_2}}{1+\dfrac{C_3}{\varTheta_g^*}+ \dfrac{C_4}{\varTheta_g^{*,2}}},\end{equation}
where 
$C_{i=1,4}= [0.3039,0.7786,-0.0021,+5.6476]$.
Open access
