2.1. Governing equations for DNS

DNS studies were performed to gain insights into the physics of natural convection in the porous medium, to determine the coefficients for the macroscopic model, as well as to obtain the validation data. The governing equations for DNS of natural convection in porous media are the Navier–Stokes equations and the species transport equation. In the flow field, the local species concentration differences are small; hence, the Boussinesq approximation is used to account for the buoyancy force (Herwig Reference Herwig2013). Using Einstein's summation convention, the governing microscopic equations for natural convection in porous media are as follows:

(2.1)\begin{gather}\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0,\end{gather}

(2.2)\begin{gather}\frac{{\partial {u_i}}}{{\partial t}} + \frac{{\partial ({u_i}{u_j})}}{{\partial {x_j}}} ={-} \frac{{\partial p}}{{\partial {x_i}}} + \nu \frac{{{\partial ^2}{u_i}}}{{\partial x_j^2}} + \beta {g_i}(c - c_{0}),\end{gather}

(2.3)\begin{gather}\frac{{\partial c}}{{\partial t}} + \frac{{\partial ({u_i}c)}}{{\partial {x_i}}} = {D_f}\frac{{{\partial ^2}c}}{{\partial x_j^2}},\end{gather}

where $\nu $, ${D_f}$, ${u_i}$, p, ${g_i}$ and c are the kinematic viscosity, the mass diffusivity, the ith component of the velocity vector, the pressure, the ith component of the gravity vector and the species concentration, respectively. The concentration expansion coefficient is defined as $\beta = \beta ({c_0}) ={-} (1/{\rho _0}){(\partial \rho /\partial c)_{{c_0}}}$ (see Herwig & Moschallski, Reference Herwig and Moschallski2009), where $\rho $ is the fluid density.

The Sherwood number $Sh$ is calculated from the DNS as the ratio of the total mass transfer rate $\dot{m}$ (by convection and diffusion) to the mass transfer rate ${\dot{m}_{diff}}$ (by diffusion only) across the lower or upper wall (Baehr & Stephan Reference Baehr and Stephan2006):

(2.4)\begin{equation}Sh = \frac{{\dot{m}}}{{{{\dot{m}}_{diff}}}} = \frac{{\displaystyle\int_w {\overline {\frac{{\partial c}}{{\partial {x_2}}}} \,\textrm{d}A} }}{{\displaystyle\int_w {{{\left. {\overline {\frac{{\partial c}}{{\partial {x_2}}}} } \right|}_{R{a_f} = 0}}\,} \textrm{d}A}},\end{equation}

where the bar denotes the time-averaging operator, while the subscript w denotes either the upper or lower wall surface.

2.2. Macroscopic equations

The macroscopic equations are obtained by averaging the Navier–Stokes equations and the species transport equations (2.1)–(2.3) over each REV (see figure 1). This method of averaging is similar to that used in de Lemos (Reference De Lemos2012); however, de Lemos (Reference De Lemos2012) carried out a time and volume averaging over each respective REV, while we performed only volume averaging. The macroscopic equations read:

(2.5)\begin{gather}\frac{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_{i}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_i}}} = 0,\end{gather}

(2.6)\begin{gather}\frac{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{t} }} + \frac{{\partial ({\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}/\phi )}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_j}}} + \frac{{\partial (\phi {{{\langle ^i}{u_i}{}^i{u_j}\rangle }^i})}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_j}}} ={-} \frac{{\partial (\phi {{\langle p\rangle }^i})}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_i}}} + \nu \frac{{{\partial ^2}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2}} - \phi \beta {g_i}(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} - {c_0}) - \phi {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i},\end{gather}

(2.7)\begin{gather}\frac{{\partial (\phi \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} )}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{t} }} + \frac{{\partial ({\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} )}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_i}}} + \frac{{\partial (\phi {{\langle {}^i{u_i}{}^ic\rangle }^i})}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_j}}} = {D_m}\frac{{{\partial ^2}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} }}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2}},\end{gather}

where the sign $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{ }$ denotes an REV-averaged quantity. The operator ${\langle \;\;\rangle ^i}$ denotes the intrinsic volume averaging in the fluid phase, which is adopted from Whitaker (Reference Whitaker1986). The left superscript i denotes the intrinsic deviation of a volume-averaged quantity, e.g. ${}^i{u_i} = {u_i} - {\langle \; {u_i}\rangle ^i}$. The porosity $\phi $ is defined as $\phi = {V_{void\ space}}/{V_{total}}$, ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i} = \phi {\langle {u_i}\rangle ^i}$ is the volume-averaged velocity, which is often referred to as the superficial velocity, and $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} = {\langle c\rangle ^i}$ is the intrinsic averaged mass concentration. The subscript m denotes an effective property in the volume-averaged equations, e.g. ${D_m}$ is the effective mass diffusivity. Simulations of small domains are needed to determine the value of ${D_m}$ for a specific pore-scale geometry (see Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020).

The terms $\phi {\langle {}^i{u_i}{}^i{u_j}\rangle ^i}$, $\phi {\langle {}^i{u_i}{}^ic\rangle ^i}$ and ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i}$ are the momentum dispersion, mass dispersion and total drag, respectively. The momentum and mass dispersion terms have been neglected in our model due to the underlying assumptions for convection in porous media with low Darcy numbers (see Appendix A1). Since the local Reynolds number $R{e_K} = \left|\boldsymbol{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}} \right|\sqrt K /\nu$ in our simulations is generally smaller than unity (Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020), the Forchheimer term in ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i}$ can also be neglected (Nield & Bejan Reference Nield and Bejan2017). The effects of the macroscopic velocity gradient on ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i}$ can be modelled with a Laplacian term, which was first proposed by Brinkman (Reference Brinkman1949) and then was extensively studied and improved; see Rao, Kuznetsov & Jin (Reference Rao, Kuznetsov and Jin2020), Zaripov, Mardanov & Sharafutdinov (Reference Zaripov, Mardanov and Sharafutdinov2019), Zhao et al. (Reference Zhao, Wang, Li, Zhang and Mahabaleshwar2018), Liu, Patil & Narusawa (Reference Liu, Patil and Narusawa2007), Valdes-Parada, Ochoa-Tapia & Alvarez-Ramirez (Reference Valdes-Parada, Ochoa-Tapia and Alvarez-Ramirez2007), Vafai (Reference Vafai2005), Starov & Zhdanov (Reference Starov and Zhdanov2001) and Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995) as examples. Here, we model the sum of the total drag ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i}$ and the diffusion term $\nu ({\partial ^2}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i}/\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2)$ in (2.6) as

(2.8)\begin{equation}\nu \frac{{{\partial ^2}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2}} + {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{R} _i} = \frac{\nu }{K}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i} + {\nu _m}\frac{{{\partial ^2}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2}},\end{equation}

where K and ${\nu _m}$ are the permeability and effective viscosity of the porous medium. Simulations of small domains are needed to determine their values a priori (see Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) for details of how they were determined). The macroscopic momentum equation (2.6) is hence simplified to

(2.9)\begin{equation}\frac{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{t} }} + \frac{{\partial ({\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}/\phi )}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_j}}} ={-} \frac{{\partial (\phi {{\langle p\rangle }^i})}}{{\partial {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}_i}}} + {\nu _m}\frac{{{\partial ^2}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}}}{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} _j^2}} - \phi \beta {g_i}(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{c} - {c_0}) - \phi \frac{\nu }{K}{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i}.\end{equation}

2.3. Two-length-scale diffusion assumption

Normalizing the governing equations (2.5), (2.9) and (2.7) using the characteristic concentration difference $\Delta c = {c_1} - {c_0}$, velocity ${u_m} = \beta \Delta cgK/\nu $, length $H$ and time ${t_m} = H/{u_m}$, the following dimensionless macroscopic equations are obtained:

(2.10)\begin{gather}\frac{{\partial {{\hat{u}}_i}}}{{\partial {{\hat{x}}_i}}} = 0,\end{gather}

(2.11)\begin{gather}\frac{{\partial {{\hat{u}}_i}}}{{\partial \hat{t}}} + \frac{{\partial ({{\hat{u}}_i}{{\hat{u}}_i}/\phi )}}{{\partial {{\hat{x}}_j}}} ={-} \frac{{\partial (\phi {{\langle \tilde{p}\rangle }^i})}}{{\partial {{\hat{x}}_i}}} + \frac{{{a_\nu }Sc}}{{{\gamma _m}Ra}}\frac{{{\partial ^2}{{\hat{u}}_i}}}{{\partial \hat{x}_j^2}} - \frac{{\phi \textrm{Sc}}}{{{\gamma _m}RaDa}}{z_i}\hat{c} - \frac{{\phi Sc}}{{{\gamma _m}RaDa}}{\hat{u}_i},\end{gather}

(2.12)\begin{gather}\frac{{\partial (\phi \hat{c})}}{{\partial \hat{t}}} + \frac{{\partial ({{\hat{u}}_i}\hat{c})}}{{\partial {{\hat{x}}_i}}} = \frac{1}{{Ra}}\frac{{{\partial ^2}\hat{c}}}{{\partial \hat{x}_j^2}}.\end{gather}

Here $\hat {}$ denotes a dimensionless volume-averaged quantity, $\hat{c}$ is the dimensionless volume-averaged species concentration defined as $\hat{c} = ({\langle c\rangle ^i} - {c_0})/({c_1} - {c_0})$, ${a_\nu } = {\nu _m}/\nu $ is the ratio of the effective viscosity ${\nu _m}$ to the molecular viscosity of the fluid $\nu $, and ${\gamma _m} = {D_m}/{D_f}$ is the ratio of the effective mass diffusivity ${D_m}$ to the mass diffusivity of the fluid ${D_f}$.

The Rayleigh number in (2.11) and (2.12) is defined by using the common definition of this parameter for natural convection in porous media, as in Nield (Reference Nield1994):

(2.13)\begin{equation}Ra \equiv \frac{{R{a_f}Da}}{{{\gamma _m}}} = \frac{{H\beta \Delta cgK}}{{{D_m}\nu }}.\end{equation}

The Schmidt number is defined as

(2.14)\begin{equation}Sc = \frac{\nu }{{{D_f}}}.\end{equation}

The Darcy number is defined as

(2.15)\begin{equation}Da = \frac{K}{{{H^2}}}.\end{equation}

By assuming that ${a_\nu }$ is independent of $Da$ and taking the leading-order terms with respect to $1/Da$ in (2.11), one obtains the well-known DOB equations. However, we reported in our recent DNS study (Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) that the boundary layer thickness is determined by the pore size, which is characterized by $\sqrt K $. In addition, similar profiles for temporally and horizontally averaged quantities are observed when the distance from the wall is normalized with the pore size. Therefore, the Laplacian term $({a_\nu }Sc/{\gamma _m}Ra)({\partial ^2}{\hat{u}_i}/\partial \hat{x}_j^2)$ in (2.11) is expected to scale as $1/K$ and should be of order $1/Da$, exactly as the Darcy term $- (\phi Sc/{\gamma _m}RaDa){\hat{u}_i}$ and the buoyancy force term $- (\phi Sc/{\gamma _m}RaDa){z_i}\hat{c}$. Hence, the Laplacian term in the macroscopic equation cannot be neglected even if the Darcy number is small.

We here propose a model for the effective viscosity ${\nu _m}$ based on a two-length-scale diffusion (TLSD) hypothesis, in which the macroscopic diffusion is determined by the pore size, characterized by $\sqrt K $, and the distance between the lower and upper boundaries H. Our TLSD hypothesis is supported by our recent DNS (Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020), where we showed that natural convection in porous media is determined by these two length scales. The pore size characterizes the boundary layer thickness and the size of proto-plumes, whereas the distance between the two walls determines the size of mega-plumes.

Based on the TLSD assumption stated above, ${a_\nu }$ is modelled as

(2.16)\begin{equation}{a_\nu } = \frac{{{\nu _m}}}{\nu } = \frac{{a_\nu ^\ast }}{{Da}},\end{equation}

where $a_\nu ^\ast $ is a constant assumed to be solely determined by the pore-scale geometry of the porous matrix. Note that the two length scales $\sqrt K $ and H are combined in $Da$. At the upper and lower walls, we imposed constant species concentrations, ${c_1}$ and ${c_0}$, respectively, and the no-slip boundary condition.

It should be noted that only the leading-order terms of $Da$ for diffusion are kept in the macroscopic equations (2.10)–(2.12). As $Da \to 0$, the macroscopic governing equations can be further simplified to

(2.17)\begin{gather}\frac{{\partial {{\hat{u}}_i}}}{{\partial {{\hat{x}}_i}}} = 0,\end{gather}

(2.18)\begin{gather}\frac{{\partial \hat{p}}}{{\partial {{\hat{x}}_i}}} + \hat{c}{z_i} + {\hat{u}_i} = \frac{{a_\nu ^\ast }}{\phi }\frac{{{\partial ^2}{{\hat{u}}_i}}}{{\partial \hat{x}_j^2}},\end{gather}

(2.19)\begin{gather}\frac{{\partial \hat{c}}}{{\partial \hat{t}^*}} + \frac{{\partial ({{\hat{u}}_i}\hat{c})}}{{\partial {{\hat{x}}_i}}} = \frac{1}{{Ra}}\frac{{{\partial ^2}\hat{c}}}{{\partial \hat{x}_j^2}},\end{gather}

where $\hat{p} = RaDa{\langle \tilde{p}\rangle ^i}/{\gamma _m}Sc$ is the normalized pressure. The dimensionless time is modified to be $\hat{t}^* = \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{t} /\phi $. These macroscopic equations (2.17)–(2.19) become the DOB equations if $a_\nu ^\ast $ is set to zero. When the DOB equations are solved, only the velocity component in the wall-normal direction is set to zero at the upper and lower walls. This boundary condition was also used in other DOB simulations; see Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) and Wen et al. (Reference Wen, Corson and Chini2015) as examples. In this paper, the macroscopic simulations were carried out by solving equations (2.10)–(2.12), so that the effect of the Darcy number can be assessed. The Sherwood number for the macroscopic model simulations is defined using the same definition as for the DNS, which is given in (2.4).

2.4. Numerical method

For the simulations, a finite-volume method (FVM) was utilized. The solvers were developed by using the open-source code package OpenFoam 6. The spatial discretization was implemented by a second-order central-difference scheme. For time derivatives, the second-order implicit backward method was used. For the correction and coupling of the pressure and velocity fields, the pressure-implicit scheme with splitting of operators (PISO) algorithm was used (Versteeg & Malalasekera Reference Versteeg and Malalasekera2007). A stabilized preconditioned (bi)conjugate gradient solver was utilized to solve the pressure field and the momentum and species concentration equations. We have performed the code validation for our DNS solver extensively in our previous studies (Jin et al. Reference Jin, Uth, Kuznetsov and Herwig2015; Uth et al. Reference Uth, Jin, Kuznetsov and Herwig2016; Jin & Kuznetsov Reference Jin and Kuznetsov2017; Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020).

3.1. Description of the test cases

We continued selected DNS cases of Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) to improve the statistics and thus to allow a more thorough validation of our hypothesis. In addition, we also computed these cases by solving the macroscopic equations (2.10)–(2.12) with our two-length-scale diffusion model. The Rayleigh numbers $Ra$ are up to $20\;000$ and the Schmidt numbers $Sc$ are $1$ and $250$. The ranges of geometrical parameters of the studied test cases are given in table 1. For both DNS and macroscopic simulation cases ${u_i} = 0$ and $c = ({c_1} - {c_0})/2$ were used as initial conditions.

Table 1. Ranges of geometrical parameters for the studied test cases.

To obtain representative statistical results, the time averaging, denoted by the bar over, of the respective variable was performed after the flow and mass concentration fields reached a statistically steady state. As an example, the time evolutions of the instantaneous Sherwood number for the DNS case with $H/s = 100$, $s/d = 1.5$, $Ra = 20\;000$ and $Sc = 250$ are shown in figure 2. The time averaging of the Sherwood number has been started after the time marked by the red dashed line. At least 200 dimensionless time units $H/{u_m}$ are calculated to obtain the statistical results.

Figure 2. The time evolution of the instantaneous Sherwood number for the DNS case with $H/s = 50$, $s/d = 1.5$, $Ra = 20\;000$ and $Sc = 250$. The dashed red line marks the time at which the time averaging is started; and $\hat{t} = t{u_m}/H$ is the dimensionless time.

3.2. Determination of the model coefficient

The coefficient $a_\nu ^\ast $ for a specific pore-scale geometry cannot be computed a priori with simulations of small domains (as for the other model parameters). Here, we empirically determine $a_\nu ^\ast $ by simulating natural convection within the specific pore-scale geometry with fixed values of $H/s$, $Sc$ and $Ra$. Since we only keep the leading-order terms of the order $Da$ in our model equations, a test case with a sufficiently small Darcy number should be used to ensure that the higher-order terms of Da are negligible. In particular, we performed a parametric study for $a_\nu ^\ast (\phi )$ while keeping $H/s = 20,\;Sc = 250$ and $Ra = 20\;000$ fixed. These parameter values were selected because the Sherwood number from DNS marginally changes as $H/s$ is increased (i.e. as the pore size is decreased); hence the effect of higher-order terms of Da on Sh can be safely neglected. The value of $a_\nu ^\ast $ is selected for each considered porosity value, so that the Sherwood number from the macroscopic simulation matches the DNS results. Figure 3 shows the dependence of $a_\nu ^\ast $ on the porosity $\phi $.

Figure 3. Dependence of the model coefficient $a_\nu ^\ast $ on the porosity $\phi $, for porous matrices composed of aligned and staggered square obstacles.

We expect that $a_\nu ^\ast $ is a geometrical parameter that is independent of $Sc$, $Ra$ and Da. This will be examined later in § 5. The value of $a_\nu ^\ast $ only mildly changes when the porous matrix is switched from aligned squares to staggered squares. According to our DNS results, $a_\nu ^\ast $ can be well correlated by

(3.1)\begin{equation}a_v^\ast (\phi ) = 7.5 \times {10^{ - 5}}{\phi ^9}.\end{equation}

This correlation has reasonable accuracy for $0.28 < \phi < 0.95$. However, the variations of pore-scale geometries used in this study are limited. In particular, the flow structures in randomly packed porous matrices may be distinctly different from those in regularly packed porous matrices (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020a). Studies with more pore-scale geometries are needed to test the generality of (3.1).

3.3. Mesh and time-step independence studies

The mesh and time-step independence studies for the DNS cases have already been performed in our previous work (see Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020). Here we focus on the influence of the mesh and time step on the macroscopic simulation results (solution of (2.10)–(2.12)). The numerical results for the Sherwood number are shown in table 2. At least 200 dimensionless time units $H/{u_m}$ were calculated to obtain the statistical results.

Table 2. Influence of mesh and time step on the Sherwood number $Sh$. The test case with $H/s = 50$, $s/d = 1.5$, $Ra = 20\;000$ and $Sc = 250$ is used in the parametric study. The cases c, d, e and f are considered to be mesh- and time-step-independent. The mesh resolution and maximum Courant number of the case f (in italic) are used in all cases of macroscopic simulation.

The results of the resolution study show that the Sherwood number is under predicted if the mesh resolution is too low (see table 2 cases a and b) or the maximum Courant number is too high (see table 2 case g). According to the mesh/time-step independence study, $C{o_{max}} = 0.8$ and mesh resolution 2000×3200 (case f) were used for all cases of macroscopic simulation. The numerical error of $Sh$ in the macroscopic simulations is estimated to be 2.8 %, which is the maximum variation of $Sh$ in the cases c, d, e and f. All simulations were performed on the clusters of the HLRN (North-German Supercomputing Alliance), using 2× Intel Cascade Lake Platinum 9242 CPUs (CLX-AP) with 96 cores per node. The DNS cases use up to $7.2 \times {10^7}$ mesh cells, which requires a parallel computing time of 1200 hours using 384 processors. The macroscopic simulation cases use up to $6.4 \times {10^6}$ mesh cells.

4.1. Budget of the macroscopic kinetic energy

The budget for the time- and line-averaged macroscopic kinetic energy ${\langle \bar{K}\rangle ^{x1}} = {\textstyle{1 \over 2}}{\langle \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i^2\rangle ^{x1}}$ was calculated from the DNS for $Ra = 20\;000,\;s/d = 1.5$, $s/d = 1.25$, $Sc = 250$ and $Da$ in the range $3.5 \times {10^{ - 7}}\textrm{ to }3.5 \times {10^{ - 5}}$. By averaging the momentum equation (2.2) over REVs, taking the dot product with the superficial velocity ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u} _i} = \phi {\langle {u_i}\rangle ^i}$, and then averaging in time and in the horizontal direction ${x_1}$, we obtained the following equation for ${\langle \bar{K}\rangle ^{x1}}$:

(4.1) \begin{align}&- {\left\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\left\langle {\frac{{\partial ({u_i}{u_j})}}{{\partial {x_j}}}} \right\rangle }^i}} \right\rangle^{x1}} - {\left\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\left\langle {\frac{{\partial p}}{{\partial {x_i}}}} \right\rangle }^i}} \right\rangle^{x1}} + {\left\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\left\langle {\nu \frac{{{\partial^2}{u_i}}}{{\partial x_j^2}}} \right\rangle }^i}}\right\rangle ^{x1}}\notag\\ &\quad + {\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\langle \beta {g_i}(c - {c_0})\rangle }^i}}\rangle^{x1}} = 0.\end{align}

Equation (4.1) shows that the budget for ${\langle \bar{K}\rangle ^{x1}}$ includes:

• • the production by the buoyancy force, ${K_{buoy}} = {\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\langle \beta {g_i}(c - {c_0})\rangle }^i}}\rangle^{x1}}$;

• • the loss due to viscous dissipation, ${K_{diff}} = {\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\langle \nu ({\partial ^2}{u_i}/\partial x_j^2)\rangle }^i}}\rangle^{x1}}$;

• • the loss due to pressure gradient, ${K_{pres}} ={-} {\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\langle \partial p/\partial {x_i}\rangle }^i}}\rangle^{x1}}$;

• • the transport due to convection, ${K_{conv}} ={-} {\langle \overline {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i}\phi {{\langle \partial ({u_i}{u_j})/\partial {x_j}\rangle }^i}}\rangle^{x1}}$.

In the DOB equations, the Darcy term (Darcy drag) is the only source of losses of macroscopic kinetic energy. The Darcy losses read

(4.2)\begin{equation}{K_{Darcy}} ={-} \phi (\nu /K){\langle \overline {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{u}_i^2}\rangle^{x1}}.\end{equation}

The budget of ${\langle \bar{K}\rangle ^{x1}}$ is studied using the test case with $s/d = 1.5\ (\phi = 0.56)$, $H/s = 20$, $Ra = 20\;000$, $Da = 8.8 \times {10^{ - 6}}$ and $Sc = 250$. Figure 4 shows the distribution of ${K_{buoy}}$, ${K_{diff}}$, ${K_{pres}}$ and ${K_{conv}}$ in the wall-normal direction. They are normalized with the characteristic kinetic energy ${K_{mean}} = {\textstyle{1 \over 2}}\; u_m^2$ or ${K_{buoy}}$. The distance from the lower wall is normalized with the pore size s. It is evident that more macroscopic kinetic energy is produced by the buoyancy force in the central region than in the region close to the wall. The transport of ${\langle \bar{K}\rangle ^{x1}}$ due to convection is much smaller than ${K_{buoy}}$, so it can be neglected. Here $- {K_{diff}}$ and $- {K_{pres}}$ are the losses of the macroscopic kinetic energy. Both $- {K_{diff}}$ and $- {K_{pres}}$ increase with increasing distance from the wall ${x_2}/s$. The loss of ${\langle \bar{K}\rangle ^{x1}}$ in the region close to the wall is mainly due to the pressure gradient $- {K_{pres}}$.

Figure 4. Distribution of the budget of the macroscopic kinetic energy ${\langle \bar{K}\rangle ^{x1}}$ in the wall-normal direction. Here $s/d = 1.5\ (\phi = 0.56)$, $H/s = 20\ (Da = 8.8 \times {10^{ - 6}})$, $Ra = 20\;000$ and $Sc = 250$.

Figure 5 shows the loss of the macroscopic kinetic energy due to the Darcy drag ${K_{Darcy}}$ (assuming that the superficial velocity calculated from the macroscopic simulation is identical to the DNS solution). The drag ${K_{Darcy}}$ is normalized by ${K_{mean}}$ or ${K_{buoy}}$. It can be seen that ${K_{Darcy}}$ is close to ${K_{buoy}}$ in the region away from the wall $({x_2}/s \gg 0)$. However, ${K_{Darcy}}/{K_{buoy}}$ is smaller than 0.85 in the first three REVs adjacent to the wall $({x_2}/s < 3)$. The DNS results confirm that the Darcy term, which accounts for the losses due to the Darcy drag, cannot account for all the losses of the macroscopic kinetic energy.

Figure 5. Distribution of the loss of the macroscopic kinetic energy ${\langle \bar{K}\rangle ^{x1}}$ due to the Darcy drag. Here $s/d = 1.5\ (\phi = 0.56)$, $H/s = 20\ (Da = 8.8 \times {10^{ - 6}})$, $Ra = 20\;000$ and $Sc = 250$.

The question arises of whether the difference between ${K_{Darcy}}$ and ${K_{buoy}}$ shown in figure 5 is because the Darcy number in the DNS case is not small enough. To answer this question, ${K_{pres}}/{K_{buoy}}$, ${K_{diff}}/{K_{buoy}}$, ${K_{conv}}/{K_{buoy}}$ and ${K_{darcy}}/{K_{buoy}}$ in the first REV cell next to the bottom wall for different Darcy numbers are compared in figure 6. It is evident from this figure that all of these quantities stay almost constant as the Darcy number is decreased from $3.5 \times {10^{ - 5}}$ to $3.5 \times {10^{ - 7}}$, suggesting that the Darcy numbers in our DNS cases are small enough for the presented analysis. The Darcy number has a noticeable effect as it is increased to ${\sim} 3 \times {10^{ - 5}}$. In this case, ${K_{pres}}/{K_{buoy}}$ and ${K_{Darcy}}/{K_{buoy}}$ become smaller, ${K_{diff}}/{K_{buoy}}$ becomes larger, whereas ${K_{conv}}/{K_{buoy}}$ is still negligibly small. We speculate that higher ${K_{diff}}$ at very large Darcy numbers is due to the mass dispersion, which is neglected in our macroscopic model (convection with very large Darcy numbers is out of the scope of this study).

Figure 6. Plots of ${K_{pres}}/{K_{buoy}}$ (a), ${K_{diff}}/{K_{buoy}}$ (b), ${K_{conv}}/{K_{buoy}}$(c) and ${K_{Darcy}}/{K_{buoy}}$ (d) in the first REV cell next to the bottom wall versus the Darcy number. Here $s/d = 1.5\ (\phi = 0.56)$ with $H/s = 10,\;20,\;50,\;100$ and $s/d = 1.25\ (\phi = 0.36)$ with $H/s = 10,\;20,\;50$, $Ra = 20\;000$ and $Sc = 250$.

Our budget analysis shows that, in the near-wall region, there is a difference between the loss due to the Darcy drag ${K_{Darcy}}$ and the overall loss, which is identical to $- {K_{buoy}}$. Since the transport of ${\langle \bar{K}\rangle ^{x1}}$ is negligibly small, this suggests that another source for the loss of ${\langle \bar{K}\rangle ^{x1}}$ should be considered in the macroscopic equations.

4.2. Sh–Da dependence

According to our hypothesis, the macroscopic diffusion, the Darcy drag and the buoyancy force are of the same order with respect to the Darcy number, so the macroscopic diffusion cannot be neglected even if the Darcy number is small. To examine our hypothesis, we investigated the relationship between the Sherwood number and the Darcy number. We varied $Da$ in the range $3.5 \times {10^{ - 7}}\textrm{ to }3.5 \times {10^{ - 5}}$ for $s/d = 1.5\ (\phi = 0.56)$ and $2.8 \times {10^{ - 7}}\textrm{ to }7 \times {10^{ - 6}}$ for $s/d = 1.25\ (\phi = 0.36)$; the corresponding $H/s$ ratios are in the range 10–100. If our hypothesis were true, the Sherwood number should gradually become independent of $Da$ and should not approach the DOB solution.

The DNS results shown in figure 7 generally support our assumption, i.e. the values of $Sh$ for $Sc = 250$ depend only weakly on $Da,$ when $Da$ is small enough. The values of $Sh$ are also different from the DOB solution. As shown in figure 7(a), the same trend is found for $s/d = 1.25\ (\phi = 0.36)$ and $Sc = 1$, where $Sh$ depends weakly on $Da$. The only exception is the case for $s/d = 1.5\ (\phi = 0.56)$ and $Sc = 1$, where $Sh$ still increases with decreasing Da (but it is still far away from the DOB result). Test cases with even smaller Darcy numbers could be computed to probe the Da dependence more thoroughly. However, the calculation of these cases would be extremely expensive and hence out of the scope of this study.

Figure 7. The $Sh(Da)$ dependence for the DNS and DOB cases. Here $s/d = 1.5\ (\phi = 0.56)$ with $H/s = 10,\;20,\;50,\;100$ and $s/d = 1.25\ (\phi = 0.36)$ with $H/s = 10,\;20,\;50$ and $Ra = 20\;000$, for (a) $Sc = 1$ and (b) $Sc = 250$.

It should be noted that the Darcy numbers for real applications are much smaller than the values used in the DNS cases. However, since our DNS results for the Sherwood number are approximately $Da$-independent, we expect that it is possible to predict the Sherwood numbers using DNS with relatively higher (computationally affordable) Darcy numbers.

In a recent DNS study, Liu et al. (Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020a) proposed the following correlation for estimating the Nusselt number (equivalent to $Sh$ in this study):

(4.3)\begin{equation}Sh \approx c{\kern 1pt} \phi {\left( {\frac{H}{l}} \right)^4}S{c^2}Re_{rms}^2Ra_f^{ - 1} + 1,\end{equation}

where c is an undetermined constant according to the work of Grossmann & Lohse (Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2004), $\; l$ is the minimum spacing between the obstacles, $R{e_{rms}} = {U_{rms}}l/\nu $ is the Reynolds number based on the volume-averaged r.m.s. velocity magnitude, and $R{a_f} = {H^3}\beta \Delta cg/\nu {D_f}$ is the Rayleigh number defined for the free fluid flow. If we set the value of c to $1250$ and determine ${U_{rms}}$ from our DNS results, the results of (4.3) are in good agreement with our DNS results for different values of $\phi $ and $Sc$ (see figure 8). It should be noted that (4.3) is proposed based on the flow condition that viscosity dominates; hence intense kinetic energy dissipation takes place within the bulk domain and turbulence is suppressed in the pore canals. For the volume- and time-averaged kinetic energy dissipation rate ${\langle {\epsilon _u}\rangle ^{v,t}}$, the following proportionality is valid: ${\langle {\epsilon _u}\rangle ^{v,t}}\sim \phi \nu U_{rms}^2/{l^2}$. This corresponds to the ∞ regime of classical Rayleigh–Bénard convection (without porous media) introduced by Grossmann & Lohse (Reference Grossmann and Lohse2001) for large $Sc$ and small $R{a_f}$. A good agreement between predictions obtained using (4.3) and our DNS results indicates the significance of macroscopic diffusion in momentum transport.

Figure 8. Sherwood number versus the Rayleigh number for $Ra$ in the range $500 - 20\;000$ compared to the correlation proposed by Liu et al. (Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020a), with (a) $Sc = 1$ and (b) $Sc = 250$.

5.1. Sherwood number

Figure 9 shows the relationship between the Sherwood number and the Rayleigh number when $H/s$ is 20 and Rayleigh numbers are up to $20\;000$. The results of our macroscopic model are compared with the correlation obtained from the DNS results (see Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) as well as the DOB results. In the DNS, $Sh$ depends on $s/d$ or $\phi $ for $Ra > 2000$. The value of $Sh$ increases as $s/d$ or $\phi $ decreases, while the difference becomes larger as $Ra$ increases. In the large-Rayleigh-number regime $(Ra \ge 5000)$, the $Sh = f(Ra)$ scaling changes from a linear scaling $Sh\sim Ra$ for $s/d = 1.25$ $(\phi = 0.36)$ to a nonlinear scaling $Sh\sim R{a^{0.8}}$ for $s/d = 1.5$ $(\phi = 0.56)$, see figure 9(a). These characteristics are not captured in DOB simulations but are well reproduced in our macroscopic model simulations.

Figure 9. Sherwood number versus the Rayleigh number with $Ra$ in the range $500 - 20\;000$ and $H/s = 20$ for three values of the Darcy number: $Da = 8.8 \times {10^{ - 6}}\ (s/d = 1.5)$, $Da = 5.4 \times {10^{ - 6}}\ (s/d = 1.4)$ and $Da = 1.8 \times {10^{ - 6}}\ (s/d = 1.25)$, with (a) $Sc = 1$ and (b) $Sc = 250$.

For the current small $H/s$ (or $Da$) value, both DNS and model results show that the Sherwood number increases as the Schmidt number is increased from 1 to 250. Similar to $Sc = 1$, the scaling for $Sc = 250$ also changes from linear $(Sh\sim Ra)$ for $s/d = 1.25$ $(\phi = 0.36)$ to nonlinear $(Sh\sim R{a^{0.8}})$ for $s/d = 1.5\ (\phi = 0.56)$. This behaviour is reproduced in our macroscopic model solution, whereas the effect of the Schmidt number is not accounted for in the DOB equations.

We neglected the high-order terms with respect to $Da$ when we proposed the TLSD model. However, since the leading-order term with respect to $Da$ is kept in the momentum equation (2.11), the effect of $Da$ on natural convection can still be accounted for when its value is small enough. The relationship between the Sherwood number and the Darcy number for $Ra = 20\;000$, $Sc = 1$ or 250, and $s/d = 1.5$ or 1.25 ($\phi = 0.56$ or 0.36) is shown in figure 10. The macroscopic model results are compared with the DNS results as well as with the DOB results. Recall that $a_\nu ^\ast $ is only related to the pore-scale geometry and is independent of the flow conditions and the Darcy number. The macroscopic model simulations are in good agreement with the DNS for $Da \le 2 \times {10^{ - 6}}$ and for different Schmidt numbers. By contrast, the DOB results are independent of the Darcy and Schmidt numbers. The DOB simulations overpredict the Sherwood number for $s/d = 1.5\ (\phi = 0.56)$ and underpredict Sh for $s/d = 1.25\ (\phi = 0.36)$.

Figure 10. Sherwood number versus the Darcy number for $s/d = 1.5\ (\phi = 0.56)$ with $H/s = 10,\;20,\;50,\;100$ and $s/d = 1.25\ (\phi = 0.36)$ with $H/s = 10,\;20,\;50$, $Ra = 20\;000$, for (a) $Sc = 1$ and (b) $Sc = 250$.

5.2. Species concentration and velocity statistics

The vertical profiles of the temporally and horizontally averaged macroscopic quantities (time- and line-averaged species concentration ${\langle \overline {\hat{c}} \rangle ^{x1}}$, species concentration fluctuation ${\langle {\hat{c}^{rms}}\rangle ^{x1}}$, and velocity fluctuations ${\langle \hat{u}_1^{rms}\rangle ^{x1}}$ and ${\langle \hat{u}_2^{rms}\rangle ^{x1}}$) for $s/d = 1.5$ ($\phi = 0.56$) and $Sc = 1$ are shown in figure 11. The Rayleigh numbers are $5000\;\textrm{and}\;20\;000$. The results of our macroscopic TLSD model are compared with the DNS results as well as the DOB results. It is evident that our macroscopic TLSD modelling results for ${\langle \overline {\hat{c}} \rangle ^{x1}}$, ${\langle \hat{u}_1^{rms}\rangle ^{x1}}$ and ${\langle \hat{u}_2^{rms}\rangle ^{x1}}$ are more accurate than the DOB results in the first REV next to the wall. The DNS results show that all statistical results can be well scaled by the pore size s and that the influence of the bounding walls is limited to within the first three REVs next to the wall (Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020). Thus, the boundary layer thickness is determined by the pore size s, instead of the Rayleigh number, as suggested in Huppert & Neufeld (Reference Huppert and Neufeld2014). These features are all well captured in the macroscopic simulation.

Figure 11. The vertical profiles of the temporally and horizontally averaged macroscopic quantities for $s/d = 1.5\ (\phi = 0.56)$, $H/s = 20\ (Da = 8.8 \times {10^{ - 6}})$ and $Sc = 1$. The Rayleigh number $Ra$ is varied. The distance from the wall is normalized by the pore size s. (a) Time- and line-averaged species concentration ${\langle \overline {\hat{c}} \rangle ^{x1}}$; (b) r.m.s. of the species concentration fluctuation ${\langle {\hat{c}^{rms}}\rangle ^{x1}}$; (c) streamwise velocity fluctuation ${\langle \hat{u}_1^{rms}\rangle ^{x1}}$; and (d) wall-normal velocity fluctuation ${\langle \hat{u}_2^{rms}\rangle ^{x1}}$.

The same statistical quantities are shown in figure 12 for $Sc = 250$. It can be seen that these statistical quantities are only marginally changed when a much higher Schmidt number is used in the simulation. Similar to the results for $Sc = 1$, the macroscopic modelling results are also in good agreement with the DNS results. The macroscopic TLSD model simulation predicts higher mass concentration ${\langle \overline {\hat{c}} \rangle ^{x1}}$ in the first REV next to the wall and higher transverse velocity fluctuation ${\langle \hat{u}_2^{rms}\rangle ^{x1}}$. One possible reason for this discrepancy is that the neglected high-order terms with respect to $Da$ may lead to modelling errors in the boundary layer region. The TLSD model accuracy can be further improved by decomposing the flow domain into a boundary layer region and a central region, so the modelling in the boundary layer region can be improved. However, this would make the model more complicated and difficult to apply. This modelling approach is not adopted in our study to achieve a compromise between the accuracy and simplicity of the macroscopic model.

Figure 12. The vertical profiles of the temporally and horizontally averaged macroscopic quantities for $s/d = 1.5\ (\phi = 0.56)$, $H/s = 20\ (Da = 8.8 \times {10^{ - 6}})$ and $Sc = 250$. The Rayleigh number $Ra$ is varied. The distance from the wall is normalized by the pore size s. (a) Time- and line-averaged species concentration ${\langle \overline {\hat{c}} \rangle ^{x1}}$; (b) r.m.s. of the species concentration fluctuation${\langle {\hat{c}^{rms}}\rangle ^{x1}}$; (c) streamwise velocity fluctuation ${\langle \hat{u}_1^{rms}\rangle ^{x1}}$; and (d) wall-normal velocity fluctuation ${\langle \hat{u}_2^{rms}\rangle ^{x1}}$.

5.3. Transient macroscopic fields

To validate the results of our macroscopic TLSD model, we first compare the transient flow fields obtained from macroscopic simulations of (2.10)–(2.12) with those obtained from the DNS results discussed in the previous section and the DOB simulations reported in Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020). For this purpose, the velocity field and the species concentration obtained with the macroscopic simulations and the DNS were volume-averaged (over each REV).

The distribution of the instantaneous $R{e_K} = (|\boldsymbol{u} |{K^{1/2}})/\nu $ for $Ra = 20\;000$, $s/d = 1.5\ (\phi = 0.56)$, $H/s = 100$ and $Sc = 250$ is shown in figure 13. The macroscopic TLSD solution (figure 13c) is qualitatively similar to the DNS solution (figure 13a) and DOB solution (figure 13b). Both the DNS solution and macroscopic solutions indicate that the local Reynolds number is $R{e_K} < 4 \times {10^{ - 3}}$. This shows that the studied parameter range is well in the Darcy regime $(R{e_K} < 1)$, hence the Forchheimer term in the momentum equation can be safely neglected. Despite the laminar flow in the pore scale, the macroscopic flow field is transient and chaotic. However, the strong spatial variation of the velocity field obtained from the DNS is captured neither in TLSD nor in the DOB simulations.

Figure 13. Instantaneous volume-averaged Reynolds number $R{e_K}$, $H/s = 100\ (Da = 3.5 \times {10^{ - 7}})$, $s/d = 1.5\ (\phi = 0.56)$, $Ra = 20\;000$ and $Sc = 250$: (a) DNS, (b) DOB and (c) TLSD.

Snapshots of the instantaneous species concentrations for $H/s = 100$, $s/d = 1.5\ (\phi = 0.56)$, $Ra = 20\;000$ and $Sc = 250$ are shown in figure 14. The DNS solution (figure 14a), TLSD solution (figure 14c) and the DOB solution (figure 14b) all exhibit large mega-plume structures in the internal region and small proto-plumes in the boundary layers. They occur due to the rising of a fluid with low species concentration and the sinking of a fluid with high species concentration, forming the instabilities in the boundary layer region (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Kränzien & Jin Reference Kränzien and Jin2019).

Figure 14. Instantaneous volume-averaged species concentration $\hat{c}$, $H/s = 100\ (Da = 3.5 \times {10^{ - 7}})$, $s/d = 1.5\ (\phi = 0.56)$, $Ra = 20\;000$ and $Sc = 250$: (a) DNS, (b) DOB and (c) TLSD.

While the macroscopic TLSD model and DOB solution exhibit relatively regular mega-plumes, in the DNS solution the mega-plumes are more irregular and chaotic. A possible reason is that the Darcy number in our simulation is still not small enough, while the TLSD model is proposed for problems with low Darcy numbers. The transient flow field from macroscopic simulation converges slower than the Sherwood number and other statistical results with decreasing $Da$.

The DNS study reported in Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) shows that the number of mega-plumes increases with the decrease of $Da$. Figure 15 shows the time-averaged fast Fourier transform (FFT) of the dimensionless mass concentration $\hat{c}$ along the centreline at ${x_2} = H/2$. The peak wavenumber calculated from the TLSD simulation is still higher than the DNS result, but it is lower than the DOB result. Figure 16 shows that the peak wavenumber from DNS approaches the TLSD or DOB results as the Darcy number approaches 0. However, DNS of natural convection with smaller Darcy numbers are still needed to confirm that the peak wavenumber from DNS will not exceed the TLSD results and approach the DOB results.

Figure 15. Average spectra of the dimensionless mass concentration, $\hat{c}$, of the DNS, DOB and TLSD results at mid-height ${x_2} = H/2$, $H/s = 100\ (Da = 3.5 \times {10^{ - 7}})$, $s/d = 1.5\ (\phi = 0.56)$, $Ra = 20\;000$ and $Sc = 250$.

Figure 16. Peak wavenumber k for the mega-plumes of the DNS, DOB and TLSD results for different Darcy numbers with $s/d = 1.5\ (\phi = 0.56)$, $Ra = 20\;000$ and $Sc = 250$.

The three-dimensional DOB simulations by Pirozzoli et al. (Reference Pirozzoli, De Paoli, Zonta and Soldati2021) revealed the supercells at the boundary, which are the footprints of mega-plumes dominating the interior part of the flow. They suggest that these supercells might lead to the nonlinear scaling of $Sh(Ra)$ in the ultimate regime of high Ra numbers. Future work needs to investigate whether these supercells will be also captured by three-dimensional TLSD simulations. Elucidating how macroscopic diffusion affects the plume structures is also a subject of future investigation.